Rotation of a fibre in simple shear flow of a dilute polymer solution

Abstract

The motion of a freely rotating prolate spheroid in a simple shear flow of a dilute polymeric solution is examined in the limit of large particle aspect ratio, $\kappa$. A regular perturbation expansion in the polymer concentration, $c$, a generalized reciprocal theorem, and slender body theory to represent the velocity field of a Newtonian fluid around the spheroid are used to obtain the $O(c)$ correction to the particle's orientational dynamics. The resulting dynamical system predicts a range of orientational behaviours qualitatively dependent upon $c\, De$ ($De$ is the imposed shear rate times the polymer relaxation time) and $\kappa$ and quantitatively on $c$. At a small but finite $c\, De$, the particle spirals towards a limit cycle near the vorticity axis for all initial conditions. Upon increasing $\kappa$, the limit cycle becomes smaller. Thus, ultimately the particle undergoes a periodic motion around and at a small angle from the vorticity axis. At moderate $c\, De$, a particle starting near the flow–gradient plane departs it monotonically instead of spirally, as this plane (a limit cycle at smaller $c\, De$) obtains a saddle and an unstable node. The former is close to the flow direction. Upon further increasing $c\, De$, the saddle node changes to a stable node. Therefore, depending upon the initial condition, a particle may either approach a periodic orbit near the vorticity axis or obtain a stable orientation near the flow direction. Upon further increasing $c\, De$, the limit cycle near the vorticity axis vanishes, and the particle aligns with the flow direction for all starting orientations.

1. Introduction

A particle-filled viscoelastic polymeric fluid undergoes simple shear flow in many industrial applications such as fibre spinning (Nakajima, Kajiwara & McIntyre 1994; Breitenbach 2002; Huang et al. 2003; Chae & Kumar 2008) and roll-to-roll manufacturing of high aspect ratio, low resistance films for flexible and transparent electronics (Yin et al. 2010; Mutiso et al. 2013). Simple shear flow occurs in the spinneret during fibre spinning and in the patterning channel during roll-to-roll manufacturing. The suspending viscoelastic fluid may include large aspect ratio particles to impart strength to the final product in fibre spinning or provide a desired anisotropy to the low resistance film. In the simple shear flow of an inertialess Newtonian fluid, a fibre/slender particle undergoes an initial condition-dependent periodic motion in orientational trajectories termed Jeffery orbits (Jeffery 1922) as shown in figure 1 for particles with aspect ratio, $\kappa =10$ and 50. However, the interaction between the polymers in a viscoelastic fluid and the fibre breaks this degenerate periodic behaviour. Previous experiments (Gauthier, Goldsmith & Mason 1971; Bartram, Goldsmith & Mason 1975; Johnson, Salem & Fuller 1990; Stover & Cohen 1990; Iso, Koch & Cohen 1996b; Gunes et al. 2008) indicate that depending on the $\kappa$, imposed shear rate and the properties of the viscoelastic fluid such as polymer concentration and relaxation time, a particle may exhibit various orientation dynamics. A slender particle, i.e. one with a large $\kappa$, may either spiral or monotonically drift towards the vorticity axis, align near the flow direction or settle somewhere within the flow–vorticity plane. Therefore, careful design and choice of flow parameters during the simple shear regime are essential for obtaining a final product with desired particle orientation and material strength. Theoretical studies are useful due to the many parameters required for characterizing a viscoelastic fluid. Polymers lead to new features in a viscoelastic fluid flow such as shear thinning or a finite first and second normal stress difference as compared with a Newtonian fluid flow. Leal (1975) predicts that a slender particle in a slow flow will spiral towards the vorticity axis due to the second normal stress difference in the fluid. Whereas, operating in a double limit of small polymer concentration and large Deborah number, $De$, (the product of the imposed shear rate and the polymer relaxation time) Harlen & Koch (1993) also predict the spiralling of the particle towards the vorticity axis, but due to first normal stress difference in the fluid. Neither of these theories captures any other orientation behaviour observed experimentally. In this paper, using a regular perturbation expansion in polymer concentration, $c$, we develop a slender body theory that spans a range of $De$. It encapsulates the $O(c)$ effect of particle–polymer interaction and qualitatively describes the rich orientation dynamics seen in previous experiments.

Figure 1. Jeffery orbits or orientation trajectories in a simple shear flow of Newtonian fluid of a prolate spheroidal particle with aspect ratio, $\kappa =$ (a) 10, and (b) 50.

Simple shear flow is ubiquitous in industrial applications as the near wall flow can always be considered locally simple shear. Many scenarios include laminar flow between two parallel walls; in these cases, the local flow is simple shear everywhere. In many industrial scenarios, particle concentration in the suspension is dilute, and due to negligible particle–particle interaction, each particle can be considered to be suspended in an unbounded fluid. In their experiments with highly elastic fluids, Iso, Cohen & Koch (1996a) observed an isolated fibre to obtain a stable orientation close to the direction of the imposed flow. They found that fibres in a moderate particle concentration suspension also obtained a highly peaked orientation distribution near the flow direction for the same parameters. Therefore, the particle–particle interaction may sometimes be ignored even with higher particle concentrations, making the studies of an isolated particle/fibre freely rotating in a simple shear flow of a viscoelastic fluid invaluable.

Based on their qualitative nature, the Jeffery orbits, illustrated for particles with aspect ratio, $\kappa =10$ and 50 in figure 1, may be classified into log-rolling, wobbling, flipping and tumbling. Log-rolling occurs when the particle is aligned with the axis of rotation of the imposed flow or the vorticity axis. Here the particle rotates about its major axis. Except in the log-rolling motion, the particle's angular velocity changes throughout its orbit. When initially placed in the flow–gradient plane (FGP), the particle remains in the plane. It rotates about its minor axis while tumbling from one side of the flow axis to the other. In the tumbling and flipping orbits, a large aspect ratio particle or fibre spends only an $O(1/\kappa )$ proportion (non-dimensionalized with shear rate) of the Jefferey time-period of $2{\rm \pi} \kappa$ away from the flow–vorticity plane indicated by dashed red lines in the plots of figure 1. In the flipping orbits (which are three-dimensional extensions of the tumbling orbit when the particle is not confined to the FGP), the particle comes within $O(1/\kappa )$ of the flow direction. In these orbits, the particle's orientation rapidly flips from being aligned with the positive to the negative flow axis. During its rapid flipping phase, a particle in a flipping orbit spans a large portion of the orientation space in the gradient direction. In wobbling orbits, which are smaller in circumference than flipping orbits, the particle does not come very close to the flow direction. In these orbits, the particle gradually wobbles in its orbit around the vorticity axis.

Gauthier et al. (1971) conducted experiments with $\kappa =16.1$ nylon rods in a viscoelastic fluid made with 0.03 wt. % polyacrylamide solution in water. They conducted experiments at a shear rate of $0.53\ {\rm s}^{-1}$ and found that a particle starting close to the FGP spirals towards the vorticity axis as it is exposed to the Couette flow. Using a similar viscoelastic fluid and a polyethylene rod of $\kappa =9.0$ Bartram et al. (1975) also found a similar behaviour with shear rates up to $5\ {\rm s}^{-1}$. Upon further increase in shear rate, they found that a $\kappa =9.0$ rod released near the gradient direction initially moves within the FGP towards an orientation near the flow direction. From here, introducing a disturbance made the particle move monotonically along the flow–vorticity plane away from the flow direction. When the particle was sufficiently close to the vorticity direction, it started to spiral towards it. Bartram et al. (1975) observed similar behaviour with a $\kappa =5.6$ rod. However, unlike the $\kappa =9.0$ rod, no disturbance was required when the particle came close to the flow direction after being placed near the gradient direction. The time period of particle rotation about the vorticity axis, that is already very large at large $\kappa$ in Newtonian fluid ($2{\rm \pi} \kappa$), is further increased in experiments with viscoelastic fluids (Gauthier et al. 1971; Bartram et al. 1975). In the experiments where the orientation of the particle centreline was found to be spiralling towards the vorticity axis, complete alignment with the vorticity axis was not shown.

Iso et al. (1996b) observed rotations of different high $\kappa$ isolated fibres in a Boger fluid consisting of 1000 ppm polyisobutylene (PIB) in polybutene (PB) in a simple shear flow with different shear rates. The polymer relaxation time and concentration, $c$, defined as the polymer-to-solvent viscosity ratio, were 3$s$ and 0.39, respectively. With a $\kappa =19$ fibre, they found the particle to be orientated very close to the vorticity axis when $De=1.5$. Increasing $\kappa$ to 34.4, or $De$ to 3.0 or both, they found that, after initial spiralling away from the FGP, the particle obtained a steady orientation between $5^\circ$ and $60^\circ$ from the vorticity axis in the flow–vorticity plane. With $\kappa =34.4$ and $De=3.0$, they report two additional observations with no initial spiralling in contrast to other experiments at identical parameters. The authors attributed different initial orientations and fluid rheology due to slight changes in room temperature as the causes of the lack of initial spiralling.

Across their two studies Iso et al. (1996a,b) also conducted experiments in a viscoelastic liquid obtained by adding a certain amount of high molecular weight polymer polyacrylamide (PAA) to a Newtonian solvent. The shear rate in these experiments was $0.5\ {\rm s}^{-1}$, and the fluid was slightly shear-thinning. They observed various behaviours as the amount of PAA was increased from 100 ppm to 2000 ppm (although the exact value of $c$ is not available, it increases with PAA amount). For 100 ppm ($\kappa =14$) and 500 ppm ($\kappa =24$) solutions of PAA at a shear rate of $0.5\ {\rm s}^{-1}$ (Iso et al. 1996a,b), fibres either end up in a trajectory where they oscillate in a small periodic orbit close to the vorticity axis or obtain a stable orientation in the flow–vorticity plane at a particular angle from the vorticity axis similar to the Boger fluid (Boger 1977) experiments at a higher shear rate of $1.0\ {\rm s}^{-1}$ (Iso et al. 1996b). With 1000 ppm PAA, the $\kappa =24$ fibre obtains a stable orientation at the flow direction or $20^\circ$ from the flow direction in the flow–vorticity plane. Irrespective of the initial condition, a $\kappa =24$ fibre in 2000 ppm PAA solution stably aligns with the flow direction. Therefore, fibres become more flow aligned with increasing elasticity or polymer concentration.

The latest available experimental results measuring the effect of viscoelasticity on the rotation of an anisotropic particle in simple shear flow are by Gunes et al. (2008). They considered hematite spheroidal particles with a much smaller aspect ratio, $\kappa$, between 2 and 7.5, than the previous experimental studies. In a 20 % hydroxypropylcellulose solution in water, they found $\kappa =3.8$ particles to be oriented close to vorticity and flow directions at low and large shear rates or $De$, respectively. At moderate shear rates, particles exhibited a bimodal orientation distribution. In a 2 % poly-(ethyleneoxide), flow alignment was obtained at a lower shear rate for a higher $\kappa$ or $c$. Most of the fluids they used were shear-thinning. They reported one experiment of a non-shear-thinning Boger fluid (Boger 1977), consisting of a 0.1 % polyisobutylene in polybutadiene solution, with $\kappa =3.8$. Here the particles were close to the vorticity axis at all shear rates reported.

Due to the variety of non-Newtonian parameters needed to fully characterize a viscoelastic fluid and several physical phenomena that polymers may undergo simultaneously, it is difficult to quantitatively compare the previous experiments and identify the source of various behaviours of particle orientation. For example, temperature changes during an experiment may change the rheology of the fluid, subsequently changing the polymer's relaxation time, or multiple relaxation times may be needed to represent the fluid fully, or adding more polymers to a solution may not only increase the polymer concentration, but it may also change the relaxation time of polymers as they entangle with one another. Thus, numerical and theoretical modelling of the relevant system is an important tool in isolating and understanding different physical mechanisms that can complement or inspire future experiments.

Recently, d'Avino et al. (2014) reported numerical simulations of a $\kappa =4.0$ prolate spheroidal particle in a simple shear flow of a Giesekus fluid that models polymer melts (Bird, Armstrong & Hassager 1987), with large polymer concentration. They used $c=10$ and found spiralling towards the vorticity axis for $De\lesssim 1$ and alignment close to the flow direction for $De\gtrsim 3$. For intermediate $De$, depending on $De$, the particle obtained either one or two stable orientations between the flow and vorticity axis. They mentioned similar observations in unreported simulations with $\kappa =8$ spheroids with the transition from vorticity-aligned to flow-aligned particle orientation occurring at a smaller $De$. Therefore, the trend of the particle's final orientation with the shear rate or $De$ and $\kappa$ between the experiments of Gunes et al. (2008) and the numerical study of d'Avino et al. (2014), both conducted at small $\kappa$, is similar. However, it is unclear from the aforementioned numerical findings if the novel orientation dynamics are due to first or second normal stress difference, shear thinning or synergistic effects of various non-Newtonian behaviours exhibited by the Giesekus fluid. Also, at intermediate shear rates or $De$, while the orientation behaviour is bimodal, i.e. either along vorticity or flow directions, in the experiments of Gunes et al. (2008), the final orientations in numerical results of d'Avino et al. (2014) lie between the flow and vorticity axes in the flow–vorticity plane. The latter is instead similar to some of the experimental observations of Iso et al. (1996a) at larger $\kappa$.

A richer orientation behaviour is illustrated in the previous experiments of Gauthier et al. (1971), Bartram et al. (1975) and Iso et al. (1996a,b) at larger $\kappa$ as compared with the more recent studies of Gunes et al. (2008) and d'Avino et al. (2014) at smaller $\kappa$. Numerical studies with large $\kappa$ in a viscoelastic fluid are expensive due to the large velocity and polymer conformation gradients near the particle surface. Resolving these large gradients and accurately modelling the shape of a slender particle requires very fine spatial resolution and hence smaller time steps to ensure numerical stability. Therefore, theoretical studies at large $\kappa$ are invaluable, and Leal (1975), Harlen & Koch (1993) and Abtahi & Elfring (2019) are such pre-existing examples.

Using the slender body theory of Batchelor (1970), Leal (1975) predicts that a fibre rotating in a slow, simple shear flow of a second-order fluid, will spiral towards the vorticity axis due to the second normal stress difference, $\psi _2$, of the fluid. Here $\psi _2$ is usually smaller than the first normal stress difference, $\psi _1$, and it is zero for Boger fluids (Boger 1977). Hence, Leal's theory predicts that a slender particle rotating in a Boger fluid undergoing a simple shear flow with a small shear rate rotates in the same manner as in a Newtonian fluid. However, the low shear rate experiments of Iso et al. (1996b) with a Boger fluid show spiralling towards vorticity. For a large Deborah number, $De\gg 1$, small polymer concentration, $c\ll 1$ and also using the slender body theory of Batchelor (1970), Harlen & Koch (1993) predict the fibre in an Oldroyd-B fluid to spiral towards the vorticity axis, but, due to $\psi _1$ (an Oldroyd-B fluid has $\psi _2=0$). Shear-thinning, another property exhibited by polymeric fluids, does not play a role in either of these theories. Abtahi & Elfring (2019) conducted a theoretical study on an asymptotically weakly shear thinning liquid and found that a prolate spheroid rotates in a closed periodic orbit but with a longer time period compared with the Jeffery orbit in a Newtonian fluid. In other words, shear thinning makes a prolate spheroid's trajectory equivalent to that of a larger aspect ratio particle in a Newtonian fluid but does not qualitatively alter the topology of the trajectories.

Spiralling towards vorticity, as indicated by the two previous theories (Leal 1975; Harlen & Koch 1993) predicting a qualitative change in particle orientation trajectory, is only one of the many qualitative influences of viscoelasticity observed in the previous simple shear experiments of Gauthier et al. (1971), Bartram et al. (1975) and Iso et al. (1996a,b) at large particle aspect ratio, $\kappa$. In this paper, assuming a small polymer concentration, we theoretically revisit the effect of viscoelasticity on a slender fibre rotating in a simple shear flow to explain the richer qualitative behaviour of a large $\kappa$ particle's orientation in a polymeric fluid observed in the experiments. We use the Oldroyd-B model to isolate the effect of elasticity from shear thinning. Also, any predicted viscoelastic effects will exclude the impact of the second normal stress difference and re-examine the fibre rotation in Boger fluid undergoing simple shear flow.

In the absence of fluid inertia, fluid stress at any point in the viscoelastic fluid surrounding a suspended particle is decomposed into three components: (a) particle motion induced solvent stress, i.e. the stress around the particle rotating in a Newtonian fluid, (b) elastic stress or the polymer stress and (c) polymer-induced solvent stress, i.e. the stress created by the perturbations in fluid velocity and pressure due to the forcing by polymer stress. Therefore, the torque acting on the particle is the sum of particle motion-induced solvent torque (MIST), elastic torque and polymer-induced solvent torque (PIST). We consider a freely rotating or torque-free particle where the sum of the three torque components is zero.

The rest of the paper is organized as follows. In § 2 we discuss different torque generating mechanisms along with the mathematical formulation relevant to any freely rotating (torque-free) particle in a viscoelastic fluid. For any polymer concentration, $c$, using a generalized reciprocal theorem, we derive the formulation where PIST on any such particle can be expressed in terms of the polymer stress instead of the polymer-induced solvent stress. After § 2 we are concerned with viscoelastic fluid with a small polymer concentration, $c$. Therefore, in § 3 we briefly describe the regular perturbation expansion of the relevant flow variables in $c$ and the procedure to obtain the $O(c)$ correction to the rotation of a torque-free particle in a low $c$ viscoelastic fluid. The formulation in § 2 that expresses PIST in terms of the polymer stress allows us to circumvent the numerical discretization of the partial differential equations to calculate the $O(c)$ polymer-induced solvent stress. Therefore, the $O(c)$ PIST (and also the elastic torque) can be evaluated using the leading order or Newtonian velocity field around the particle. In § 4, we use the matched asymptotic/slender body theory (SBT) solution for the Newtonian flow field around a slender prolate spheroid to calculate the torques and obtain the $O(c)$ correction to the Jeffery rotation (Jeffery 1922) rates for large aspect ratio prolate spheroids due to viscoelasticity. In SBT, in the inner region close to the particle, owing to a large $\kappa$, the velocity field is obtained by assuming the flow to vary slowly along the length of the particle compared with its variation perpendicular to the particle surface. This solution is taken from Cox (1970). Farther from the particle surface, in the outer region, the particle centre line is replaced by a line of Stokeslets and doublets. The velocity disturbance created by these are taken from Batchelor (1970) and Cox (1971), respectively. In the SBT (Cox 1970, 1971; Batchelor 1970), the inner and the outer solution approach one another in the matching region, i.e. in the dual limit of the radial distance from the particle centreline written in inner and outer variables approaching infinity and zero, respectively. In § 5 we study the dynamical system defined by these equations for different $c$ and $De$ and illustrate various orientation dynamics of the spheroid predicted by this system. We also provide a qualitative comparison of the theoretical orientation trajectories with the previous experimental observations. Finally, we conclude our findings in § 6.

2. Mathematical formulation and different torque generating mechanisms in viscoelastic fluids

In the absence of inertia, equations governing the conservation of mass and momentum in the viscoelastic fluid surrounding a particle are

(2.1a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}=0, \end{equation}

where $\boldsymbol {u}$ and $\boldsymbol {\sigma }$ are the fluid velocity vector and stress tensor field. We consider a particle with its centre of mass at the origin of a simple shear flow such that it freely rotates with an angular velocity $\boldsymbol {\omega }_p$, but does not translate. Therefore, the no-slip boundary condition on the particle surface and the far-field (as $|\boldsymbol {r}|\rightarrow \infty$) imposed flow boundary condition are

(2.2)\begin{equation} {\boldsymbol{u}=\boldsymbol{\omega}_p\times\boldsymbol{r}, \quad\text{on particle surface}, \quad \text{and,}\quad \boldsymbol{u}=\boldsymbol{r}\boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{u})_\infty, \quad \text{as }|\boldsymbol{r}|\rightarrow\infty,} \end{equation}

where $(\boldsymbol {\nabla }\boldsymbol {u})_\infty$ is the velocity gradient tensor of the imposed flow. The equations are non-dimensionalized with the particle's maximum length and the inverse of the imposed shear rate as length and time scales. The fluid stress,

(2.3)\begin{equation} \boldsymbol{\sigma}=\boldsymbol{\tau}+\boldsymbol{\varPi}, \end{equation}

is the sum of the solvent stress, $\boldsymbol {\tau }=-p\boldsymbol {I}+\boldsymbol {\nabla } \boldsymbol {u}+(\boldsymbol {\nabla } \boldsymbol {u})^{\rm T}$ and the polymeric stress, $\boldsymbol {\varPi }$. The solvent viscosity is the viscosity scale in our non-dimensionalization. In the solvent stress, $p$ is the fluid pressure. We use the Oldroyd-B model (Bird et al. 1987) where the polymers are modelled as dumbbells consisting of two beads attached to a linearly elastic spring. The polymeric stress is

(2.4)\begin{equation} \boldsymbol{\varPi}=\frac{c}{De}(\boldsymbol{\varLambda}-\boldsymbol{I}),\end{equation}

where $c$ is the polymer concentration, $De$ is a non-dimensional product of the polymer relaxation time and imposed shear rate, $\boldsymbol {\varLambda }$ is the polymer conformation (defined as the average over possible polymer conformations of the outer product of the polymer end to end vector with itself) and $\boldsymbol {I}$ is the identity tensor. Here $\boldsymbol {\varLambda }$ is affected by convection due to the fluid velocity, stretching and rotation by the velocity gradients, and relaxation of the polymer to its equilibrium orientation, $\boldsymbol {\varLambda }=\boldsymbol {I}$. It is governed by

(2.5)\begin{equation} \frac{\partial \boldsymbol{\varLambda}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varLambda}=(\boldsymbol{\nabla} \boldsymbol{u})^{\rm T}\boldsymbol{\cdot}\boldsymbol{\varLambda}+ \boldsymbol{\varLambda}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}+\frac{1}{De}(\boldsymbol{I}-\boldsymbol{\varLambda}).\end{equation}

Due to the absence of any nonlinear term explicitly involving the velocity and pressure variables in the mass and momentum conservation equations, we can split (2.1a,b) and its boundary conditions into two parts. The first part is the same as the flow around a particle rotating with an angular velocity $\boldsymbol {\omega }_p$ in an imposed simple shear flow of a Newtonian fluid. It is governed by

(2.6a,b)\begin{equation} {\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{M}=0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau}^{M}=0, } \end{equation}

where the particle motion induced solvent stress, $\boldsymbol {\tau }^{M}=-p^{M}\boldsymbol {I}+\boldsymbol {\nabla }\boldsymbol {u}^{M}+(\boldsymbol {\nabla } \boldsymbol {u}^{M})^{\rm T}$, is the sum of the pressure and rate of strain tensor in the solvent due to the particle motion in an inertialess Newtonian fluid. These equations are subject to boundary conditions

(2.7)\begin{equation} \boldsymbol{u}^{M}=\boldsymbol{\omega}_p\times\boldsymbol{r}, \quad \text{on particle surface},\quad \text{and,}\quad \boldsymbol{u}^{M}=\boldsymbol{r}\boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{u})_\infty, \quad\text{as }|\boldsymbol{r}|\rightarrow\infty.\end{equation}

The second part is the balance of the divergence of the sum of the polymeric stress and $\boldsymbol {\tau }^{P}$,

(2.8a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{P}=0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau}^{P}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varPi}=0, \end{equation}

where the solvent stress generated due to the forcing by the polymeric stress is $\boldsymbol {\tau }^{P}=-p^{P}\boldsymbol {I}+\boldsymbol {\nabla } \boldsymbol {u}^{P}+(\boldsymbol {\nabla } \boldsymbol {u}^{P})^{\rm T}$. Here $p^{P}$ and $(\boldsymbol {\nabla } \boldsymbol {u}^{P}+(\boldsymbol {\nabla } \boldsymbol {u}^{P})^{\rm T})/2$ are the modification of pressure and rate of strain tensor by the polymers. Boundary conditions for (2.8a,b) are

(2.9)\begin{equation} \boldsymbol{u}^{P}=0, \quad\text{on particle surface},\quad \text{and,}\quad \boldsymbol{u}^{P}=0, \quad \text{as }|\boldsymbol{r}|\rightarrow\infty.\end{equation}

The two parts are coupled via (2.5), i.e. the polymer constitutive equation, where

(2.10)\begin{equation} \boldsymbol{u}=\boldsymbol{u}^{M}+\boldsymbol{u}^{P},\end{equation}

and the torque-free condition (2.17). In this framework, the total fluid stress,

(2.11)\begin{equation} \boldsymbol{\sigma}=\boldsymbol{\tau}^{M}+\boldsymbol{\tau}^{P}+\boldsymbol{\varPi}, \end{equation}

and the total torque acting on the particle surface, $\boldsymbol {\sigma }$,

(2.12)$$\begin{gather} \boldsymbol{G}(\boldsymbol{\sigma})=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times( \boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n}), \end{gather}$$

(2.13)$$\begin{gather}\boldsymbol{G}(\boldsymbol{\sigma})=\boldsymbol{G}(\boldsymbol{\tau}^{M})+ \boldsymbol{G}(\boldsymbol{\tau}^{P})+\boldsymbol{G}(\boldsymbol{\varPi})=\boldsymbol{G}^{MIST}+ \boldsymbol{G}^{PIST}+\boldsymbol{G}^{Elastic} \end{gather}$$

are decomposed into three underlying mechanisms, where

(2.14)$$\begin{gather} \boldsymbol{G}^{MIST}=\boldsymbol{G}({\boldsymbol{\tau}^{M}})=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times \boldsymbol{\tau}^{M}\boldsymbol{\cdot} \boldsymbol{n}, \end{gather}$$

(2.15)$$\begin{gather}\boldsymbol{G}^{PIST}=\boldsymbol{G}({\boldsymbol{\tau}^{P}})=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times \boldsymbol{\tau}^{P}\boldsymbol{\cdot} \boldsymbol{n} \end{gather}$$

and

(2.16)\begin{equation} \boldsymbol{G}^{Elastic}=\boldsymbol{G}(\boldsymbol{\varPi})=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times\boldsymbol{\varPi}\boldsymbol{\cdot} \boldsymbol{n} \end{equation}

are the particle motion induced solvent, the polymer-induced solvent, and the elastic or polymeric torques, respectively. The angular velocity, $\boldsymbol {\omega }_p$, of a freely rotating particle, is determined by the torque-free condition

(2.17)\begin{equation} \boldsymbol{\omega}_p=\{\boldsymbol{\omega}_p: \boldsymbol{G}=\boldsymbol{G}^{MIST}+\boldsymbol{G}^{PIST}+\boldsymbol{G}^{Elastic}=0\}.\end{equation}

We consider the motion of a freely rotating particle in a viscoelastic fluid. Our main motivation is to study the behaviour of a prolate spheroid in simple shear flow. Due to symmetry, this particle has zero net hydrodynamic force if it moves with the local velocity of the imposed flow. Its physically motivated decomposed components (particle motion induced solvent, polymer-induced solvent and elastic forces) are individually zero by symmetry. However, the force balance can be considered similar to the torque balance discussed above for determining the motion of a freely translating particle (of any shape) in a general linear flow or a particle sedimenting under gravity (where the net hydrodynamic force must balance the force due to gravity).

2.1. Using a generalized reciprocal theorem to obtain ${\boldsymbol {G}^{PIST}}$ in terms of $\boldsymbol {\varPi }$

In its original form, $\boldsymbol {G}^{PIST}=\int _{\boldsymbol {r}_{p}} {\rm d}S\,\boldsymbol {r}\times \boldsymbol {\tau }^{P}\boldsymbol{\cdot} \boldsymbol {n}$ is a function of $\boldsymbol {\tau }^{P}$ which in turn is driven by $\boldsymbol {\varPi }$ through (2.8a,b). In this section, using the balance of the divergence of the polymeric stress ($\boldsymbol {\varPi }$) and $\boldsymbol {\tau }^{P}$ from (2.8a,b) and using a generalized reciprocal theorem we derive an expression to obtain $\boldsymbol {G}^{PIST}$ directly from $\boldsymbol {\varPi }$ without the need to compute $\boldsymbol {\tau }^{P}$.

The auxiliary or comparison problem in a generalized reciprocal theorem must be chosen based on the surface integral one is interested in evaluating. The effect of the torque is to rotate the particle. Hence, we consider the Stokes flow around a particle rotating in a quiescent Newtonian fluid. We define the following auxiliary Stokes problem for a ‘velocity’ field consisting of a rank-two pseudotensor $\boldsymbol {b}$, a ‘pressure’ field that is a pseudovector, $\boldsymbol {q}$, and a ‘fluid stress’ field that is a rank-three pseudotensor $\boldsymbol {B}$,

(2.18a–c)\begin{equation} \frac{\partial B_{ijk}}{\partial r_i}=0,\quad \frac{\partial b_{ij}}{\partial r_i}=0,\quad B_{ijk}=-\delta_{ij}q_k+\frac{\partial b_{jk}}{\partial r_i}+\frac{\partial b_{ik}}{\partial r_j},\end{equation}

with boundary condition

(2.19)\begin{equation} b_{ij}=\epsilon_{ijk}r_k, \quad \text{on particle surface}, \quad \text{and}, \quad b_{ij}=0, \quad \text{as } |\boldsymbol{r}|\rightarrow\infty, \end{equation}

where $\epsilon _{ijk}$ is the permutation tensor. Here $\boldsymbol {b}\boldsymbol{\cdot}\boldsymbol {\omega }_{auxiliary}$ is the velocity field around a particle rotating with an angular velocity $\boldsymbol {\omega }_{auxiliary}$ in a quiescent Newtonian fluid. Using the definitions of $\boldsymbol {\tau }^{P}$ and $\boldsymbol {B}$ in terms of the respective velocities, $\boldsymbol {u}^{P}$ and $\boldsymbol {b}$, incompressibility of the velocities in (2.8a,b) and (2.18a–c), $\boldsymbol {\nabla }\boldsymbol{\cdot}\boldsymbol {B}=0$ from (2.18a–c) and the symmetry of $B_{lki}$ and ${\tau }^{P}_{lk}$ about the $l$ and $k$ indices we obtain

(2.20)\begin{equation} {{\tau}^{P}_{lk}\frac{\partial b_{ki}}{\partial r_l}=\frac{\partial B_{lki} {u}^{P}_k}{\partial r_l}.}\end{equation}

Using the balance of the divergence of $\boldsymbol {\varPi }$ and $\boldsymbol {\tau }^{P}$ from (2.8a,b), the volume integral of (2.20) in a region bounded by the particle surface, $\boldsymbol {r}_{p}$, and a far-field spherical surface at $|\boldsymbol {r}|\rightarrow \infty$ is

(2.21)\begin{equation} {\int_{Fluid} {\rm d}V\,\frac{\partial}{\partial r_l}[{\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k]=- \int_{Fluid} {\rm d}V\,b_{ki}\frac{\partial{\varPi}_{lk}}{\partial r_l}.} \end{equation}

Using the divergence theorem, the left-hand side of the above equation can be written as

(2.22)\begin{align} \int_{Fluid} {\rm d}V\,\frac{\partial}{\partial r_l}[{\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k]& = \int_{|\boldsymbol{r}|\rightarrow\infty} {\rm d}S\,n_l[{\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k]\nonumber\\ &\quad -\int_{\boldsymbol{r}_{p}}{\rm d}S\,n_l[{\tau}^{P}_{lk}b_{ki}- B_{lki}{u}^{P}_k], \end{align}

where the surface normal $n_l$ points into the fluid (away from particle region) on the particle surface. The fluid velocity due to a particle that exerts a force (force dipole) on the fluid decays as $1/r$ ($1/{r}^2$) in the far field. Hence, the velocities $\boldsymbol {u}^{P}$ and $\boldsymbol {b}$ scale as $1/r$ and $1/{r}^2$, respectively, and the stresses $\boldsymbol {\tau }^{P}$ and $\boldsymbol {B}$ scale as $1/{r}^2$ and $1/{r}^3$. Therefore, the first surface integral in (2.22) vanishes. Using the boundary conditions for $\boldsymbol {u}^{P}$ and $\boldsymbol {b}$ from (2.9) and (2.19) in the second surface integral in (2.22) we obtain

(2.23)\begin{equation} \int_{Fluid} {\rm d}V\,\frac{\partial}{\partial r_l}({\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k)= \int_{\boldsymbol{r}_{p}}{\rm d}S\,n_l{\tau}^{P}_{lk}\epsilon_{kim}r_m=-{G}^{PIST}_i.\end{equation}

Therefore, from (2.21) and (2.23), for a particle of any shape,

(2.24)\begin{equation} \boldsymbol{G}^{PIST}=\int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\varPi}}) \boldsymbol{\cdot}\boldsymbol{b}.\end{equation}

This completes the derivation expressing $\boldsymbol {G}^{PIST}$ in terms of the polymer stress $\boldsymbol {\varPi }$ and a quasisteady two-tensor field $\boldsymbol {b}$ that is dependent on the particle shape and is a solution of the auxiliary Stokes problem defined by (2.18a–c) and (2.19).

In a linear imposed flow such as a simple shear, the undisturbed polymer stress, $\boldsymbol {\varPi }^U$, i.e. the polymer stress without the particle, is spatially constant. Hence,

(2.25)\begin{equation} \boldsymbol{G}^{PIST}=\int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}( \boldsymbol{\varPi}-\boldsymbol{\varPi}^U))\boldsymbol{\cdot}\boldsymbol{b}.\end{equation}

Using the chain rule and divergence theorem,

(2.26)\begin{align} {G}^{PIST}_i&=\int_{|\boldsymbol{r}|\rightarrow\infty} {\rm d}S({\varPi}_{lk}-{\varPi}^{U}_{lk}) b_{ki} n_l-\int_{\boldsymbol{r}_{p}} {\rm d}S({\varPi}_{lk}-{\varPi}^{U}_{lk}) b_{ki}n_l\nonumber\\ &\quad -\int_{Fluid} {\rm d}V({\varPi}_{lk}-{\varPi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_l}. \end{align}

The disturbance of the polymer stress created by the particle, ${\varPi }_{lk}-{\varPi }^{U}_{lk}$, also decays as $1/{r}^2$ in the far-field since it is forced by the disturbance to the far-field velocity gradients. In the case of Oldroyd-B fluids or other dumbbell models such as FENE-P and Giesekus (Bird et al. 1987), this far-field scaling of ${\varPi }_{lk}-{\varPi }^{U}_{lk}$ is ascertained by linearizing the polymer constitutive equation (for example (2.4) and (2.5) for the Oldroyd-B model) about the far-field velocity and polymer conformation to obtain a governing equation for ${\varPi }_{lk}-{\varPi }^{U}_{lk}$. Therefore, the first surface integral in the above equation vanishes, and using the surface boundary condition of (2.19) leads to

(2.27)\begin{equation} {G}^{PIST}_i=-\int_{\boldsymbol{r}_{p}} {\rm d}S({\varPi}_{lk}-{\varPi}^{U}_{lk}) \epsilon_{kif}r_fn_l-\int_{Fluid}{\rm d}V({\varPi}_{lk}-{\varPi}^{U}_{lk}) \frac{\partial b_{ki}}{\partial r_l}, \end{equation}

and using ${G}^{Elastic}_i=\int _{\boldsymbol {r}_{p}} {\rm d}S\,{\varPi }_{lk} \epsilon _{kif}r_fn_l$,

(2.28)\begin{equation} {G}^{PIST}_i+{G}^{Elastic}_i=\int_{\boldsymbol{r}_{p}} {\rm d}S\,{\varPi}^{U}_{lk} \epsilon_{kif}r_fn_l-\int_{Fluid} {\rm d}V({\varPi}_{lk}-{\varPi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_l}.\end{equation}

The first term on the right-hand side is the torque on a particle about its centre of mass due to (constant) undisturbed polymer stress acting on its surface. It is zero by symmetry for a constant density fore-and-aft and axisymmetric particle such as a slender prolate spheroid, and hence for such a particle,

(2.29)\begin{equation} {\boldsymbol{G}^{PIST}}+{\boldsymbol{G}^{Elastic}}=-\int_{Fluid} {\rm d}V (\boldsymbol{\varPi}-\boldsymbol{\varPi}^{U}):\boldsymbol{\nabla}\boldsymbol{b}\end{equation}

and

(2.30)\begin{equation} {\boldsymbol{G}^{Elastic}}=-\int_{Fluid}{\rm d}V\,\boldsymbol{\nabla}\boldsymbol{\cdot} ((\boldsymbol{\varPi}-\boldsymbol{\varPi}^{U})\boldsymbol{\cdot}\boldsymbol{b}). \end{equation}

Equation (2.28) (or (2.29) for a fore-and-aft and axisymmetric particle) represents the total torque acting on the particle due to the presence of the polymers (including effects from both the polymeric stress and the polymer-induced solvent stress), as a function of polymer stress, $\boldsymbol {\varPi }$, and its undisturbed value, $\boldsymbol {\varPi }^U$. For the Oldroyd-B fluid in a simple shear flow with 1, 2 and 3 as flow, gradient and vorticity directions, respectively, in Cartesian coordinates, ${\varPi }^{U}_{ij}=c(2De\delta _{i1}\delta _{j1}+(\delta _{i2}\delta _{j1}+\delta _{i1}\delta _{j2}))$.

3. Regular perturbation expansion for small polymer concentration

The leading-order solution in a regular perturbation expansion for $c\ll 1$ corresponds to a freely rotating particle in simple shear flow of a Newtonian fluid. At this order the stresses and hence the respective torques arising due to the polymers, i.e. ${\boldsymbol {G}^{PIST}}$ and ${\boldsymbol {G}^{Elastic}}$ are zero and the particle rotates with an angular velocity that satisfies ${\boldsymbol {G}^{MIST}}=0$ (2.17). This is simply the Jeffery rotation. At the leading order the polymer configuration is driven by the leading-order velocity field (2.5) and this leads to an $O(c)$ polymer stress (2.4). Therefore, the torques ${\boldsymbol {G}^{PIST}}$ and ${\boldsymbol {G}^{Elastic}}$ are $O(c)$ ((2.25), (2.14), (2.15), (2.16)). Hence, the particle rotation must be modified at $O(c)$ such that the sum of all three torques ${\boldsymbol {G}^{MIST}}$, ${\boldsymbol {G}^{PIST}}$ and ${\boldsymbol {G}^{Elastic}}$ is zero at $O(c)$. The regular perturbation expansion of the relevant flow variables is

(3.1)$$\begin{gather} \boldsymbol{u}^{M}={\boldsymbol{u}^{M}}^{(0)}+c{\boldsymbol{u}^{M}}^{(1)}+O(c^2), \end{gather}$$

(3.2)$$\begin{gather}p^{M}={p^{M}}^{(0)}+c{p^{M}}^{(1)}+O(c^2), \end{gather}$$

(3.3)$$\begin{gather}\boldsymbol{\tau}^{M}={\boldsymbol{\tau}^{M}}^{(0)}+c{\boldsymbol{\tau}^{M}}^{(1)}+O(c^2), \end{gather}$$

(3.4)$$\begin{gather}\boldsymbol{u}^{P}=c{\boldsymbol{u}^{P}}^{(1)}+O(c^2), \end{gather}$$

(3.5)$$\begin{gather}p^{P}=c{p^{P}}^{(1)}+O(c^2), \end{gather}$$

(3.6)$$\begin{gather}\boldsymbol{\tau}^{P}=c{\boldsymbol{\tau}^{P}}^{(1)}+O(c^2), \end{gather}$$

(3.7)$$\begin{gather}\boldsymbol{\varLambda}=\boldsymbol{\varLambda}^{(0)}+c\boldsymbol{\varLambda}^{(1)}+O(c^2), \end{gather}$$

(3.8)$$\begin{gather}\boldsymbol{\varPi}=c\boldsymbol{\varPi}^{(1)}+O(c^2), \end{gather}$$

(3.9)$$\begin{gather}\boldsymbol{\omega}_p={\boldsymbol{\omega}^{(0)}_p}+c{\boldsymbol{\omega}^{(1)}_p}+O(c^2), \end{gather}$$

(3.10)$$\begin{gather}\boldsymbol{G}^{MIST}=c{\boldsymbol{G}^{MIST}}^{(1)}+O(c^2), \end{gather}$$

(3.11)$$\begin{gather}\boldsymbol{G}^{PIST}=c{\boldsymbol{G}^{PIST}}^{(1)}+O(c^2), \end{gather}$$

(3.12)$$\begin{gather}\boldsymbol{G}^{Elastic}=c{\boldsymbol{G}^{Elastic}}^{(1)}+O(c^2). \end{gather}$$

In an inertialess Newtonian fluid undergoing simple shear, a particle rotating at an angular velocity $\boldsymbol {\omega }^{(0)}_p$ generates the velocity field ${\boldsymbol {u}^{M}}^{(0)}$. As mentioned earlier, the leading-order angular velocity, ${\boldsymbol {\omega }^{(0)}_p}$, is the Jeffery rotation. It allows the leading-order torque-free condition, $\int _{\boldsymbol {r}_{p}} {\rm d}S\,\boldsymbol {r}\times {\boldsymbol {\tau }^{M}}^{(0)}\boldsymbol{\cdot} \boldsymbol {n}=0$ to be satisfied. The leading-order polymer constitutive equation is

(3.13)\begin{equation} \frac{\partial \boldsymbol{\varLambda}^{(0)}}{\partial t}+{\boldsymbol{u}^{M}}^{(0)}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varLambda}^{(0)}=(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)})^{\rm T}\boldsymbol{\cdot}\boldsymbol{\varLambda}^{(0)}+ \boldsymbol{\varLambda}^{(0)}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}^{M}}^{(0)}+ \frac{1}{De}(\boldsymbol{I}-\boldsymbol{\varLambda}^{(0)}).\end{equation}

Solving this equation, one obtains the $O(c)$ polymer stress,

(3.14)\begin{equation} \boldsymbol{\varPi}^{(1)}=\frac{1}{De}(\boldsymbol{\varLambda}^{(0)}-\boldsymbol{I}),\end{equation}

and hence the $O(c)$ polymer-induced solvent and elastic torques,

(3.15)$$\begin{gather} {\boldsymbol{G}^{PIST}}^{(1)}=\int_{Fluid} {\rm d}V\,\boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{\varPi}^{(1)}-\boldsymbol{\varPi}^{U}/c) \boldsymbol{\cdot} \boldsymbol{b}, \end{gather}$$

(3.16)$$\begin{gather}{\boldsymbol{G}^{Elastic}}^{(1)}=\int_{\boldsymbol{r}_{p}} {\rm d}S\, \boldsymbol{r}\times\boldsymbol{\varPi}^{(1)}\boldsymbol{\cdot} \boldsymbol{n}. \end{gather}$$

The $O(c)$ angular velocity, ${\boldsymbol {\omega }^{(1)}_p}$ is the one that allows the $O(c)$ torque-free condition to be satisfied,

(3.17)\begin{equation} {\boldsymbol{G}^{MIST}}^{(1)}=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times {\boldsymbol{\tau}^{M}}^{(1)}\boldsymbol{\cdot} \boldsymbol{n}=-{\boldsymbol{G}^{PIST}}^{(1)}-{\boldsymbol{G}^{Elastic}}^{(1)},\end{equation}

where ${\boldsymbol {\tau }^{M}}^{(1)}=-{p^{M}}^{(1)}\boldsymbol {I}+\boldsymbol {\nabla } {\boldsymbol {u}^{M}}^{(1)}+(\boldsymbol {\nabla }{\boldsymbol {u}^{M}}^{(1)})^{\rm T}$ is obtained from the $O(c)$ equations,

(3.18a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}^{M}}^{(1)}=0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\tau}^{M}}^{(1)}=0,\end{equation}

subject to the velocity boundary conditions

(3.19)\begin{equation} {\boldsymbol{u}^{M}}^{(1)}={\boldsymbol{\omega}^{(1)}_p}\times\boldsymbol{r}, \quad\text{on particle surface},\quad\text{and,}\quad {\boldsymbol{u}^{M}}^{(1)}=0, \quad\text{as }|\boldsymbol{r}|\rightarrow\infty.\end{equation}

Equations (3.18a,b) and (3.19) represent the rotation of the particle in a quiescent Newtonian fluid, i.e. a Stokes flow. Due to the linearity of the Stokes flow, ${\boldsymbol {G}^{MIST}}^{(1)}$ is a linear function of ${\boldsymbol {\omega }^{(1)}_p}$. The latter can be viewed as the additional angular velocity of the particle that generates a large enough viscous torque, ${\boldsymbol {G}^{MIST}}^{(1)}$, to balance the sum of polymer-induced solvent and elastic torque.

The formulation in this section is valid for any polymer constitutive model and particle shape. By using the formulation of $\boldsymbol {G}^{PIST}$ in (2.25) to express its $O(c)$ value in (3.16), we have avoided dealing with the $O(c)$ equation for the balance of the divergence of the polymer stress and the polymer induced solvent stress (obtained by regularly expanding (2.8a,b) in $c$). Otherwise, obtaining $\boldsymbol {G}^{PIST}$ would have required a numerical solution via discretization of the governing partial differential equations.

4. Rotation of a fibre due to simple shear flow in viscoelastic fluid with small polymer concentration, $c$

In this section, we calculate torques from the various physical mechanisms discussed in § 2 on a large aspect ratio prolate spheroid (considered a slender fibre) freely rotating in a simple shear flow of a viscoelastic fluid with a small polymer concentration, $c$, using the procedure indicated in § 3. These torques are then balanced to obtain the $O(c)$ correction to the particle's rotation rate due to the presence of the polymers. The Jeffery orbit period of a slender fibre with aspect ratio $\kappa$ is $2{\rm \pi} \kappa$ and the proportion of time spent outside $|p_2|>O(1/\kappa )$ is only $O(1/\kappa )$, where $p_2=0$ defines the flow–vorticity plane. Therefore, a slender fibre suspended in a Newtonian fluid spends most of its Jeffery orbit close to the flow–vorticity plane, where the particle rotation rate is very small. Hence, most of the elasticity influence arises when the particle is in this orientation, and polymer conformation when the fibre is close to the flow–vorticity plane is quasisteady. Thus, the polymer constitutive equation (3.13) is simplified to

(4.1)\begin{align} \frac{\partial \boldsymbol{\varLambda}^{(0)}}{\partial t}\approx 0\rightarrow {\boldsymbol{u}^{M}}^{(0)}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varLambda}^{(0)}\approx(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)})^{\rm T}\boldsymbol{\cdot}\boldsymbol{\varLambda}^{(0)}+ \boldsymbol{\varLambda}^{(0)}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}^{M}}^{(0)}+ \frac{1}{De}(\boldsymbol{I}-\boldsymbol{\varLambda}^{(0)}).\end{align}

The Stokes flow solution of a Newtonian fluid around a large aspect ratio or a slender prolate spheroid is analytically solvable via a matched asymptotic expansion in $1/\kappa$ called SBT. In SBT, the fluid velocity close to the particle or in the inner region is considered quasi-two-dimensional. It is solved by ignoring the end effects and treating the slender particle as an infinite cylinder. In terms of the radial distance from the centreline (non-dimensionalized with the major radius of the particle) $\rho$, the inner region is defined by $\rho \ll 1$. Away from the particle in the outer region, defined by $\rho \gg 1/\kappa$ ($1/\kappa$ is the minor radius), the particle is assumed to be a line of point forces (Batchelor 1970; Cox 1970) and force doublets (Cox 1971). These singularity solutions are used to represent the velocity and pressure fields. For a large enough $\kappa$ (which is required for SBT to be valid), an intermediate or matching region exists which is defined by $1/\kappa \ll \rho \ll 1$. In this region, the $\rho \rightarrow \infty$ asymptote of the flow in the inner region and $\rho \rightarrow 0$ asymptote of the flow in the outer region are matched. See previous SBT calculations in Cox (1970, 1971) and Batchelor (1970) for more details. In the rest of this section, we will use these SBT results to obtain the effect of viscoelasticity on the rotation of a slender prolate spheroid. In contrast to a fibre with blunt ends, such as a cylinder, a slender prolate spheroid is more convenient for analysis due to the absence of localized forces at the ends.

4.1. Flow of Newtonian fluid around a slender fibre and particle rotation rates

Due to the relative velocity between the particle's centreline and the imposed flow, the flow disturbance at $O(1/\log (\kappa ))$ and higher orders in $1/\log (\kappa )$ can be considered to be generated by a line of point forces or force Stokeslets located at the particle centreline. This flow is considered by the general SBT of Cox (1970) up to $O(1/\log (\kappa )^2)$. An equivalent theory by Batchelor (1970) is valid at all orders in $1/\log (\kappa )$ for a prolate spheroid. For quantitative accuracy, it is advantageous to use the Stokeslet distribution, defined as $\boldsymbol {h}^{(0)}$ below, from the SBT of Batchelor (1970) instead of Cox (1970). For a prolate spheroidal particle fixed in the flow–vorticity plane of a simple shear flow or rotating about its centreline in a quiescent fluid, the force Stokeslet from these two theories is zero. The flow driven by the local velocity gradients, which acts at $O(1/\kappa ^2)$, is dominant in these cases. Considering the velocity gradients in the far-field relative to the particle, Cox (1971) provides the solution for this flow. It is a combination of flows driven by $O(1/\kappa ^2)$ force Stokeslets ($\boldsymbol {f}^{(0)}$) and doublets ($\boldsymbol {\mathcal {G}}^{(0)}$).

The Newtonian or the leading-order (in $c$) outer flow generated by a fibre freely rotating with an angular velocity $\boldsymbol {\omega }^{(0)}_p$ in an imposed shear flow of Newtonian fluid with 1, 2 and 3 as the flow, gradient and vorticity directions, when in orientation $\boldsymbol {p}$,

(4.2)\begin{equation} \boldsymbol{p}=\left[\,p_1\quad p_2\quad p_3\right]^{\rm T}, \end{equation}

close to the flow–vorticity plane is obtained by superposition of flows from Batchelor (1970), and Cox (1971) with the assumption $p_2\ll 1$ and $\kappa \gg 1$. The outer flow in the presence of the fibre is

(4.3)\begin{equation} \boldsymbol{u}_{out}^{{M}{(0)}}(\boldsymbol{r};\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})=\boldsymbol{r}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{u})_\infty+{\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}} (\boldsymbol{r};\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p}).\end{equation}

The disturbance velocity field,

(4.4)\begin{align} {\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}}(\boldsymbol{r};\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})&= \frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\left\{\frac{\boldsymbol{h}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_p, \boldsymbol{p})}{\log(2\kappa)-1.5}+\frac{\boldsymbol{f}^{(0)}(\lambda;\boldsymbol{p})}{\kappa^2}\right\}\nonumber\\ &\quad \boldsymbol{\cdot}\left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}\nonumber\\ &\quad +\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\,\frac{1}{\kappa^2}\boldsymbol{\mathcal{G}}^{(0)}( \lambda;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p}):\boldsymbol{\nabla} \left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}, \end{align}

is the sum of flows due to the point forces distributions $\boldsymbol {h}^{(0)}$ and $\boldsymbol {f}^{(0)}$, and the force doublet distribution $\boldsymbol {\mathcal {G}}^{(0)}$ discussed above. The different source distributions are

(4.5)\begin{align} \left.\begin{gathered} \boldsymbol{h^{(0)}}(\lambda;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})=-4 {\rm \pi}\lambda [\boldsymbol{p}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})_\infty-\boldsymbol{\omega}^{(0)}_p\times\boldsymbol{p}]\boldsymbol{\cdot} \left(\frac{1}{2}\boldsymbol{pp}+(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p})\frac{\log(2\kappa)-1.5}{ \log(2\kappa)-0.5}\right),\\ \boldsymbol{f}^{(0)}(\lambda;\boldsymbol{p})=-4{\rm \pi}\lambda \left[0\quad p_1\quad 0\right]^{\rm T} \left(1-\frac{1}{\log(2\kappa)-0.5}\right)+O(\,p_2),\\ \boldsymbol{\mathcal{G}}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})=2{\rm \pi} (1-\lambda^2)\left(\begin{bmatrix} 0 & 1+\dfrac{p_3^2}{2} & 0\\ \dfrac{p_3^2}{2} & 0 & -\dfrac{1}{2}p_1p_3\\ 0 & -\dfrac{1}{2}p_1p_3 & 0 \end{bmatrix}+(\boldsymbol{\omega}^{(0)}_p\boldsymbol{\cdot} \boldsymbol{p})\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{p}\right)+O(\,p_2). \end{gathered}\right\} \end{align}

A torque-free spheroid with aspect ratio, $\kappa$ rotating in a simple shear flow has the exact result for the temporal evolution of the orientation vector (Jeffery 1922; Kim & Karrila 2013),

(4.6)\begin{equation} \dot{\boldsymbol{p}}^{(0)}=\boldsymbol{\omega}^{(0)}_p\times\boldsymbol{p}=\boldsymbol{\omega}_\infty \times\boldsymbol{p}+\frac{\kappa^2+1}{\kappa^2-1}(\boldsymbol{E}_\infty\boldsymbol{\cdot} \boldsymbol{p})\boldsymbol{\cdot}(\boldsymbol{I}-\boldsymbol{pp}), \end{equation}

where $\boldsymbol {\omega }_\infty$ and $\boldsymbol {E}_\infty$ are the vorticity vector and strain rate tensor of the imposed simple shear,

(4.7a,b)\begin{equation} \boldsymbol{\omega}_\infty=-\tfrac{1}{2}\left[ 0\quad 0\quad 1 \right]^{\rm T},\quad \boldsymbol{E}_\infty=\tfrac{1}{2}\begin{bmatrix}0 & 1 & 0\\1 & 0 & 0\\0 & 0 & 0\end{bmatrix}. \end{equation}

For a slender spheroid close to the flow–vorticity plane,

(4.8)\begin{equation} \boldsymbol{\omega}^{(0)}_p\boldsymbol{\cdot}\boldsymbol{p}=-\frac{p_3}{2}+O(\,p_2).\end{equation}

Using $p_3=1$ in the above equation, we obtain the log-rolling velocity (when the particle is oriented with the vorticity axis) of the particle in a Newtonian fluid and find that in this orientation, a slender particle rotates at the angular velocity of the fluid, i.e. half of the shear rate. $\boldsymbol {\omega }^{(0)}_p\boldsymbol{\cdot}\boldsymbol {p}$ from (4.8) and $\boldsymbol {\omega }^{(0)}_p\times \boldsymbol {p}$ from (4.6) are used to obtain $\boldsymbol {\mathcal {G}}^{(0)}$ and $\boldsymbol {h^{(0)}}$, respectively, in (4.5).

The flow due to $\boldsymbol {h}^{(0)}$ is taken from Batchelor (1970). It is obtained by an expansion in $1/\log (\kappa )$ for a slender particle. But for a slender prolate spheroid, this expansion terminates such that the flow due to $\boldsymbol {h}^{(0)}$ captures the disturbance created by the prolate spheroid at all orders in $1/\log (\kappa )$ when expressed as in (4.4) and (4.6). The flow generated at the next order in $\kappa$ arises at $O(\kappa ^{-2})$ and is generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$. It is taken from Cox (1971) where only the $O(\kappa ^{-2})$ flow is available. Thus, the disturbance velocity in the outer region given by (4.4) has an overall error of $O(\kappa ^{-3})$. If only $\boldsymbol {h}^{(0)}$ is considered, the error is of $O(\kappa ^{-2})$. The flow generated by $\boldsymbol {h}^{(0)}$ is proportional to $p_2$ and therefore, no flow is produced by $\boldsymbol {h}^{(0)}$ when the particle is aligned in the flow–vorticity plane. In this plane $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ capture the (highest) $O(\kappa ^{-2})$ disturbance created by the particle. Accounting for the flow generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ allows us to consider the influence of elasticity within the flow–vorticity plane.

As mentioned just before (4.1) we expect most of the changes in the particle's rotation rate due to elasticity to arise when the particle is near the flow–vorticity plane, i.e. when $|p_2|\le O(1/\kappa )$. Therefore, in $\boldsymbol {h}^{(0)}$ we only consider the flow at $O(\,p_2)$, i.e. a flow of $O(\,p_2/\log (\kappa ))$ with an error of $O(\,p_2^2/\log (\kappa ))$. Since the primary purpose of using the $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ flow is to capture the finite effect of elasticity when the particle is in the flow–vorticity plane, we use $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ with $p_2=0$ (as expressed in (4.5)) which leads to a flow of $O(1/\kappa ^2)$ with an error of $O(\,p_2/\kappa ^2)$. When $|p_2|/\log (\kappa )$ is more than $1/\kappa ^2$, $\boldsymbol {h}^{(0)}$ driven flow dominates. In the $|p_2|\le O(1/\kappa )$ regime considered, the errors in $\boldsymbol {h}^{(0)}$ driven flow, i.e. $O(\,p_2^2/\log (\kappa ))$, are always smaller than the flow generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$, i.e. $O(1/\kappa ^2)$. As $p_2$ approaches zero, the $\boldsymbol {h}^{(0)}$ driven flow and the associated errors fall rapidly to zero making $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ driven flow at $O(1/\kappa ^2)$ the dominating one. Hence, the error terms arising from either of the $\boldsymbol {h}^{(0)}$ or $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ driven flow are always lower than the actual flow when $|p_2|\le O(1/\kappa )$. For $|p_2|>O(1/\kappa )$, the Newtonian rotation rate of the fibre dominates over the changes due to elasticity, and we consider the exact Jeffery rotation to account for it.

We have considered spheroidal particles in our theory. However, in experiments with slender particles (Gauthier et al. 1971; Bartram et al. 1975; Iso et al. 1996a,b) it is convenient to fabricate cylindrical particles. The forces generated by the blunt ends of a slender cylinder lead to an additional torque of $O(1/\kappa ^2)$ (Cox 1971) rendering the $O(1/\kappa ^2)$ flow generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ inaccurate. Instead of being valid at all orders in $1/\log (\kappa )$, the $\boldsymbol {h}^{(0)}$ generated flow may be considered up to order $1/\log (\kappa )^3$ (Batchelor 1970). Taking the $O(\,p_2)$ terms with this velocity field can still allow us to consider the flow accurately up to $O(\,p_2/\log (\kappa ))$. Thus, the upcoming theoretical development for the orientation dynamics of a slender prolate spheroid leading up to (4.46) can still be used for a slender cylinder while ignoring the $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ flow and using only the appropriate terms up to $1/\log (\kappa )^3$ instead of the factor $1/(2\log (2\kappa )-3)$ appearing in the following text. Most of the features described in § 5 while analysing the theoretical prediction of the influence of viscoelasticity on the orientation of a slender spheroid will be qualitatively valid for a slender cylinder or a general slender particle.

The Newtonian flow at $O(c)$ is due to a fibre rotating with the perturbed angular velocity, $\boldsymbol {\omega }^{(1)}_p$, in a quiescent fluid, i.e. (3.18a,b)–(3.19). In the outer region, the fluid velocity is

(4.9) \begin{align} {{\boldsymbol{u}_{out}^{M}}^{(1)}}(\boldsymbol{r};\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})& =\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\,\frac{\boldsymbol{h}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_p, \boldsymbol{p})}{\log(2\kappa)-0.5}\boldsymbol{\cdot} \left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}\nonumber\\ &\quad +\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\,\frac{1}{\kappa^2}\boldsymbol{\mathcal{G}}^{(1)}( \lambda;\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p}):\boldsymbol{\nabla} \left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}, \end{align}

where

(4.10a,b)\begin{align} \boldsymbol{h^{(1)}}(\lambda;\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})=4{\rm \pi}\lambda \boldsymbol{\omega}^{(1)}_p\times\boldsymbol{p},\quad \boldsymbol{\mathcal{G}}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})=-2{\rm \pi} (1-\lambda^2)(\boldsymbol{\omega}^{(1)}_p\boldsymbol{\cdot} \boldsymbol{p})\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{p}.\end{align}

The torque generated by this flow is

(4.11)\begin{align} &{\boldsymbol{G}^{MIST}}^{(1)}(\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})\nonumber\\ &\quad =\int_{-1}^1{\rm d}\lambda\left(\frac{1}{\kappa^2}\boldsymbol{\boldsymbol{\epsilon}}:\boldsymbol{\mathcal{G}}^{(1)}- \frac{\lambda\boldsymbol{p}\times\boldsymbol{h}^{(1)}}{\log(2\kappa)-0.5}\right)=\frac{8{\rm \pi}}{3} \left(\frac{(\boldsymbol{\omega}^{(1)}_p\boldsymbol{\cdot}\boldsymbol{p}) \boldsymbol{p}}{\kappa^2}-\frac{\boldsymbol{\omega}^{(1)}_p\boldsymbol{\cdot} (\boldsymbol{I}-\boldsymbol{pp})}{\log(2\kappa)-0.5}\right). \end{align}

Taking the cross product of this equation with orientation vector $\boldsymbol {p}$ leads to the $O(c)$ rotation rate

(4.12)\begin{equation} \dot{\boldsymbol{p}}^{(1)}=\boldsymbol{\omega}^{(1)}_p\times\boldsymbol{p}=-\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5}) ({\boldsymbol{G}^{MIST}}^{(1)}(\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})\times\boldsymbol{p}). \end{equation}

Using the torque balance at $O(c)$ from (3.17) we obtain

(4.13)\begin{equation} \dot{\boldsymbol{p}}^{(1)}=\dot{\boldsymbol{p}}^{(1)}_{PIST}+\dot{\boldsymbol{p}}^{(1)}_{Elastic},\end{equation}

where

(4.14)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}=\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5})({\boldsymbol{G}^{PIST}}^{(1)}\times\boldsymbol{p}) \end{equation}

is the effect of polymer-induced solvent stresses on the particle rotation rate, and

(4.15)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{Elastic}=\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5})({\boldsymbol{G}^{Elastic}}^{(1)} \times\boldsymbol{p})\end{equation}

is the effect of the elastic stress on the particle rotation rate.

4.2. Rotation due to polymer-induced solvent stress

From (2.4), (3.14), (3.16) and (4.14), the rotation rate due to the polymer-induced solvent stress is

(4.16)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}=-\frac{3}{8{\rm \pi} De}({\log(2\kappa)-0.5})\boldsymbol{p}\times\int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}^{(0)}}-{\boldsymbol{\varLambda}^U})) \boldsymbol{\cdot}\boldsymbol{b}.\end{equation}

We remind the reader that here $\boldsymbol {b}$ is the auxiliary ‘velocity’ field used in the reciprocal theorem for deriving the polymer-induced solvent torque ((2.24) or (3.16)). Here $\boldsymbol {b}\boldsymbol{\cdot}\boldsymbol {\omega }$ corresponds to the fluid velocity around a fibre rotating with an angular velocity $\boldsymbol {\omega }$ in a quiescent fluid. Therefore,

(4.17)\begin{equation} \boldsymbol{b}=\boldsymbol{\nabla}_{\boldsymbol{\omega}^{(1)}_p} {\boldsymbol{u}^{M}}^{(1)}, \end{equation}

where ${\boldsymbol {u}^{M}}^{(1)}$ is the solution to (3.18a,b) and (3.19). The volume integral in (4.16) is approximated from the outer region of slender body theory, where the particle is replaced with a line of point forces and doublets. We find the integral from the outer region to converge, and the contribution from the inner region is expected to be small due to the smaller volume of that region. Therefore,

(4.18)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}\approx-\frac{3}{8{\rm \pi} De}({\log(2\kappa)-0.5})\boldsymbol{p}\times\int {\rm d}^3\boldsymbol{r}\,\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}_{out}^{(0)}}-{ \boldsymbol{\varLambda}^U})\boldsymbol{\cdot}\boldsymbol{b}_{out},\end{equation}

where the volume integral is taken over the entire space, ${\boldsymbol {\varLambda }_{out}^{(0)}}$ is the polymer conformation in the outer region and from (4.9), (4.10a,b) and (4.17),

(4.19)\begin{equation} \boldsymbol{b}_{out}=\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\, \frac{4{\rm \pi}\lambda\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\epsilon}}{\log(2\kappa)-0.5} \boldsymbol{\cdot}\left(\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda \boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right)+O(\kappa^{-2}). \end{equation}

It is straightforward to show

(4.20)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}\approx\tfrac{3}{2}(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p}) \boldsymbol{\cdot}\int_{-1}^1{\rm d}\lambda\,\lambda\boldsymbol{u}_{out}^{polymer}(\lambda\boldsymbol{p}), \end{equation}

where

(4.21)\begin{align} \boldsymbol{u}_{out}^{polymer}(\lambda\boldsymbol{p})=\frac{1}{8{\rm \pi}}\int {\rm d}^3\boldsymbol{r}\,\frac{1}{De}\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}_{out}^{(0)}}-{ \boldsymbol{\varLambda}^U})\boldsymbol{\cdot}\left(\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda \boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda\boldsymbol{p})(\boldsymbol{r}-\lambda \boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right) \end{align}

is the velocity field evaluated on the particle centreline that is produced by polymeric force, $({1}/{De})\boldsymbol {\nabla }\boldsymbol{\cdot}({\boldsymbol {\varLambda }_{out}^{(0)}}-{\boldsymbol {\varLambda }^U})$, acting in the outer region. This form of the rotation rate, (4.20) and (4.21), was considered by Harlen & Koch (1993). However, they did not account for all the relevant terms in calculating ${\boldsymbol {\varLambda }_{out}^{(0)}}$ as we show below.

As discussed in § 3 and shown by (3.13) and (4.1), the polymer conformation, ${\boldsymbol {\varLambda }^{(0)}}$ depends on the leading-order Newtonian velocity. In the outer region, the velocity disturbance created by the particle is $O(\max [\,p_2/\log (\kappa ),1/\kappa ^2])$ smaller than the velocity of the imposed simple shear ((4.3), (4.4) and (4.5)). Thus, we linearize (4.1) in the outer region about the imposed flow field, $\boldsymbol {r}\boldsymbol{\cdot}(\boldsymbol {\nabla } \boldsymbol {u})_\infty$, and obtain the governing equation for the disturbance in the polymer conformation from its undisturbed value,

(4.22)\begin{equation} {\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}}=\boldsymbol{\varLambda}_{out}^{(0)}-\boldsymbol{\varLambda}_{out}^U, \end{equation}

to be

(4.23) \begin{align} \left(\boldsymbol{r}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})_\infty\boldsymbol{\cdot}\boldsymbol{\nabla}+\frac{1}{De}\right){\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}}&- (\boldsymbol{\nabla} \boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}}-{\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\nonumber\\ &=(\boldsymbol{\nabla}{\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}}^{\rm T}\boldsymbol{\cdot} \boldsymbol{\varLambda}^U+\boldsymbol{\varLambda}^U\boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}}. \end{align}

It is more convenient to solve (4.23) in Fourier space,

(4.24)\begin{align} \left(-k_1\frac{\partial}{\partial k_2}+\frac{1}{De}\right){\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}}}}-(\boldsymbol{\nabla} \boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}}}} & -{\hat{\boldsymbol{\varLambda}}_ {out}^{(0)^{{\prime}}}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\nonumber\\ &={\rm i}({{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{{\prime}}}}}}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{\varLambda}^U+\boldsymbol{\varLambda}^U\boldsymbol{\cdot}\boldsymbol{k}{{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{{\prime}}}}}}), \end{align}

where ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}=\mathcal {F}({{\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}})$ and ${{\hat {\boldsymbol {u}}_{out}^{M^{(0)^{\prime}}}}}=\mathcal {F}({{{\boldsymbol {u}}_{out}^{M^{(0)^{\prime}}}}})$ are the Fourier transforms of the disturbance of polymer conformation and fluid velocity in the outer region. The rotation rate in (4.20) is expressed as

(4.25)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}\approx\tfrac{3}{2}(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p}) \boldsymbol{\cdot}\int_{-1}^1{\rm d}\lambda\int {\rm d}^3\boldsymbol{k}\,\lambda\exp(2{\rm \pi} {\rm i}\lambda\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p})\hat{\boldsymbol{u}}_{out}^{polymer},\end{equation}

where from the convolution theorem,

(4.26)\begin{equation} \hat{\boldsymbol{u}}_{out}^{polymer}={\rm i}\frac{1}{De}\boldsymbol{k}\boldsymbol{\cdot} {\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}}}}(\boldsymbol{k};\boldsymbol{p})\boldsymbol{\cdot} \hat{\boldsymbol{J}}(\boldsymbol{k}),\quad\text{with } \hat{\boldsymbol{J}}(\boldsymbol{k})=\mathcal{F}\left(\frac{1}{8{\rm \pi}} \left(\frac{\boldsymbol{I}}{r}-\frac{\boldsymbol{rr}}{r^3}\right)\right)= \frac{1}{k^2}\left(\boldsymbol{I}-\frac{\boldsymbol{kk}}{k^2}\right). \end{equation}

We solve the polymer constitutive equation (4.24) in Fourier space using the method of characteristics to obtain ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$. The characteristics are the streamlines of the imposed shear flow (in the Fourier space) and so only integration along the $k_2$ direction is needed. The limits of the integral are $\infty \textrm {sgn}(k_1)$ and $k_2$ and the boundary condition is ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}(\infty \textrm {sgn}(k_1))=0$. We use computer algebra for this calculation. Here ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$, hence obtained, allows us to evaluate the integral in (4.25). The expression used for the disturbance in the Newtonian velocity field, ${\boldsymbol {u}_{out}^{{M}{(0)^{\prime}}}}$, has errors of $O(\textrm {max}[\,p_2^2/\log (\kappa ),p_2/\kappa ^2,1/\kappa ^3])$. Due to the linearization of the constitutive equation, the error in evaluation of rotation rate is $O(\epsilon _{\dot {p}})$, where

(4.27)\begin{equation} \epsilon_{\dot{p}}=\text{max}\left[\frac{p_2^2}{\log(\kappa)},\frac{p_2}{\kappa^2},\frac{1}{\kappa^3}\right].\end{equation}

Due to fortuitous cancelling of terms, a simple expression for the second (gradient direction) component of the rotation rate is obtained,

(4.28)\begin{equation} \dot{{p}}^{(1)}_{2,PIST}=-\frac{De p_1^2 p_2}{2\log(2\kappa)-3}+O(\epsilon_{\dot{p}}),\end{equation}

valid at all $De$. Equivalently, from the outer region integral of (3.16), i.e.

(4.29)\begin{equation} {\boldsymbol{G}^{PIST}}^{(1)}\approx \frac{1}{De}\int {\rm d}^3\boldsymbol{r}\,\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}_{out}^{(0)}}- {\boldsymbol{\varLambda}^U})\boldsymbol{\cdot}\boldsymbol{b}_{out}, \end{equation}

we obtain the first and third components of the polymer-induced solvent torque,

(4.30)\begin{equation} \frac{{{G}_1^{PIST}}^{(1)}}{p_3}=-\frac{{{G}_3^{PIST}}^{(1)}}{p_1}=\frac{4De^2 p_1^2 p_2}{3(\log(2\kappa)-0.5)(\log(2\kappa)-1.5)}+O\left(\frac{\epsilon_{\dot{p}}}{\log(\kappa)}\right). \end{equation}

Since, $\dot {\boldsymbol {p}}^{(1)}_{PIST}={3}/{(8{\rm \pi} )}({\log (2\kappa )-0.5})({\boldsymbol {G}^{PIST}}^{(1)} \times \boldsymbol {p})$ (4.14) we may obtain $\dot {{p}}^{(1)}_{2,PIST}=-De{p_1^2 p_2}/({2\log (2\kappa )-3})$. However, the equations for determining the torque component in the gradient direction, ${{G}_2^{PIST}}^{(1)}$ or rotation rates in flow and vorticity direction, $\dot {{p}}^{(1)}_{1,PIST}$ and $\dot {{p}}^{(3)}_{1,PIST}$ are not tractable for a general $De$. Therefore, to obtain these components we consider small and large $De$ limits separately in §§ 4.2.1 and 4.2.2.

We find the particle rotation rate due to polymer-induced solvent stress in gradient direction, $\dot {{p}}^{(1)}_{2, PIST}$ to be dependent on the first normal stress difference of the Oldroyd-B fluid. This is because $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ in (4.25) is directly proportional to ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$ that is in turn dependent on $\boldsymbol {\varLambda }^U$ via (4.24). This first normal stress difference dependence is further confirmed by repeating the above calculation by artificially setting ${\varLambda }^U_{ij}=De^2\delta _{i1}\delta _{j1}$, such that the only non-zero non-Newtonian property of the fluid is a finite first normal stress difference equivalent to that of an Oldroyd-B fluid. Similarly, we find the two remaining rotation rate components in the limit of large $De$ calculation of § 4.2.1 to arise from the first normal stress difference of the Oldroyd-B fluid. The rotation rate for a second-order fluid ($O(De)$ rotation rate in the small $De$ calculation of § 4.2.2) is also due to the first normal stress difference. An Oldroyd-B fluid has no second normal stress difference. In addition to its appropriate modelling of the simple shear flow of a dilute polymeric liquid, the Oldroyd-B model the advantages over other models such as FENE-P or Giesekus (Bird et al. 1987) of simplicity and hence better analytical tractability. As mentioned in § 1, the second normal stress difference for most polymeric fluids is much smaller than the first normal stress difference. Hence, the effect of the first normal stress difference as ascertained from the Oldroyd-B model is likely to be the most important contribution in determining the influence of polymers on the orientational dynamics of a prolate spheroid in simple shear flow.

4.2.1. Large $De$

In the large $De$ regime, the relaxation of the disturbance in polymer conformation in the outer region is much slower than its convection and stretching by the imposed velocity field. Therefore, when $De\gg 1$, (4.24) simplifies to

(4.31)\begin{align} -k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}-(\boldsymbol{\nabla} \boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}- {\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})_\infty={\rm i}({{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{\prime}}}}} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\varLambda}^U+\boldsymbol{\varLambda}^U\boldsymbol{\cdot} \boldsymbol{k}{{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{\prime}}}}}).\end{align}

This equation is solved in a similar way as (4.24) described earlier. Using ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$ obtained from this equation we find the particle rotation rate due to polymer-induced solvent stresses at large $De$ to be

(4.32)\begin{align} \lim_{De\gg 1}\dot{\boldsymbol{p}}^{(1)}_{PIST}=De\begin{bmatrix} -p_1p_3^2 \alpha_{De\gg 1}\\- \dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\gg 1} \end{bmatrix},\quad \text{with } \alpha_{De\gg 1}=\frac{-4}{\kappa^2}+\frac{3{\rm \pi}}{4\log(2\kappa)-6}p_2\tan^{-1}\frac{p_1}{|p_3|}.\end{align}

The large $De$ approximation can be viewed as the leading $O(De)$ term in an expansion in powers of $1/De$. Thus, the error due to the expansion in $1/De$ is $O(1)$ in the first and third components of $\lim _{De\gg 1}\dot {\boldsymbol {p}}^{(1)}_{PIST}$. There is no error due to expansion in $1/De$ for $\dot {{p}}^{(1)}_{2,PIST}$ as same value is obtained via $De\gg 1$ approximation as that for a general $De$ in (4.28). As mentioned earlier, $\dot {{p}}^{(1)}_{2,PIST}$ has an error of $O(\epsilon _{\dot {p}})$ (4.27). Here $\dot {{p}}^{(1)}_{1,PIST}$ and $\dot {{p}}^{(1)}_{3,PIST}$ have an error of $O(De\epsilon _{\dot {p}})$. To the best of our knowledge, the only previous theoretical study concerning the rotation of a slender particle in a viscoelastic fluid at high $De$ was conducted by Harlen & Koch (1993). However, they did not account for the stretching of the conformation disturbance, ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$, by the mean velocity gradients, i.e. they omitted the $(\boldsymbol {\nabla } \boldsymbol {u})_\infty ^T\boldsymbol{\cdot}{\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}+ {\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}\boldsymbol{\cdot}(\boldsymbol {\nabla } \boldsymbol {u})_\infty$ term in (4.31).

4.2.2. Small $De$

When $De\ll 1$, a solution of (4.23) or (4.24) is obtained via a regular perturbation expansion of ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$ in the powers of $De$

(4.33)\begin{equation} {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}}=\varSigma_{n=0}^\infty De^n {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}^{(n)}}}}}.\end{equation}

The first two terms in this expansion are

(4.34)\begin{align} \left.\begin{gathered} {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}^{(0)}}}}}=\boldsymbol{0},\\ {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime^{(1)}}}}}={\rm i}(\boldsymbol{k}{\hat{\boldsymbol{u}}^{M^{(0)}}} +{\hat{\boldsymbol{u}}^{M^{(0)}}}\boldsymbol{k}). \end{gathered}\right\} \end{align}

The polymer conformation at higher order in $De$ is obtained from the following equations:

(4.35)\begin{equation} \left.\begin{aligned} {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}} & =(\boldsymbol{\nabla}\boldsymbol{u})_\infty^T \boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(1)}}}+{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(1)}}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\\ & \quad +{\rm i}({\hat{\boldsymbol{u}}^{M^{(0)}}}\boldsymbol{k}\boldsymbol{\cdot} {\boldsymbol{\varLambda}^{U^{(1)}}}+{\boldsymbol{\varLambda}^{U^{(1)}}}\boldsymbol{\cdot} \boldsymbol{k}{\hat{\boldsymbol{u}}^{M^{(0)}}})+k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(1)}}},\\ {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(3)}}} & = (\boldsymbol{\nabla}\boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}} + {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\\ & \quad +{\rm i}({\hat{\boldsymbol{u}}^{M^{(0)}}}\boldsymbol{k}\boldsymbol{\cdot}{\boldsymbol{\varLambda}^{U^{(2)}}}+ {\boldsymbol{\varLambda}^{U^{(2)}}}\boldsymbol{\cdot}\boldsymbol{k}{\hat{\boldsymbol{u}}^{M^{(0)}}})+ k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}},\\ {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(n)}}} & =( \boldsymbol{\nabla}\boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(n-1)}}} +{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(n-1)}}}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{u})_\infty+k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{ (n-1)}}},\quad n\ge 4. \end{aligned}\right\}\end{equation}

Therefore, in the limit of low $De$, $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ can be obtained up to arbitrary order in $De$,

(4.36)\begin{equation} \lim_{De\ll 1}\dot{\boldsymbol{p}}^{(1)}_{PIST}=-\frac{3}{8{\rm \pi} }({\log(2\kappa)-0.5})\boldsymbol{p}\times\varSigma_{n=1}^\infty De^{n-1} \int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}({{{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}}^{(n)})) \boldsymbol{\cdot}\boldsymbol{b}. \end{equation}

We will discuss the contribution of the $n=1$ term, i.e. $O(1)$ term in $De$ expansion of $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ in the above equation in § 4.4, along with the similar term in the expansion of $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$. Higher-order contributions to $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ can be evaluated in Fourier space. From (4.33)–(4.35), using the expansion of ${{\boldsymbol {\varLambda }}_{out}^{(0)}}'$ up to $O(De^{8})$ (even higher orders can be easily evaluated), we obtain

(4.37) \begin{align} \lim_{De\ll 1} \dot{\boldsymbol{p}}^{(1)}_{PIST}&=-\frac{3}{8{\rm \pi} }({\log(2\kappa)-0.5}) \boldsymbol{p}\times \int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}({{{\boldsymbol{\varLambda}}_{out}^{(0)}}'}^{(1)})) \boldsymbol{\cdot}\boldsymbol{b}\nonumber\\ &\quad +De\begin{bmatrix} -p_1p_3^2 \alpha_{De\ll 1}\\ -\dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\ll 1} \end{bmatrix}, \end{align}

where

(4.38)\begin{align} \alpha_{De\ll 1}=-\frac{4}{\kappa^2}+\frac{ p_1p_2 }{8\log(2\kappa)-12} \left(5 De-\frac{3}{2}(2+p_3^2)De^3+\frac{13}{10}(8+p_3^2(4+3p_3^2))De^5\right).\end{align}

As expected, we obtain the same expression for $\dot {{p}}^{(1)}_{2,PIST}$ from this perturbation expansion in $De$ valid at small $De$ as was obtained directly for a general $De$ in (4.28). The error in the rotation rate in (4.37) is $O(De\epsilon _{\dot {p}})$ for all the components ($\epsilon _{\dot {p}}$ is given in (4.27)).

4.3. Rotation due to elastic stress

The elastic torque, ${\boldsymbol {G}^{Elastic}}^{(1)}$, is due to the polymer stress on the particle surface ((3.14) and (3.16)). This is evaluated through the polymer constitutive equations on the particle surface written in the frame of reference moving with the particle surface. Due to the absence of polymer convection relative to the particle on the latter's surface, the quasisteady constitutive equation (4.1) simplifies to a set of coupled algebraic equations,

(4.39)\begin{equation} (\boldsymbol{\nabla}{\boldsymbol{u}_{in}^{M}}^{(0)}(\boldsymbol{r}_p; \boldsymbol{\omega}^{(0)}_p,\boldsymbol{p}))^{\rm T}\boldsymbol{\cdot} \boldsymbol{\varLambda}^{(0)}_p+\boldsymbol{\varLambda}^{(0)}_p \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}_{in}^{M}}^{(0)}( \boldsymbol{r}_p;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})+ \frac{1}{De}(\boldsymbol{I}-\boldsymbol{\varLambda}^{(0)}_p)=0, \end{equation}

at each point on the surface. Here $\boldsymbol {r}_p$ is the position vector of a point on the surface. $\boldsymbol {\nabla }{\boldsymbol {u}_{in}^{M}}^{(0)}(\boldsymbol {r}_p;\boldsymbol {\omega }^{(0)}_p,\boldsymbol {p})$ and $\boldsymbol {\varLambda }^{(0)}_p$ represent the surface velocity gradient of the inner velocity field and the surface polymer conformation, respectively, at $\boldsymbol {r}_p$. The inner velocity field is obtained from Cox (1970, 1971). These velocity fields include the effects of point force Stokeslets and doublets that correspond to the outer velocity field accurate up to $O(\,p_2/\log (\kappa )^2)$ and $O(1/\kappa ^2)$. Solving (4.39) via an asymptotic expansion in $1/\kappa$ and ignoring higher-order terms leads to the elastic torque

(4.40)\begin{equation} {\boldsymbol{G}}^{(1)}_{Elastic}=\frac{8{\rm \pi} c}{3\kappa^2}\begin{bmatrix} p_1 p_3\\ 0\\-p_1^2 \end{bmatrix}+O(\max[\,p_2\kappa^{-2},\kappa^{-3}]),\end{equation}

and the corresponding rotation rate

(4.41)\begin{align} \dot{\boldsymbol{p}}^{(1)}_{Elastic}&=-\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5}) \boldsymbol{p}\times\int {\rm d}S\left(\boldsymbol{r}\times\frac{1}{De}(\boldsymbol{\varLambda}^{(0)}_p- \boldsymbol{I})\boldsymbol{\cdot}\boldsymbol{n}\right)\nonumber\\ &=\frac{{\log(2\kappa)-0.5}}{\kappa^2}p_1\begin{bmatrix} 0\\ 1\\0 \end{bmatrix}+O(\log(\kappa)\max[\,p_2\kappa^{-2},\kappa^{-3}]). \end{align}

Therefore, for all $De$, in limits of large $\kappa$ and small $p_2$, the $O(c)$ elastic torque is independent of the polymer relaxation time, $De$.

4.4. Net rotation due to the polymers

Combining the results of the previous two sections, we obtain the net change in the fibre's rotation rate due to polymers.

For $De\ll 1$, we can obtain the $O(1)$ term in the $De$ expansion of $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ and $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$ for all particle orientations, without resorting to the outer region approximation for the former. The first two terms in the polymer configuration, everywhere (inner and outer region) is similar to that mentioned earlier in (4.34), i.e.

(4.42)\begin{equation} {{{{\boldsymbol{\varLambda}}^{(0)}}}=\boldsymbol{I}+De(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)}+(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)})^{\rm T})+O(De^2).}\end{equation}

From the combined effect of polymers from (3.16), (4.13), (4.14) and (4.15), along with using the value of ${{{\boldsymbol {\varLambda }}^{(0)}}}$ from (4.42), divergence theorem, the leading-order momentum equation (i.e. $\boldsymbol {\nabla }{p^{P}}^{(0)}=\boldsymbol {\nabla }\boldsymbol{\cdot}[\boldsymbol {\nabla } {\boldsymbol {u}^{P}}^{(1)}+(\boldsymbol {\nabla } {\boldsymbol {u}^{P}}^{(1)})^\textrm {T}]$) and the auxiliary Stokes problem ((2.18a–c) and (2.19)), we obtain

(4.43) \begin{align} \lim_{De\ll 1}{\dot{\boldsymbol{p}}^{(1)}}_{Elastic}&+{\dot{\boldsymbol{p}}^{(1)}}_{PIST}\nonumber\\ &=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times[(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)}+(\boldsymbol{\nabla}{\boldsymbol{u}^{M}}^{(0)})^{\rm T}-{p^{M}}^{(0)} \boldsymbol{I})\boldsymbol{\cdot}\boldsymbol{n}]+O(De)\nonumber\\ &= \int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times {\boldsymbol{\tau}^{M}}^{(0)}+O(De)\boldsymbol{\cdot}\boldsymbol{n}=O(De). \end{align}

The $O(1)$ term in this sum is equivalent to the rotation due to the leading-order torque (in $c$) acting on the particle. This is the torque on a freely rotating particle in a Newtonian fluid and is hence zero. Thus, for a freely rotating particle for $De\ll 1$, the $O(1)$ rotation rates (in the $De$ expansion) due to the torque generated by the elastic and polymer induced solvent stress cancel each other and lead to no change in the net particle rotation. The change in the particle rotation due to the polymers arises at $O(De)$. Close to the flow–vorticity plane, the rotation rate due to the elastic torque was shown to be independent of $De$ for all $De$ in (4.41), which we have just shown balances the equivalent rotation rate due to polymer-induced solvent torque for $De\ll 1$. Thus, the $O(De)$ and higher effect of polymers on the particle's rotation rate arises entirely due to the polymer-induced solvent torque (4.37) and is given by

(4.44) \begin{equation} \lim_{De\ll 1}(\kern0.7pt\dot{\boldsymbol{p}}^{(1)}_{Elastic}+\dot{\boldsymbol{p}}^{(1)}_{PIST})=De\begin{bmatrix} -p_1p_3^2 \alpha_{De\ll 1}\\ -\dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\ll 1} \end{bmatrix},\end{equation}

where $\alpha _{De\ll 1}$ is given in (4.38) and the errors are $O(De\epsilon _{\dot {p}})$ for all the components (where $\epsilon _{\dot {p}}$ is given in (4.27)).

For $De\gg 1$, the total effect of viscoelasticity on the particle rotation also arises mainly from the polymer-induced solvent stress. But, unlike the $De\ll 1$ regime, here the reason is that the effect of the elastic stress is $O(De)$ smaller. The net rotation rate due to the polymers is

(4.45)\begin{equation} \lim_{De\gg 1}(\kern0.7pt\dot{\boldsymbol{p}}^{(1)}_{Elastic}+\dot{\boldsymbol{p}}^{(1)}_{PIST})=De\begin{bmatrix} -p_1p_3^2 \alpha_{De\gg 1}\\ -\dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\gg 1} \end{bmatrix},\end{equation}

where $\alpha _{De\gg 1}$ is given by (4.32). The error in $\dot {{p}}^{(1)}_{1,PIST}$ and $\dot {{p}}^{(1)}_{3,PIST}$ is of $O(De\epsilon _{\dot {p}})$ (4.27). Neglecting the elastic torque leads to an additional error of $O(\log (\kappa )/\kappa ^2)$ in the second component and $O(\log (\kappa )\max [\,p_2\kappa ^{-2},\kappa ^{-3}])$ in the first and third component of the net rotation rate due to the polymers (4.41).

4.5. Equations of motion of a freely rotating fibre in low $c$ viscoelastic fluid at all $De$

Equations (4.44) and (4.45) encompass our main result for the effect of viscoelasticity on the rotation of a particle suspended in simple shear flow at large and small $De$, respectively. In the limit of large $De$ and for a second-order fluid ($O(De)$ rotation rate in small $De$ limit), the effect of viscoelasticity is due to the first normal stress difference of the fluid as discussed in § 4.2. The governing equation for the orientation dynamics of a slender prolate spheroid that includes the viscoelastic effects near the flow–vorticity plane is

(4.46)\begin{equation} \dot{\boldsymbol{p}}=\frac{{\rm d}\boldsymbol{p}}{{\rm d}t}=\begin{bmatrix} \dfrac{p_2}{2}+\dfrac{\kappa^2-1}{\kappa^2+1}\left(\dfrac{p_2}{2}-p_1^2p_2\right)-c\,Dep_1p_3^2 \alpha\\ - \dfrac{p_1}{2}+\dfrac{\kappa^2-1}{\kappa^2+1}\left(\dfrac{p_1}{2}-p_1p_2^2\right)-c\,De \dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ -\dfrac{\kappa^2-1}{\kappa^2+1}p_1 p_2 p_3+c\,Dep_1^2p_3\alpha \end{bmatrix},\end{equation}

where in the large $De$ limit $\alpha$ is $\alpha _{De\gg 1}$ (4.32). In the small $De$ limit we consider $\alpha _{De\ll 1}$ (4.38) up to $O(De)$ and interpolating between these two limits of $De$, we obtain the following uniformly valid approximation for $\alpha$ at all $De$:

(4.47)\begin{equation} \alpha=-\frac{4}{\kappa^2}+\frac{p_2}{4\log(2\kappa)-6}\frac{2.5p_1De}{1+2.5p_1De\left/\left(3{\rm \pi}\tan^{-1}\dfrac{p_1}{|p_3|}\right)\right.}.\end{equation}

Equation (4.46) introduces no errors in the Newtonian rotation rate of a prolate spheroidal particle. The slender body theory provides a good approximation for the Newtonian velocity field for $\kappa \gtrsim 10$. The errors in these equations due to the viscoelastic terms are of $O(c\,De\epsilon _{\dot {p}})$ for all the components ($\epsilon _{\dot {p}}$ is given in (4.27)) in the $De\ll 1$ limit. For $De\gg 1$, the errors are of $O(c\,De\epsilon _{\dot {p}})$ in the first and third component and $O(c\max [\log (\kappa )/\kappa ^2,\epsilon _{\dot {p}}])$ in the second component. Before analysing the influence of polymers on particle orientation as suggested by this equation we compare our theory at low $De$ with that of Leal (1975).

4.6. Comparison of second-order fluid result with Leal (1975)

Leal (1975) considered the motion of a fibre in a second-order fluid and found

(4.48)$$\begin{gather} \dot{p}_2^{Leal}=- p_1p_2 (\,p_2+Vp_1 (1-2p_2^2)) \end{gather}$$

(4.49)$$\begin{gather}\dot{p}_3^{Leal}=- p_1 p_2 p_3(1-2V p_1p_3), \end{gather}$$

where $V=-({3 \lambda \gamma }/{16 \log (\kappa )})M_1(1+2\epsilon _1)$, $\lambda =\varPhi _3 U/\mu l$ is the polymer relaxation time ($\mu$ is the zero-shear-rate viscosity, $l$ the particle half-length and $U$ a characteristic velocity scale) and $\epsilon _1=\varPhi _2/\varPhi _3$,

(4.50a,b)\begin{equation} \varPhi_2=-\lim_{\gamma\rightarrow 0}\frac{\sigma_{11}-\sigma_{22}}{2\gamma^2} \quad\text{and}\quad \varPhi_3=\lim_{\gamma\rightarrow 0}\frac{\sigma_{11}-\sigma_{33}}{\gamma^2}, \end{equation}

where $\sigma _{11}$, $\sigma _{22}$, $\sigma _{33}$ and $\gamma$ are the first, second and third normal stresses and the shear rate. and $M_1$ is a positive number depending upon the shape of the particle. In non-dimensional terms $\lambda =2De$, $V=({-3De}/{8\log (\kappa )})M_1(1+2\epsilon _1)$. Here $(1+2\epsilon _1)$, and hence $V$ in Leal's theory, is proportional to the second normal stress difference. Up to $O(De p_2)$ Leal's results are, therefore,

(4.51)$$\begin{gather} \dot{p}_2^{Leal}=-p_1p_2^2+\frac{3De}{8\log(\kappa)}M_1 (1+2\epsilon_1)p_1^2p_2 +O(De^2,p_2^2De), \end{gather}$$

(4.52)$$\begin{gather}\dot{p}_3^{Leal}=-p_1p_2p_3+O(De^2,p_2^2De). \end{gather}$$

By taking only up to $O(De)$ terms from (4.46), our results up to this order for a prolate spheroidal particle in a second-order fluid are

(4.53)$$\begin{gather} \dot{p}_2=-p_1p_2^2-\frac{c\,De}{2\log(2\kappa)-3}p_1^2p_2 +O(c\,De^2,cp_2^2De), \end{gather}$$

(4.54)$$\begin{gather}\dot{p}_3=-p_1p_2p_3+O(c\,De^2,cp_2^2De). \end{gather}$$

Hence, our theory's first viscoelastic effects at low $De$ have the same functional dependence of rotation rates on the particle orientation as that of Leal's second-order fluid theory. However, our theory predicts these effects to arise from the first normal stress difference. In contrast, Leal's theory predicts these to emerge from the second normal stress difference. According to our theory, a prolate spheroid rotating in a simple shear flow of a Boger fluid that has zero second normal stress difference (Magda et al. 1991), will exhibit a different rotational motion at a finite but small $c\, De$ as compared with $c=0$ (Newtonian fluid). Specifically, a particle spirals towards the stable limit cycle close to the vorticity, as discussed in the next section. However, Leal's theory predicts Jeffery rotations or no effect of viscoelasticity. Brunn (1977) considered the motion of rigid particles in second-order fluid. While Brunn (1977) does not calculate the rotation rates for a rod-like particle such as prolate spheroid, for a dumbbell representing two spheres joined by a thin, rigid rod, he also obtains the same functional dependence of rotations rates on the particle orientation. However, instead of the factor $M_1(1+2\epsilon _1)$, Brunn (1977) obtains $(1+4\epsilon _1)$ so that the rotation rate of a dumbbell in a second-order fluid with zero second normal stress difference is finite.

5. Analysis of particle orientation dynamics

Before a more detailed analytical and numerical treatment of (4.46), we qualitatively describe the changes in orientation dynamics introduced by viscoelasticity as given by this equation. We obtain three primary regions in $c\, De- \kappa$ space where viscoelasticity leads to different final fibre orientation behaviours. These regions, $R_{vort}=R_{vort}^{(1)}\cup R_{vort}^{(2)}$, $R_{flow}=R_{flow}^{(1)}\cup R_{flow}^{(2)}$ and $R_{flow-vort}=R_{flow-vort}^{(1)}\cup R_{flow-vort}^{(2)}\cup R_{flow-vort}^{(3)}$, separated by the curves $\lambda _{2,flow}^-=0$ and $c\, De=c\, De|_{cut-off}^{Limit}$ are shown in figure 2 for $c=0.005$. In $R_{vort}$, i.e. the region to the left of the curve $\lambda _{2, flow}^-=0$ in figure 2, irrespective of the initial orientation, the particle is attracted towards a stable limit cycle near the vorticity axis where the particle revolves in a small periodic orbit around the vorticity axis. In $R_{flow}$, i.e. the region to the right of the curve $c\, De=c\, De|_{cut-off}^{Limit}$ in figure 2, a particle aligns close to the flow direction for all initial orientations as a stable fixed point at this location is the only stable attractor in the orientation space. This flow alignment behaviour is similar to the particle trajectories observed by Iso et al. (1996a) at large elasticity (large $c$). In $R_{flow-vort}$, defined as the region between $\lambda _{2,flow}^-=0$ and $c\, De=c\, De|_{cut-off}^{Limit}$ in figure 2, depending upon the initial orientation, the particle can either obtain a final orientation within the FGP (it may either obtain a stable orientation or rotate within the plane) or rotate in a small periodic orbit around the vorticity axis.

Figure 2. The $\kappa$–$c\, De$ parameter space, for $c=0.005$, divided into regions with different qualitative behaviour of a particle's orientation dynamics. Here $b^2$ and $\zeta$ are given in (5.4) and (5.12a,b), respectively; $\lambda _{2,{flow}}^{0-}$ and $\lambda _{2,{flow}}^{0--}$ are in (5.14). The procedure for numerically obtaining $c\, De|_{cut\textrm -off}^{Limit}$ is described in § 5.4.1.

There are subdivisions of the regions $R_{vort}$, $R_{flow-vort}$ and $R_{flow}$ based on the behaviour of the particle's orientation trajectory near the FGP and the vorticity axis. In $R_{vort}^{(1)}$ (figure 2), a particle starting near the FGP spirals away from the plane and towards the limit cycle near the vorticity axis. In contrast, in $R_{vort}^{(2)}$ (figure 2), once a particle comes close to the flow direction, it drifts along the flow–vorticity plane in a monotonic fashion. $R_{vort}^{(1)}$ trajectories are reminiscent of the observations of Gauthier et al. (1971) at low shear rates (low $De$) and also some of the particle trajectories of Iso et al. (1996a,b) at low to medium elasticity (small to medium $c$). The trajectories in $R_{vort}^{(2)}$, on the other hand, are similar to that observed by Bartram et al. (1975) at larger shear rates or $De$ (than the earlier low shear rate experiments with the same fluid reported by the same laboratory in Gauthier et al. (1971)). Within the region $R_{flow}$ shown in figure 2, trajectories in $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ have different behaviour near the vorticity axis. In the former, to the left of the curve $\kappa =c\, De$ in figure 2, the vorticity axis is an unstable spiral. Thus, the particle leaves the region close to the vorticity axis in a spiral motion. In the latter, the vorticity axis is an unstable node, and the particle leaves the vicinity of the vorticity axis in a monotonic fashion. The trajectories near the vorticity axis in also in $R_{vort}$ and $R_{flow-vort}$ are similar to those in $R_{flow}^{(1)}$. As mentioned above, in the region of multiple final particle orientations, $R_{flow-vort}$, in addition to a stable limit cycle/periodic orbit close to the vorticity axis, a stable attractor exists in the FGP. $R_{flow-vort}$ can be further subdivided into three subregions. In $R_{flow-vort}^{(1)}$ (figure 2), the particle spirals towards the FGP and then tumbles within the plane. In this case, the particle slows down close to the flow direction before speeding up again as it departs from this orientation. In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ (figure 2), there is a stable fixed point very close to the flow direction. Hence a particle starting close to the FGP ends up being flow-aligned. In $R_{flow-vort}^{(2)}$, there is an unstable node within the FGP (farther away from the flow direction than the stable fixed point). In contrast, in $R_{flow-vort}^{(3)}$, the unstable node is replaced with a saddle point with its stable manifold perpendicular to the FGP. Therefore, the basin of attraction of the stable fixed point (with the stable limit cycle being the other stable attractor) is larger for $R_{flow-vort}^{(3)}$ than for $R_{flow-vort}^{(2)}$.

The above is a qualitative description of all the different behaviours exhibited by the dynamical system of (4.46). However, the predictions for regions of large $c\, De$ (such as $R_{flow}^{(2)}$) are speculative. This is because the first normal stress difference in the polymer stress is proportional to $c\, De$ and our regular perturbation theory is based on the polymer stress being smaller than the Newtonian stress. We nevertheless include these regions here to provide a complete description of the dynamical system and later numerical studies of the original governing equations described in § 2 may be useful in determining the quantitative and qualitative validation of the theory considered here.

In the rest of this section we will derive the boundaries that determine the aforementioned divisions of $c\, De- \kappa$ space and provide a more detailed analytical and numerical treatment of (4.46). The $b^2=0$ and $\kappa =2c\, De$ boundaries shown in figure 2 will be determined in §§ 5.1 and 5.2, respectively. The $\lambda _{2,{flow}}^{0-}=0$, $\lambda _{2,{flow}}^{0--}=0$ and $\zeta =0$ boundaries are derived in § 5.3. Here $c\, De|_{cut\textrm -off}^{Limit}$, is numerically obtained (the other boundaries of figure 2 are analytical) in § 5.4.1. To numerically integrate the orientation trajectory we transform equation (4.46) into $\theta -\phi$ coordinates, where

(5.1)\begin{equation} \boldsymbol{p}=\begin{bmatrix} p_1\\p_2\\p_3 \end{bmatrix}=\begin{bmatrix} \sin(\theta)\sin(\phi)\\ \sin(\theta)\cos(\phi)\\ \cos(\theta) \end{bmatrix}\rightarrow\begin{array}{l} \dfrac{{\rm d}\theta}{{\rm d} t}=-\dfrac{1}{\sin(\theta)}\dfrac{\partial p_3}{\partial t}, \\ \dfrac{{\rm d}\phi}{{\rm d}t}=\dfrac{1}{\sin^2(\theta)}\left(\dfrac{\partial p_2}{\partial t}\cos(\phi)- \dfrac{\partial p_1}{\partial t}\sin(\phi)\right) \end{array}. \end{equation}

Integrating this $\theta -\phi$ system instead of the equivalent equation (4.46) directly for $\boldsymbol {p}$ numerically preserves $||\boldsymbol {p}||_2=1$ constraint.

The boundaries in the $c\, De-\kappa$ space shown in figure 2 are obtained for $c=0.005$. The primary $c$ dependence of the various boundaries is in the form of $c\, De$ and similar partitioning of the $c\, De- \kappa$ space is found with different values of $c$, albeit with small quantitative changes. The boundaries associated with $b^2=0$ (5.4) and $\kappa =2c\,De$ only depend on $c\, De$ for a given $\kappa$. We find (§ 5.4.1) the $c\, De|_{cut\textrm -off}^{Limit}$ boundary to be insensitive to changes in $c$ at constant $c\, De$ for $c\lessapprox 0.1$. However, increasing $c$ moves the $\lambda _{2,{flow}}^{0-}=0$, $\lambda _{2,{flow}}^{0--}=0$ and $\zeta =0$ boundaries in the $c\, De-\kappa$ kappa space to the right or larger $c\, De$ (not shown) even at small $c$. The $\lambda _{2,{flow}}^{0--}=0$ curve moves more than $\lambda _{2,{flow}}^{0-}=0$ upon increasing $c$. Thus, in the $c\, De-\kappa$ space, upon increasing $c$, $R_{vort}^{(1)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ are unchanged, $R_{vort}^{(2)}$ and $R_{flow-vort}^{(2)}$ are enlarged, and, $R_{flow-vort}^{(1)}$ and $R_{flow-vort}^{(3)}$ are reduced in size.

5.1. Low dimensional orbits: log-rolling and tumbling

We begin our analysis with two convenient orientational states: (a) log-rolling with the particle aligned with the vorticity axis, and, (b) tumbling in the FGP ($p_3=0$). Due to symmetry, and as suggested by (4.46), viscoelasticity does not change the particle's orientation from the log-rolling state. Also, in the log-rolling state the theory predicts that the polymer-induced solvent and elastic torques due to viscoelasticity, i.e. ${{G}_3^{PIST}}^{(1)}$ from (4.30) and ${{G}_3^{Elastic}}^{(1)}$ from (4.40), are both zero. Thus, the polymers do not change the log-rolling angular velocity of the particle at $O (c)$.

In the FGP, close to the flow direction $p_1\approx 1$, the equation governing the particle orientation is

(5.2)\begin{equation} \dot{p}_2\approx-p_2^2-\frac{1}{\kappa^2}-c\,De\frac{p_2}{2\log(2\kappa)-3}.\end{equation}

Therefore, we notice that viscoelasticity ($c\, De>0$) reduces the rotation rate of the particle as compared with the Newtonian value ($c\, De=0$). This equation has an analytical solution

(5.3)\begin{equation} p_2=-b\tan({(t-t_0)b})\frac{c\,De}{4\log(2\kappa)-6},\end{equation}

where

(5.4)\begin{equation} b^2={\frac{1}{\kappa^2}-\left(\frac{c\,De}{4\log(2\kappa)-6}\right)^2}.\end{equation}

At small values of $c\, De$ when $b^2>0$, the solution is periodic with a time-period $2{\rm \pi} /b$. Upon increasing $c\, De$, $b=0$ when $c\,De_{crit}=(4\log (2\kappa )-6)/\kappa$ and an infinite period bifurcation occurs. Further increasing $c\, De$ lead to $b^2<0$ and two fixed points appear at $[\,p_{2,{flow}}^{0-},0]$ and $[\,p_{2,{flow}}^{0--},0]$ within the FGP, where

(5.5a,b)\begin{equation} p_{2,{flow}}^{0-}=-\frac{c\,De}{4\log(2\kappa)-6}+\sqrt{-b^2},\quad p_{2,{flow}}^{0--}=-\frac{c\,De}{4\log(2\kappa)-6}-\sqrt{-b^2}.\end{equation}

There are two additional fixed points at $-p_{2,{flow}}^{0-}$ and $-p_{2,{flow}}^{0--}$ that can be mathematically obtained by repeating the above analysis near $p_1=-1$. In the sense of Jeffery orbits, the fixed points are downstream of the flow–vorticity plane and very close to the flow direction when $\kappa$ is large, and $c\, De$ is small; $p_{2,{flow}}^{0--}$ is farther downstream. Within the FGP, particle trajectories approach $\pm p_{2,{flow}}^{0-}$, while they depart $\pm p_{2,{flow}}^{0--}$. Therefore, in the $b^2<0$ regime, a particle placed in the FGP approaches a steady state orientation $p_2=\pm p_{2,{flow}}^{0-}$. We will later observe in § 5.3 that in the orientation space near the flow-direction off the FGP, i.e. for a finite $0< p_3\ll 1$, the trajectories may either approach or leave these two fixed points along the vorticity direction.

5.1.1. More accurate location of fixed points on FGP in $b^2<0$ regime

Above, we analysed the equations near the flow direction under the assumption $p_1\approx 1$ and found two fixed points in the $b^2<0$ regime at orientations given by (5.5a,b). These expressions are only valid when the fixed points are near the flow direction. However, as $c\, De$ is increased at a given $\kappa$, the fixed points separate from each other in the flow gradient plane. The fixed point at $p_{2,{flow}}^-$ moves closer to the flow direction, while $p_{2,{flow}}^{--}$ moves away from the flow direction. To better estimate the latter's location we relax the assumption of $p_1\approx 1$ in the expression for ${\textrm {d}p_2}/{\textrm {d}t}$ and obtain the improved expressions

(5.6) \begin{align} \left.\begin{gathered} \tilde{p}_{2,{flow}}^{0-}=-\sqrt{\frac{-f^2/\kappa^2+{c^2De^2}/2-fc \, De\sqrt{-\tilde{b}^2}}{{c^2De^2+f^2}}},\\ \tilde{p}_{2,{flow}}^{0--}=-\sqrt{\frac{-f^2/\kappa^2+{c^2De^2}/2+fc\,De\sqrt{-\tilde{b}^2}}{{c^2De^2+f^2}}},\end{gathered}\right\} \end{align}

where $\tilde {b}^2={{(1+\kappa ^2)}/{\kappa ^4}-({c\,De}/{(2f)})^2}\approx b^2$ and $f=2\log (2\kappa )-3$. From figure 3, as $c\, De$ or $\kappa$ are varied, we observe the qualitative behaviours of $\tilde {p}_{2,{flow}}^{0-}$ and $\tilde {p}_{2,{flow}}^{0--}$ are the same as those of ${p}_{2,{flow}}^{0-}$ and ${p}_{2,{flow}}^{0--}$ mentioned above. Here $\tilde {p}_{2,{flow}}^{0-}$ is shown with dashed lines and $\tilde {p}_{2,{flow}}^{0--}$ with solid lines in the region $\widetilde {b^2}\approx {b^2}<0$. At the beginning of the dashed and solid lines, $\widetilde {b^2}\approx {b^2}=0$, for a given $\kappa$, $\tilde {p}_{2,{flow}}^{0-}=\tilde {p}_{2,{flow}}^{0--}$, because these fixed points arise out of an infinite period bifurcation (Strogatz 2018), as mentioned earlier.

Figure 3. Variation of the fixed points’ location on the FGP with $c\, De$ for various $\kappa$ in the $\widetilde {b^2}\approx {b^2}<0$ regime.

Therefore, while polymers have no influence on the log-rolling motion, they slow down or stop the tumbling motion of a prolate spheroidal particle within the FGP. We now move on to consider three-dimensional orbits by first observing the particle's orientation trajectories close to the vorticity axis in § 5.2 and then the trajectories near the FGP in § 5.3.

5.2. Effect of viscoelasticity near the vorticity direction, $p_3\approx 1$

Equation (4.46) has a fixed point on the vorticity axis,

(5.7)\begin{equation} \boldsymbol{p}^0_{vort}=\left[0\quad 0\quad 1 \right]. \end{equation}

Due to the constraint $||\boldsymbol {p}||_2=1$, we consider the linear stability of the $p_1-p_2$ dynamical system using $p_3=\sqrt {1-p_1^2-p_2^2}$. In the $p_1-p_2$ coordinate system, at the fixed point $[0\ 0]$, the eigenvalues are

(5.8a,b)\begin{align} \lambda_{1,{vort}}=2c\,De/\kappa^2-1/\kappa\sqrt{4c^2De^2/\kappa^2-1},\quad \lambda_{2,{vort}}=2c\,De/\kappa^2+1/\kappa\sqrt{4c^2De^2/\kappa^2-1},\end{align}

with the corresponding eigenvectors

(5.9a,b)\begin{equation} \boldsymbol{v}_{1,{vort}}=[ \lambda_{1,{vort}}\quad 1 ],\quad \boldsymbol{v}_{2,{vort}}=[ 1\quad \lambda_{2,{vort}} ]. \end{equation}

The fixed point at the vorticity axis undergoes a Hopf bifurcation at $c\, De=0$ and is unstable for all finite $c\, De$. The vorticity axis is a centre (imaginary eigenvalues) when $c\, De=0$. In the presence of viscoelasticity, $c\, De\ne 0$, it becomes a hyperbolic fixed point (eigenvalues with non-zero real part). Therefore, the linear stability analysis provides qualitative insight into the effect of viscoelasticity on the behaviour of fibre orientation governed by the full nonlinear equation (4.46) near the vorticity axis when $c\, De\ne 0$. This follows from the Hartman–Grobman theorem (see Guckenheimer & Holmes 2013), which guarantees the homeomorphism between the linearized and full nonlinear system near hyperbolic fixed points, while also preserving the time parametrization. The instability of the vorticity axis arises from the polymer conformation driven by the force doublet and $O(1/\kappa ^2)$ Stokeslet discussed in § 4.1. When $0< c\,De<\kappa /2$ a small perturbation leads to a particle departing the vorticity axis in a spiral fashion (as the vorticity fixed point is an unstable spiral). However, for $c\,De>\kappa /2$ a particle departs the vorticity axis monotonically (as the vorticity fixed point is an unstable fixed point).

The linear stability analysis is confirmed by the full numerical integration of (4.46) for $\kappa =50$ at $c\, De=24$ (also for $c\, De=5$ and 6) versus $c\, De=26$ in figure 4(a) zoomed near the vorticity axis. At $\kappa =50$, $c\, De=26$ is a point in $R_{flow}^{(2)}$, $c\, De=24$ and 6 are points in $R_{flow}^{(1)}$ and $c\, De=5$ is a point in $R_{flow-vort}^{(3)}$. In $R_{flow}^{(2)}$, the particle drifts out of the vorticity axis monotonically, and in the rest, it spirals out. We will discuss the dynamical system's features of figure 4(b) in more detail later in § 5.4.1 where we will find the unstable spiral at the vorticity axes to be surrounded by a stable limit cycle up to a cut-off $c\, De$. Hence, the Hopf bifurcation occurring at $c\, De=0$ at the vorticity axis is supercritical (Strogatz 2018). In $R_{flow-vort}^{(3)}$ ($c\, De=5$), the particle starting very close to the vorticity axis spirals outwards into the stable limit cycle, where it then moves periodically around the vorticity axis. The limit cycle does not exist in $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$, so that, after the particle either spirals or monotonically drifts out of the vorticity axis, it ends up in the flow aligned state as it approaches the stable fixed point ($\pm \tilde {p}_{2,{flow}}^{0-}$) near the flow direction in the FGP for $c\, De=6$, 24 and 25.

Figure 4. Various orientation behaviours near the vorticity axis: trajectories of particle orientation starting very close to the vorticity axis at different $c\, De$ in $R_{flow}^{(1)}$ ($c\, De=6$, 24), $R_{flow}^{(2)}$ ($c\, De=26$) and $R_{flow-vort}^{(3)}$ ($c\, De=5$) at $\kappa =50$. All trajectories start at the same point. Panel (a) is the same as panel (b) (showing complete particle trajectory) but zoomed near the vorticity axis ($p_3=1$). The grey surface is the unit sphere, i.e. the orientation space.

The identification of the stable limit cycle near the vorticity axis is made possible by including the flow generated by force dipoles per unit length (flow from Cox (1971)) along the fibre in addition to the flow generated by the force per unit length (flow from Batchelor (1970)). If this dipole-generated flow is neglected from our theory the stable limit cycle will not be predicted and instead when the theory predicts particle to be repelled from the FGP it will approach the vorticity axis. A previous slender body theory for a second-order fluid by Férec et al. (2021), relying on only the flow generated by the force per unit length, predicts the fibre to align with the vorticity axis. When the axis of symmetry of an axisymmetric slender fibre is aligned in the flow–vorticity plane of the imposed flow no force per unit length is exerted by the fluid (Newtonian or polymeric). This is because in the flow–vorticity-aligned state the imposed flow has no variation along the fibre axis. However, according to (29) and (35) of Férec et al. (2021) when a fibre is in the flow–vorticity plane, the force per unit length is non-zero and proportional to the polymer relaxation time. This implies an error in the expression for the tension force in (29) of Férec et al. (2021).

Figure 3(a) of Wang, Yu & Lin (2019) clearly shows that a $\kappa =4.0$ spheroid started away from the vorticity axis in a plane Couette flow at $De=0.1$ reaches a closed orbit around the vorticity axis, instead of approaching the axis. Figures 2(a) and 2(b) of d'Avino et al. (2014) at $De=1.0$ and 2.0, respectively, for a $\kappa =4.0$ spheroid in an unbounded simple shear flow show drift towards the vorticity axis, but do not show the particle approaching the axis. Although these studies were conducted at high polymer concentration, $c$, and small $\kappa$, outside the formal range of validity of our theory, they indicate the possibility of a stable limit cycle around the vorticity axis. Furthermore, in the second-order fluid regime with $De=0.1$, careful observation of figure 8(a) of the theoretical investigation of Wang, Tai & Narsimhan (2020) shows a $\kappa =3.0$ spheroid approaching a stable limit cycle, instead of reaching the vorticity axis of an unbounded parabolic slit flow ($u=1-y^2$). Lastly, while the boundary element formulation aided study of Phan-Thien & Fan (2002) in an Oldroyd-B fluid shows the $\kappa =2.0$ particle drifting towards the vorticity axis, the simulation (see figures 7 and 8 of their paper) stops before we can conclude if it will approach the axis or a stable limit cycle.

5.3. Effect of viscoelasticity near the FGP and flow direction

We described the effect of viscoelasticity on the particle's motion near the flow direction within the FGP in § 5.1. Here, we consider the particle motion near (but not exactly on) the FGP. The analysis of § 5.1 for the $p_2$ (gradient) direction is valid even outside the FGP. Near the flow direction, $p_1\approx 1$, the orientation dynamics and solution in the $p_2$ direction are given by (5.2) and (5.3). The orientation dynamics in the $p_3$ (vorticity) direction for $p_1\approx 1$ is governed by the simplified equation (from (4.46))

(5.10)\begin{equation} \dot{p}_3\approx-p_2p_3+c\,Dep_3\left(-\frac{4}{\kappa^2}+\frac{p_2}{8\log(2\kappa)-12}\frac{15{\rm \pi}^2De}{3{\rm \pi}^2+5De}\right)+O(\,p_3^2).\end{equation}

Its closed form solution is

(5.11)\begin{equation} p_3=C\exp(c\,De\zeta (t-t_0))\cos({(t-t_0)b})^{-1+c\,De\gamma/(8\log(2\kappa)-12)}\end{equation}

where $b^2$ is given in (5.4), and

(5.12a,b)\begin{equation} \gamma=\frac{15{\rm \pi}^2De}{3{\rm \pi}^2+5De},\quad \zeta=\frac{1}{4\log(2\kappa)-6}-c\,De\frac{\gamma}{2(4\log(2\kappa)-6)^2}-\frac{4}{\kappa^2}.\end{equation}

When $b$ is real-valued, i.e. $b^2>0$ (5.4), there are no fixed points on the FGP. Thus, the particle undergoes periodic motion in the FGP, as discussed in § 5.1. For a particle perturbed slightly away from the FGP, from (5.11) we note that the sign of $\zeta$ determines whether the particle will be attracted to or repelled from the FGP ($p_3=0$). Here $\zeta >0$ represents the region $R_{vort}^{(1)}$ where the particle spirals away from the FGP and $\zeta <0$ represents the region $R_{flow-vort}^{(1)}$ where the particle spirals into the FGP and continues tumbling within the plane (albeit with a larger orbit period $2{\rm \pi} /b$ than in Newtonian case). In the $b^2>0$ regime, the phase portrait of the dynamical system of (4.46) (through the simplified equations (5.2) and (5.10)) projected in the $p_2\unicode{x2013} p_3$ plane near $p_1=1$ is shown in figure 5 for $\zeta >0$ and $\zeta <0$. The phase flow in the $p_2\unicode{x2013} p_3$ plane in the case of a Newtonian fluid ($c\, De=0$) is similar to that in the region $R_{vort}^{(1)}$ (figure 5a), but it is symmetric about $p_2=0$ (not shown). Therefore, a particle in a Newtonian fluid continues in a particular (periodic) Jeffery orbit before and after $p_2=0$. At small but finite $c\, De$, i.e. in the region $R_{vort}^{(1)}$, represented figure 5(a), this symmetry about $p_2=0$ is broken and a particle comes out of the $p_2=0$ plane with a larger $\dot {p}_3$ velocity than it enters the plane. This explains the drift towards the vorticity axis (greater $p_3$). In the region $R_{flow-vort}^{(1)}$ represented by figure 5(b) we can observe that the phase flow points towards $p_3=0$, indicating migration of a particle towards the FGP. The FGP is a stable limit cycle in $R_{flow-vort}^{(1)}$ and an unstable limit cycle in $R_{vort}^{(1)}$.

Figure 5. Phase portraits of (5.2) and (5.10) in the $b^2>0$ regime in the gradient ($p_2$)–vorticity ($p_3$) plane with $\zeta >0$ (a) and $\zeta <0$ (b). When $b^2>0$ and $\zeta >0$, i.e. in (a) representing $R_{vort}^{(1)}$ close to the flow gradient plane ($p_3=0$), the $p_3$ component of the phase velocity changes sign from negative to positive along with an increase in magnitude downstream of the $p_2=0$ plane indicating a departure of a particle from the flow gradient plane at a rate higher than it approaches the plane. However, when $b^2>0$ and $\zeta <0$, i.e. in (b) representing $R_{flow-vort}^{(1)}$ the $p_3$ component of the phase velocity is negative for all $p_2$ indicating an approach towards the flow gradient plane. Since $p_2$ never approaches zero on the flow gradient plane for $b^2>0$ it is an unstable and stable limit cycle for $R_{vort}^{(1)}$ ($\zeta >0$, a) and $R_{flow-vort}^{(1)}$ ($\zeta <0$, b), respectively.

The numerically integrated trajectories for parameters chosen in $R_{vort}^{(1)}$ and $R_{flow-vort}^{(1)}$ for $c=0.005$ are shown in figures 6 and 7, respectively. The prediction of the spiral exit of the orientation trajectories from the FGP in $R_{vort}^{(1)}$ and the spiral approach towards it in $R_{flow-vort}^{(1)}$ mentioned above is confirmed from these figures. The spiralling rate in $R_{vort}^{(1)}$ increases with $c\, De$ and $\kappa$ (not shown).

Figure 6. In $R_{vort}^{(1)}$ (shown here for $c=0.005$, $\kappa =50, c\, De=0.01$), trajectories starting near the FGP (exemplified here with the blue trajectory) spiral out of the plane. Globally they approach the same stable limit cycle as the trajectories starting near the vorticity axis (exemplified here with the orange trajectory). The blue trajectory starting near the flow direction spans a larger portion of phase space in $p_2$, but we show the region close to the flow–vorticity plane to highlight the limit cycle.

Figure 7. In $R_{flow-vort}^{(1)}$ (shown here for $c=0.005$, $\kappa =10, c\, De=0.48$), trajectories of particle orientation starting near the FGP (blue) spiral into the plane. Globally they emanate from an unstable limit cycle – the boundary between blue and green trajectories (that are started very close to each other in this numerical integration). There is a stable limit cycle above this unstable limit cycle at the boundary between green and orange trajectories.

As $b^2$ is reduced by increasing $c\, De$ or $\kappa$, the time period, $2{\rm \pi} /b$, increases. The effect of viscoelasticity (increasing $c\, De$) to increase the time period is consistent with the experimental observations of Gauthier et al. (1971) and Bartram et al. (1975). As mentioned earlier in § 5.1, when $b^2=0$, due to an infinite period bifurcation (Strogatz 2018) two fixed points emerge on the FGP, and we discuss the trajectories off the FGP in the $b^2<0$ regime next.

5.3.1. Monotonic behaviour near FGP: ${b^2<0}$

Equation (5.2) has two fixed points at $p_{2,{flow}}^{0-}$ and $p_{2,{flow}}^{0--}$ from (5.5a,b) (or their more accurate values in (5.6)). For the dynamical system in the $p_2\unicode{x2013} p_3$ plane near the flow direction, defined by the system of (5.2) and (5.10), these fixed points are $[\,p_{2,{flow}}^{0-},0]$ and $[\,p_{2,{flow}}^{0--},0]$. The eigenvectors at each of the fixed points are the same, i.e.

(5.13a,b)\begin{equation} \boldsymbol{v}_{1,{flow}}=\left[1\quad 0\right],\quad \boldsymbol{v}_{2,{flow}}=\left[0\quad 1\right]. \end{equation}

The corresponding eigenvalues at $[\,p_{2,{flow}}^{0-},0]$ and $[\,p_{2,{flow}}^{0--},0]$ are $\lambda _{i,{flow}}^{0-}$ and $\lambda _{i,{flow}}^{0--}$, where $i\in [1,2]$,

(5.14)\begin{equation} \left.\begin{gathered} \lambda_{1,{flow}}^{0-}=-2\sqrt{-b^2},\quad \lambda_{2,{flow}}^{0-}=\zeta c\, De -(\zeta+4/\kappa^2) (4\log(2\kappa)-6)\sqrt{-b^2}, \\ \lambda_{1,{flow}}^{0--}=2\sqrt{-b^2},\quad \lambda_{2,{flow}}^{0--}=\zeta c\, De +(\zeta+4/\kappa^2) (4\log(2\kappa)-6)\sqrt{-b^2}. \end{gathered}\right\}\end{equation}

Both the fixed points are hyperbolic in the $b^2<0$ regime. Thus, similar to the fixed point at the vorticity axis using the Hartman–Grobman theorem, the linear stability of the fixed points allows qualitative insight into the complete nonlinear behaviour of particle orientation. Both eigenvalues of both the fixed points are real in the $b^2<0$ regime. The first eigenvalues ($\lambda _{1,{flow}}^{0-}$ and $\lambda _{1,{flow}}^{0--}$) of both the fixed points do not change sign in the $b^2<0$ regime, i.e. $\lambda _{1,{flow}}^{0-}<0$ and $\lambda _{1,{flow}}^{0--}>0$. But, the second eigenvalues ($\lambda _{1,{flow}}^{0-}$ and $\lambda _{1,{flow}}^{0--}$), that are initially positive in the $b^2<0$ regime, become negative at different parameter ($c$, $De$ and $\kappa$) values. Here $\lambda _{2,{flow}}^{0-}$ becomes negative first. $\lambda _{2,{flow}}^{0-}=0$ is the bifurcation boundary where $[\,p_{2,{flow}}^{0-},0]$ changes from a saddle point (stable manifold along gradient direction and unstable along vorticity) to a stable fixed point. Here $\lambda _{2,{flow}}^{0--}=0$ is the bifurcation boundary where $[\,p_{2,{flow}}^{0--},0]$ changes from an unstable node to a saddle point (stable manifold along vorticity direction and unstable along gradient). Therefore, three new regimes arise within the $b^2<0$ region that are labelled in figure 2 as (1) $R_{vort}^{(2)}$ ($\lambda _{1,{flow}}^{0-}>0$ and $\lambda _{1,{flow}}^{0--}>0$), (2) $R_{flow-vort}^{(2)}$ ($\lambda _{1,{flow}}^{0-}<0$ and $\lambda _{1,{flow}}^{0--}>0$) and (3) $R_{flow-vort}^{(3)}+R_{flow}$ ($\lambda _{1,{flow}}^{0-}<0$ and $\lambda _{1,{flow}}^{0--}<0$). The phase flow close to the fixed points in the gradient–vorticity plane for these three cases is shown in figure 8. We can observe the phase flow approaching $[\,p_{2,{flow}}^{0-},0]$ and leaving $[\,p_{2,{flow}}^{0--},0]$ along the gradient direction ($p_2$) for all three cases. This implies that a particle with an initial orientation close to the FGP will approach the flow direction (specifically the fixed point $[\,p_{2,{flow}}^{0-},0]$) while slowing down. For the parameters within $R_{vort}^{(2)}$, the particle will then monotonically drift away from the FGP. This is similar to the experimental observation of Bartram et al. (1975) discussed in § 1. For a given $\kappa$ and $c$, at larger $c\, De$ than $R_{vort}^{(2)}$, within the regions $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the particle starting near the FGP (and for all starting orientations in $R_{flow}$) achieves a stable orientation near the flow direction, similar to the experiments of Iso et al. (1996a).

Figure 8. Phase portraits of the system of (5.2) and (5.10) in the $b^2<0$ regime. In this regime, two fixed points exist on the FGP ($p_3=0$) close to the flow direction. Both the fixed points are downstream of the flow–vorticity plane ($p_2=0$). An unstable (red marker) and a saddle (green marker) node with its unstable manifold along the $p_3$ axis in (a) ($R_{vort}^{(2)}$) indicates a monotonic drift of the particle away from the flow gradient plane. A stable fixed point (blue marker) near the flow direction ($p_2\approx 0, p_3=0$) in (b) ($R_{flow-vort}^{(2)}$) and (c) ($R_{flow-vort}^{(3)}$) panels indicate that particles with starting orientation near the FGP may align near the flow direction. The presence of an unstable point in $R_{flow-vort}^{(2)}$ (b) instead of a saddle point with its stable manifold perpendicular to the flow gradient plane in $R_{flow-vort}^{(2)}$ (c) in addition to the stable point in these cases indicates a lower proportion of trajectories leading to flow alignment in $R_{flow-vort}^{(2)}$ than in $R_{flow-vort}^{(3)}$.

The locations of fixed points corresponding to $\tilde {p}_{2,{flow}}^{0-}$ and $\tilde {p}_{2,{flow}}^{0--}$ in the orientation space are $\left [ \pm \sqrt {1-(\tilde {p}_{2,{flow}}^{0-})^2}\ \pm \tilde {p}_{2,{flow}}^{0-} \ 0 \right ]$ and $\left [\pm \sqrt {1-(\tilde {p}_{2,{flow}}^{0--})^2} \ \pm \tilde {p}_{2,{flow}}^{0--} \ 0 \right ]$. These locations match the corresponding fixed points in the numerically integrated trajectories in $R_{vort}^{(2)}$, $R_{flow-vort}^{(2)}$, and, $R_{flow-vort}^{(3)}+R_{flow}$ shown in the plots of figures 9, 10 and 11, respectively. We mark the analytical locations of these fixed points with two different coloured markers on the FGP ($p_3=0$) in figures 9, 10 and 11 and show the nearby trajectories to approach/leave these locations in the fashion described by the linear stability theory discussed above.

Figure 9. In $R_{vort}^{(2)}$ (shown here for $c=0.005$, $\kappa =50, c\, De=0.3$), a particle's orientation trajectories approach the stable limit cycle around and near the vorticity axis. In contrast to $R_{vort}^{(1)}$ (figure 6) where the trajectories spiral away from the FGP, here they leave the FGP monotonically. We show the region close to the flow–vorticity plane as this is where the different attractors lie. The grey surface is the unit sphere, i.e. the orientation space.

Figure 10. In $R_{flow-vort}^{(2)}$ (shown here for $c=0.005$, $\kappa =100$ and $c\, De=0.9$) a particle's orientation trajectories that start close to the FGP (solid lines) approach the flow direction either on the same or the opposite side of the gradient–vorticity plane. Trajectories farther away from the FGP (dashed lines) approach the stable limit cycle near the vorticity axis. Due to a saddle and stable node close to each other on the FGP and a stable limit cycle near the vorticity axis, another saddle node emerges on the faster eigendirection of the stable node. The grey surface is the unit sphere, i.e. the orientation space.

Figure 11. In $R_{flow-vort}^{(3)}$ (shown here for $c=0.005$, $\kappa =100$ and $c\, De=1.2$) the behaviour of a particle's orientation trajectories is similar to that in $R_{flow-vort}^{(2)}$ shown in figure 10. The primary difference between $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ is that, in the former, the fixed point on the FGP farther from the flow direction is an unstable fixed point (figure 10), while it is a saddle node in the latter (shown here). Additionally, in $R_{flow-vort}^{(3)}$, there is an unstable node at the intersection of the stable manifold of the saddle points in FGP and flow–vorticity plane. Trajectories with solid lines end up at one of the fixed points near the flow direction, and those with dashed lines end in the stable limit cycle near the vorticity axis.

Numerical integration of the governing equation shown in figure 9 confirms the existence of the stable fixed point and unstable node predicted by linear stability analysis. In agreement with the $p_2-p_3$ phase plane analysis above, the particle starting near the flow gradient plane in $R_{vort}^{(2)}$ approaches the flow direction. It then monotonically departs the FGP along the flow–vorticity plane. The primary difference between $R_{vort}^{(1)}$ and $R_{vort}^{(2)}$ is that the particle leaves the FGP spirally in the former (figure 6) and monotonically in the latter (figure 9). As mentioned earlier, these trajectories are reminiscent of that in the experimental observations of Bartram et al. (1975). We can only make qualitative comparisons with Bartram et al. (1975) due to unknown non-Newtonian properties of their viscoelastic fluid. Furthermore, the velocity field around the $\kappa =9.1$ particles used in their experiments is likely to have a quantitative difference from the SBT approximations used in our theory.

In $R_{flow-vort}^{(2)}$, the numerically integrated trajectories shown in figure 10 reveal the presence of stable and unstable nodes in the FGP at the locations predicted by the linear stability analysis. Similarly, the existence of a stable node close to the flow direction and a saddle point farther away on the FGP is confirmed from the numerical integration shown in figure 11 for the parameters chosen in $R_{flow-vort}^{(3)}$. In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ the trajectories starting close to the FGP approach a stable orientation near the flow direction (stable fixed point at $\tilde {p}_{2,{flow}}^{0-}$).

5.4. Global particle orientation dynamics

In § 5.1 we analysed the equations for the rotational motion of a particle at the vorticity axis and within the FGP. We found that viscoelasticity does not lead to a non-zero $\dot {\boldsymbol {p}}$ at the vorticity axis and it does not change the log-rolling rotational velocity. In the FGP, we observed that a small amount of viscoelasticity (small $c\, De$) leads to a larger orbit time period. Further increasing $c\, De$ leads to an infinite period bifurcation and birth of two fixed points, thereby leading to particle migration towards the flow direction when in the FGP. In § 5.2, by linearizing the governing equation at the vorticity axis, we showed that the vorticity axis is an unstable spiral for $c\, De<\kappa$. For $c\, De>\kappa$, it is an unstable node. Therefore, a small perturbation of the particle orientation from the vorticity axis leads to departure from the axis. This was corroborated by the orientation trajectories obtained by numerical integration of the system of governing equations. By analysing the system and its fixed points (when they exist) near the flow direction in the FGP and performing numerical integration of the trajectories starting near but not on the FGP in § 5.3, we determined the orientational dynamics in the region nearby the flow gradient plane. The particle either spirals away from or towards the FGP when there are no fixed points in FGP. When the fixed points in FGP arise, the one closer to the flow direction is either a saddle point (with its stable manifold in the FGP) or a stable node. Therefore, the particle first migrates towards the flow direction nearly parallel to the FGP. Then it may either settle near the flow direction or monotonically drift away from the FGP along the flow–vorticity plane. The boundaries in the $c\, De-\kappa$ space separating these different qualitative behaviours were shown in figure 2 for $c=0.005$. The qualitative nature of these boundaries is not sensitive to the exact value of $c$. As we observed in figures 6, 7, 9, 10 and 11, there are other invariant features (stable limit cycle, unstable limit cycle, saddle point, and unstable node) that arise in the orientation space off the FGP and the vorticity axis. These features may be viewed as a result of the orientation space being restricted to a unit sphere.

In $R_{vort}^{(1)}$ (figure 2), which extends to arbitrarily small (but finite) viscoelasticity (small $c\, De$), the vorticity axis is a spiral source, and the FGP is an unstable limit cycle. Therefore, at least one stable attractor must exist between these two regions on the unit sphere. In the case of a Newtonian fluid, the particle follows an initial-condition-dependent Jeffery orbit, i.e. one of a concatenation of neutrally stable periodic orbits centred at the vorticity axis. Therefore, at small (but finite) $c\, De$ the simplest change in phase-plane dynamics leads to a stable limit cycle around and near one of these limit cycles, as shown in figure 6. As shown by linear stability analysis at the vorticity direction in § 5.2 the rate of deviation of trajectories away from the vorticity axis is driven by the $O(1/\kappa ^{2})$ flow generated by the force doublet and $O(1/\kappa ^{2})$ Stokeslet. Thus, a slender particle starting at any orientation between the vorticity axis and the FGP leads to a final orientation behaviour where the particle undergoes periodic motion very close to and around the vorticity axis.

In $R_{flow-vort}^{(1)}$ (figure 2), the FGP changes from being an unstable to a stable limit cycle. Therefore, an unstable invariant object or a repeller in the form of an unstable limit cycle exists between the stable limit cycle near the vorticity direction and the FGP, as shown in figure 7. Thus, in $R_{flow-vort}^{(1)}$, a particle with an initial orientation between the unstable limit cycle and FGP spirals to the FGP, where it undergoes a tumbling motion. A particle with an initial orientation between the unstable limit cycle above the FGP and the stable limit cycle near the vorticity axis spirals towards the stable limit cycle near the vorticity axis. A particle released at an arbitrarily small (but finite) angle from the vorticity axis spirals outwards. In both these cases, eventually, the particle undergoes perpetual periodic motion on the stable limit cycle close to the vorticity axis.

In $R_{vort}^{(2)}$ (figure 2), the FGP contains two fixed points as shown in figure 9. One is an unstable node, and the other is a saddle point with its unstable manifold perpendicular to the FGP. Thus, no other invariant feature (attractor or repeller) is needed between the FGP and the stable limit cycle near the vorticity axis, which is the only stable invariant object on the orientation space. Hence a particle starting at an arbitrarily small angle from the vorticity axis or the FGP ends up in a periodic motion very close to the vorticity axis.

Similar to $R_{vort}^{(1)}$, $R_{vort}^{(2)}$ and $R_{flow-vort}^{(1)}$ discussed above, there is a stable limit cycle around the vorticity axis in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$. Therefore, a particle with an initial orientation arbitrarily close to the vorticity direction (but not on the vorticity axis) eventually undergoes a periodic motion around the vorticity axis. This stable limit cycle exists for all regions except for $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ and it will be further analysed in § 5.4.1.

In $R_{flow-vort}^{(2)}$ (figure 2) the fixed point in the FGP closest to the flow direction is a stable node and the other fixed point in the FGP is an unstable node as shown in figure 10. Hence, a saddle node exists near the flow–vorticity plane to repel the particle orientation trajectories towards the stable limit cycle near the vorticity axis on one side and the stable fixed point near the flow direction on the FGP on the other. This saddle point receives trajectories emanating from the unstable point on the FGP. Thus, depending upon the initial orientation of the particle, it may eventually either obtain a stable orientation near the flow direction or undergo a periodic motion close to the vorticity axis.

The dynamics in $R_{flow-vort}^{(3)}$ (figure 2) are similar to those in $R_{flow-vort}^{(2)}$, but the unstable fixed point in the FGP in $R_{flow-vort}^{(2)}$ changes to a saddle node in $R_{flow-vort}^{(3)}$ with no qualitative change to the other invariant objects present in $R_{flow-vort}^{(2)}$. Therefore, to allow smooth and continuous phase flow, an unstable node occurs at the intersection of the saddle points in the FGP and near the flow–vorticity plane, as shown in figure 11. At a given $\kappa$, the value of $c\, De$ in $R_{flow-vort}^{(3)}$ is larger than that in $R_{flow-vort}^{(2)}$ (figure 2) and the separation between two fixed points increases with increasing $c\, De$ as shown by figure 3. Due to this increased separation and the unstable node above the flow gradient plane, more of the phase flow is directed towards the FGP in $R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$. Thus, the basin of attraction for the stable fixed point near the flow direction is larger in $R_{flow-vort}^{(3)}$ as compared with $R_{flow-vort}^{(2)}$. In other words, more of the initial particle orientations lead to a stable final orientation near the flow direction in the $R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$ as observed by comparing the proportion of solid lines in figures 10 and 11. When the stable limit cycle near the vorticity axis exists, the saddle point on the flow–vorticity plane is important in determining the proportions of initial orientations leading to the flow alignment and towards the periodic motion around the vorticity axis. Its location will be further analysed in § 5.4.2.

The $c\, De|_{cut\textrm -off}^{Limit}$ boundary in figure 2, beyond which the stable limit cycle near the vorticity direction does not exist, may be viewed as arising due to the interaction between the saddle node near the flow–vorticity plane and the stable limit cycle near the vorticity direction as these move in the orientation space upon increasing $c\, De$. We will discuss this in § 5.4.3.

5.4.1. Stable limit cycle around the vorticity axis

The previous study of Leal (1975) concerning a slender particle in simple shear flow of a second-order fluid predicts that viscoelasticity (although incorrect in attributing it to the second instead of the first normal stress difference) drives the particle orientation towards the vorticity axis. Leal (1975) only considered the Stokeslet fluid flow (of order $O(1/\log (\kappa ))$ and $O(1/\log (\kappa )^2)$) from the slender body theory. This velocity field is proportional to $p_2$ and hence has no effect when the particle is in the flow–vorticity plane. However, our study shows that a stable limit cycle exists around the vorticity axis due to the competition between polymer-driven torque arising from the doublet flow causing the particle to spiral away from the vorticity axis and that from the Stokeslet flow at $O(1/\log (\kappa ))$ causing it move towards the region around vorticity.

The real part of the eigenvalues of the fixed point at the vorticity axis in $R_{vort}^{(1)}$ is $2c\,De/\kappa ^2$ (5.8a,b), and the rate at which trajectories leave the FGP is $c\,De\zeta$ (5.11). Figure 12 shows the variation of $2c\,De/\kappa ^2$ and $c\,De\zeta$ with $\kappa$ for a few values of $c\, De$ at $c=0.005$. From this figure, we conclude that for $\kappa \gtrapprox 15$, trajectories leave the flow gradient plane much faster than the outward spiralling rate from the vorticity axis. Hence, the stable limit cycle shifts closer to the vorticity axis as $\kappa$ increases. Furthermore, as $\kappa$ increases beyond $\kappa \approx 20$, $c\,De\zeta$ remains almost constant while the rate of spiralling from the vorticity axis, $\kappa ^{-2}$, decreases. Therefore, the stable limit cycle approaches the vorticity axes at larger $\kappa$.

Figure 12. The variation with $\kappa$ of the rate of deviation of orientation trajectories away from the vorticity axis, $2c\,De/\kappa ^2$ (orange), and the flow gradient plane, $c\,De\zeta$ (black), in $R_{vort}^{(1)}$ ($b^2>0$) for a few values of $c\, De$ at $c=0.005$. On each curve, $b^2>0$ is to the left of the solid markers. An increasing gap between orange and black curves, with almost horizontal black curves, with $\kappa$ indicates a larger increase in the departure rate of the orientation trajectories from the flow gradient plane than that from the vorticity axis, implying a shift of the stable limit cycle closer to the vorticity axis.

We quantify the location of the limit cycle observed in the numerical integration of (4.46). For this purpose we transform this equation into the $C-\tau$ coordinate system (Leal & Hinch 1971), using

(5.15a,b)\begin{equation} C=\frac{\sqrt{p_2^2+p_1^2/\kappa^2}}{p_3}, \quad\tau=\frac{p_1}{p_2}, \end{equation}

to obtain

(5.16)$$\begin{gather} \frac{{\rm d} C}{{\rm d}t}=c\,De\tau^2\left(\frac{4(1+C^2)}{C\kappa^4}-\frac{C}{(4\log(2\kappa)-3)(\kappa^2+\tau^2)}\right)p_2^2+c\,DeO(\,p_2^4), \end{gather}$$

(5.17)$$\begin{gather}\frac{{\rm d}\tau}{{\rm d}t}=1+c\,De\tau p_2^2\left(\frac{\tau^2}{4\log(2\kappa)-3}+4\frac{\kappa^2+\tau^2}{C^2\kappa^4}\right)+O(c\,De). \end{gather}$$

In the Newtonian limit, ${\textrm {d}C}/{\textrm {d}t}=0$ and the particle trajectory is determined by its initial condition $C(t)=C(0)=C_N$. Each $C_N=[0,\infty ]$ represents a different Jeffery orbit. Here $C=0$ represents the log-rolling motion, i.e. the particle rotating about its major axis which is aligned with the vorticity axis, and $C=\infty$ represents the major axis rotating within the FGP. Here $\tau =t$ represents the phase on the Jeffery orbit. A periodic orbit can thus be ascertained using a single parameter $C$ that repeats periodically. In the range, $10\le \kappa \le 200$, for $c\, De\ge 0.01$ we evolve the system of equations for an initial condition starting near the vorticity direction, at $C=10^{-7}$. We evolve the equations for each choice of $c\, De$ and $\kappa$ until the change in the period-averaged $C$ across 100 successive Jeffery periods is less than $10^{-5}$. Finally, we define $\bar {C}_{Limit}^{numerical}$ as the average $C$ over the last 100 Jeffery periods. The $\bar {C}_{Limit}^{numerical}$ versus $c\, De$ for a few values of $\kappa \in [10,200]$ is shown in figure 13(a).

Figure 13. (a) The location of the stable limit cycle at different $c\, De$ and $\kappa$ at $c=0.005$ (nearly identical curves are obtained for other polymer concentrations in the range $0.01\le c\le 0.2$) indicated by the average $C=\sqrt {p_2^2+p_1^2/\kappa ^2}/{p_3}$, $\bar {C}_{Limit}^{numerical}$ on the limit cycle. Here $C=0$ is the vorticity axis and $C=\infty$ is the FGP. (b) The $\bar {C}_{Limit}^{numerical}$ on the limit cycle versus $\kappa$ at low $c\, De$.

Now, $\bar {C}_{Limit}^{numerical}$ is smaller than $0.1$ for $\kappa >10$, decreases with $\kappa$ and varies slightly with $c\, De$ at each $\kappa$. For $\kappa >50$ $\bar {C}_{Limit}^{numerical}$ is less than $0.01$ for a range of $c\, De$. Therefore, the limit cycle is close to the vorticity axis. Beyond a certain value of $c\, De=c\, De|_{cut\textrm -off}^{Limit}$, the limit cycle does not exist. This marks the boundary between $R_{flow-vort}^{(2)}$ and $R_{flow}^{(1)}$ for a given $\kappa$, as shown in figure 2. Now, $c\, De|_{cut\textrm -off}^{Limit}$ increases with $\kappa$ and for $10\le \kappa \le 200$ we find

(5.18)\begin{equation} c\, De|_{cut\text-off}^{Limit}\approx 0.1056\kappa-0.5183. \end{equation}

At larger $\kappa$, the stable limit cycle is closer to the vorticity axes and exists for a larger range of $c\, De$. Hence, marked as $R_{flow-vort}^{(1)}$ and $R_{flow-vort}^{(2)}$ in figure 2, there is a significant range of $c\, De$ at large $\kappa$ where two possible types of orientation dynamics are possible: a steady state orientation close to the flow direction and a periodic orbit near the vorticity direction. For $\kappa \lessapprox 20$, increasing $c\, De$ from very small values does not affect the position of the limit cycle much. However, as $c\, De$ approaches $c\, De|_{cut\textrm -off}^{Limit}$, the limit cycle moves away from the vorticity axis. For $\kappa \gtrapprox 20$, the limit cycle moves towards the vorticity axis upon increasing $c\, De$ from very small values. However, beyond a given $c\, De$, up to $c\, De|_{cut\textrm -off}^{Limit}$ the limit cycle's angular position is not influenced by $c\, De$. These features are reflected by the trends in $\bar {C}_{Limit}^{numerical}$ in figure 13(a). At low $c\, De$, the limit cycle is very close to one of the degenerate Jeffery orbits of the Newtonian case, i.e. $C$ does not change much on the limit cycle. The variation of $\bar {C}_{Limit}^{numerical}$ with $\kappa$ at low $c\, De=0.01, 0.02$ and 0.03 is shown in figure 13(b), where we can observe $\bar {C}_{Limit}^{numerical}$ to scale as $\kappa ^{-1.5}\log (2\kappa -3)$ for $10\le \kappa \le 100$ and approximately $\kappa ^{-1}$ for larger $\kappa$.

We find the particle orientation trajectories approaching this stable limit cycle close to the vorticity axis instead of spiralling towards the vorticity axis as concluded from the experiments of Gauthier et al. (1971) and Iso et al. (1996b). Gauthier et al. (1971) claim the approach towards the vorticity axis by extrapolating the observed data. In some of the experimental observations of Iso et al. (1996b) at small elasticity after undergoing initial spiralling, the particle oscillates around the vorticity axis without settling at a steady orientation. This may indicate the limit cycle discussed above. The theoretical rate of spiralling away from the vorticity axis ($2c\,De/\kappa ^2$) is small (figure 12), and small imperfections from the Couette cell might have led to secondary flows in the experiments. Thus, the existence and size of the stable limit cycle around the vorticity axis must be further tested in future experiments and numerical simulations.

5.4.2. Saddle node near the flow–vorticity plane in $R_{flow-vort}^{(2)}$, $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$

From the trajectories of particle orientation in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ shown in figures 10 and 11 we find the stable limit cycle around the vorticity axis and the stable node on the flow gradient plane near the flow axis to be separated by a saddle point. The saddle point arises on the boundary between $R_{vort}^{(2)}$ and $R_{flow-vort}^{(2)}$ when the saddle point near the flow direction on the FGP changes into a stable node. It also arises on the boundary between $R_{flow-vort}^{(3)}$ and $R_{flow-vort}^{(1)}$ when the unstable limit cycle in the FGP vanishes. The saddle point persists in the region $R_{flow-vort}^{(3)}$ when approaching it by increasing $c\, De$ from $R_{flow-vort}^{(2)}$. In the regions $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the unstable direction of the saddle point leads to the limit cycle on one side, and the stable node on the other. Based on numerical evidence from figures 10 and 11 we assume the saddle point, $\boldsymbol {p}_{saddle}=[{p}_{1,saddle}\ {p}_{2,saddle}\ {p}_{3,saddle}]$, to approximately lie at the same $p_2$ coordinate as the stable node in the flow gradient plane:

(5.19a,b)\begin{equation} {p}_{2,saddle}={p}_{2,{flow}}^{0-}+{p}_{2,saddle}',\quad {p}_{2,saddle}' \ll {p}_{2,{flow}}^{0-} \end{equation}

where ${p}_{2,{flow}}^{0-}$ is given by (5.5a,b). From $\dot {p}_{2,saddle}=0$, we obtain

(5.20)\begin{equation} {p}_{2,saddle}=\frac{\sqrt{-b^2}c\, De+(4\log(2\kappa)-6)b^2}{c\, De(1-p_{1,saddle})-\sqrt{-b^2}(4\log(2\kappa)-6)}.\end{equation}

From (4.46), $\dot {p}_{1,saddle}=\dot {p}_{3,saddle}=0$ is equivalent to

(5.21)\begin{equation} p_{1,saddle} \alpha|_{p_1=p_{1,saddle},p_3=p_{3,saddle}}= p_{2,saddle}/(c\, De). \end{equation}

A general solution to the above equation is intractable, but we can solve it in the limit of small $p_{1,saddle}$ or small $p_{3,saddle}$. We denote the value of $p_{1,saddle}$ in the limit of small $p_{1,saddle}$ as $\tilde {p}_{1,saddle}$ and the value of $p_{3,saddle}$ in the limit of small $p_{3,saddle}$ as $\tilde {p}_{3,saddle}$. We obtain

(5.22) \begin{equation} \tilde{p}_{1,saddle}=\frac{\begin{array}{l}-8 (c\, De)^2 f^2+16 c \,De f^3 \sqrt{-b^2}\\ -\sqrt{\begin{array}{l} 4c \,De f \left(-(c \,De)^2 \kappa^2+2 (c \,De) f \kappa^2 \sqrt{-b^2}+4 f^2\right) \\ 4\left(-3 {\rm \pi}(c \,De)^2 \kappa^2+2 c \,De f \left(8 f+3 {\rm \pi}\kappa^2 \sqrt{-b^2}\right)+12 {\rm \pi}f^2\right)\\ +\,128 c \,De f^3 \left((c \,De)-2 f \sqrt{-b^2}\right)^2 \end{array}}\end{array}} {(c \,De) \left(3 {\rm \pi} (c \,De)^2 \kappa^2-2 (c \,De) f \left(8 f+3 {\rm \pi} \kappa^2 \sqrt{-b^2}\right)-12 {\rm \pi} f^2\right)},\end{equation}

and

(5.23) \begin{align} \tilde{p}_{3,saddle}=-\frac{\begin{array}{l} 3 {\rm \pi} (c \,De) \left((c \,De)^2 \kappa^2-2 c \,De f \kappa^2 \sqrt{-b^2}-4 f^2\right)\\ +\sqrt{\begin{array}{l} 9 {\rm \pi}^2 (c \,De)^2 \left(-(c \,De)^2 \kappa^2+2 c \,De f \kappa^2 \sqrt{-b^2}+4 f^2\right)^2\\ +\dfrac{1}{2}( c \,De)^3 \kappa^2 \left(3 {\rm \pi}^2 c \,De-2 f \left(3 {\rm \pi}^2 \sqrt{-b^2}+2\right)\right)\\ +\dfrac{1}{2} c \,Def^2\left(\,f \left(16-64 c \,De \sqrt{-b^2}\right)+c \,De \left(8 \kappa^2 \sqrt{-b^2}-12 {\rm \pi}^2\right)\right) \\ \left(3 {\rm \pi}^2 (c \,De)^2 \kappa^2-2 c \,De f \left(16 f+3 {\rm \pi}^2 \kappa^2 \sqrt{-b^2}\right)-4 f^2 \left(16 f \sqrt{-b^2}+3 {\rm \pi} ^2\right)\right)\end{array}}\end{array}}{2 \left(\dfrac{3}{4} {\rm \pi}^2 c \,De \left((c \,De)^2 \kappa^2-2 c \,De f \kappa^2 \sqrt{-b^2}-4 f^2\right)-8 c \,De f^2 \left((c \,De)+2 f \sqrt{-b^2}\right)\right)},\end{align}

where $f=2\log (2\kappa )-3$ and $\sqrt {-b^2}$ is in (5.12a,b). The variation of $\tilde {p}_{3,saddle}$ and $\sqrt {1-\tilde {p}_{1,saddle}^2-{p}_{2,saddle}^2}$ with $c\, De$ for a few $\kappa$ is shown in figure 14(a). From either of these approximations, we find that the saddle point moves closer to the vorticity axis as $c\, De$ is increased or $\kappa$ is reduced. Both approximations yield a similar location for the saddle point. The best approximation for the location of the saddle point for any $c\, De$ and $\kappa$, is $\boldsymbol {p}_{saddle}=[{p}_{1,saddle}\ {p}_{2,saddle}\ {p}_{3,saddle}]$, where

(5.24)\begin{gather} \left.\begin{gathered} {p}_{1,saddle}=\tilde{p}_{1,saddle},\quad {p}_{3,saddle}=\sqrt{1-\tilde{p}_{1,saddle}^2-{p}_{2,saddle}^2}, \quad\text{if }|\tilde{p}_{1,saddle}|\le|\tilde{p}_{3,saddle}|\\ {p}_{1,saddle}=\sqrt{1-{p}_{2,saddle}^2-\tilde{p}_{3,saddle}^2},\quad {p}_{3,saddle}=\tilde{p}_{3,saddle},\quad \text{if } |\tilde{p}_{1,saddle}|>|\tilde{p}_{3,saddle}| \end{gathered},\right\}\end{gather}

and $p_2$ is from (5.20). We show this saddle point's more accurate location of $p_3$ in figure 14(b). The location of this analytically predicted saddle point, shown with a green marker in figures 10, 11 and 15, agrees well with that inferred from the of numerically integrated trajectories.

Figure 14. The $p_3$ coordinate of the saddle point in the flow–vorticity plane in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$. (a) Solid lines (the graph of $\sqrt {1-\tilde {p}_{1,saddle}^2-p_{2,saddle}^2}$ from (5.22) and (5.20)) provide the more accurate $p_3$ location of the saddle point when it is closer to the vorticity axis ($p_{3,saddle}\approx 1$) and dashed lines (the graph of (5.23)) represent the more accurate $p_3$ location of the saddle point when it is closer to the flow axis ($p_{3,saddle}\approx 0$). (b) Most accurate location of the saddle point constructed from (a) using the crossover point between the solid and dashed lines for each $\kappa$, i.e. using (5.24).

Figure 15. Trajectories at large $c\, De$ for $\kappa =20$ in $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ zoomed near the flow–vorticity plane. At very large $c\, De$ (corresponding to regions $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$), the stable limit cycle does not exist around the vorticity axis and a particle starting anywhere apart from vorticity axis ends up being nearly flow aligned.

In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the unstable manifold of this saddle point leads to two stable attractors, i.e. either a limit cycle near the vorticity axis or a stable node near the flow direction. Therefore, the location of the saddle point partially dictates the relative sizes of the basins of attraction of these stable attractors. As $c\, De$ increases, due to the changing location of the saddle point, $\boldsymbol {p}_{saddle}$, the particle is more likely to attain a final orientation towards the flow axis as compared with approaching the limit cycle around the vorticity axis (see figure 14 and compare the proportion of the solid lines in 10 and 11). The saddle point also persists in the large $c\, De$ regions $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ as shown in the different plots in figure 15.

5.4.3. The $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$: very large $c\, De$ guarantees flow alignment

In § 5.4.1, we found that for $\kappa \gtrapprox 20$, the stable limit cycle moves closer to the vorticity axis upon increasing $c\, De$. The saddle point near the flow–vorticity plane that exists in regions shown in figure 2 with $c\, De$ equal and larger than that in $R_{flow-vort}^{(2)}$ also moves towards the vorticity axis upon increasing $c\, De$ as discussed in § 5.4.2. The saddle point moves towards the vorticity axis faster than the limit cycle. Beyond a $\kappa$-dependent value of $c\, De$, when the saddle point is too close to the limit cycle, the latter ceases to exist, making the stable fixed point near the flow direction the only stable attractor. This boundary is marked as $c\, De|_{cut\textrm -off}^{Limit}$ in figure 2. This implies that a particle with any initial orientation other than the vorticity axis ends up being aligned near the flow direction. Plots corresponding to $c\, De=1.2$ and 1.8 for $\kappa =20$ in figure 15 show the vanishing of the stable limit cycle near the vorticity axis upon going from $R_{flow-vort}^{(3)}$ to $R_{flow}^{(1)}$. For $c\, De=1.2$ ($R_{flow-vort}^{(3)}$) the saddle point is just outside the limit cycle, whereas for $c\, De=1.8$ ($R_{flow}^{(1)}$), corresponding to a location close to the $R_{flow-vort}^{(3)}-R_{flow}^{(1)}$ boundary, the limit cycle does not exist and instead the stable manifolds of the saddle point originate from the unstable spiral at the vorticity axis. The unstable manifold of the saddle point leads to the stable node near the flow direction. The saddle point persists for larger $c\, De$ in $R_{flow}^{(1)}$ and also $R_{flow}^{(2)}$ as show by the $c\, De=4$ and 12 plots at $\kappa =20$, respectively, in figure 15. As also discussed in § 5.2 (figure 4a) in $R_{flow}^{(2)}$, a particle leaving the vorticity axis changes from spiralling to monotonic trajectories as the vorticity axis changes from an unstable spiral to an unstable node when going from $R_{flow}^{(1)}$ to $R_{flow}^{(2)}$. From figure 15(c) ($c\, De=4$ at $\kappa =20$) and figure 15(d) ($c\, De=12$ at $\kappa =20$), the proximity of the saddle point to the vorticity axis can also explain the breakdown of spiral motion near the vorticity axis in $R_{flow}^{(2)}$.

6. Conclusion

Using a regular perturbation expansion in polymer concentration, $c$, the balance of torques on the particle surface at each order in $c$, and the velocity field generated by a slender prolate spheroid (obtained from the slender body theories of Batchelor (1970) and Cox (1971)), we develop a theory to characterize the orientation dynamics of a freely rotating (torque-free) particle in simple shear flow of a viscoelastic fluid with small polymer concentration, $c$, and an arbitrary polymer relaxation time in the absence of inertia. This theory predicts a wide variety of orientation behaviours that are qualitatively similar to previous experimental observations of Gauthier et al. (1971), Bartram et al. (1975), and Iso et al. (1996a,b).

In the absence of inertia, for a viscoelastic fluid, where the fluid stress is a sum of the solvent and the polymer or non-Newtonian stress, the momentum equation can be decomposed into a Newtonian and a non-Newtonian part. The Newtonian part comprises the particle's motion in a Newtonian fluid undergoing the imposed flow and leads to the Newtonian stress field. The non-Newtonian part has zero velocity at the boundaries and comprises the balance of the divergence of the polymer-induced solvent stress and the divergence of the polymer stress. Thus, three physically distinct stress mechanisms impacting the particle surface are the elastic or polymer stress, the polymer-induced solvent stress, and the Newtonian stress. The decomposition of torques in this way can be useful in obtaining further insights into the particle motion in a viscoelastic or a non-Newtonian fluid where the fluid stress can be decomposed into a Newtonian and a non-Newtonian part. In addition to the non-Newtonian momentum equation, different stress components are coupled with the torque-free boundary condition and the polymer constitutive equation. The balance of torques generated by the three stresses leads to the appropriate angular velocity to allow torque-free particle motion. The boundary conditions for the motion of the particle and the imposed flow are contained in the Newtonian part of the momentum equation. Hence, the torque generated from the corresponding stress is termed the MIST. The polymer constitutive equation is driven by the sum of the velocity field from the Newtonian and non-Newtonian momentum components. The elastic torque is a function of the polymer stress on the particle surface. By definition, the PIST is the antisymmetric first moment of the polymer-induced solvent stress over the particle surface. However, using a generalized reciprocal theorem and the non-Newtonian momentum equation, we can express the PIST directly as a volume integral of the polymer stress.

Using a regular perturbation in the polymer concentration, $c$, we obtain the equations for particle's orientation dynamics in a small $c$ viscoelastic fluid subjected to a simple shear flow. At $O(1)$, the elastic and the polymer-induced solvent stress are zero. Thus, the particle orientation dynamics are the same as that in the simple shear flow of a Newtonian fluid, i.e. the particle undergoes Jeffery rotations. At $O(c)$, all three mechanisms mentioned above lead to a finite torque. The $O(c)$ Newtonian stress is that on a particle rotating (at $O(c)$ velocity) in a quiescent Newtonian fluid. Therefore, the $O(c)$ rotation rate is the one that generates enough $O(c)$ MIST to balance the PIST and the elastic torque. The total $O(c)$ rotation rate may be decomposed as the sum of $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$ and $\dot {\boldsymbol {p}}^{(1)}_{PIST}$, i.e. the $O(c)$ elastic and PIST generated rotation rates. The leading-order velocity field is taken from the slender body theory of Batchelor (1970) and Cox (1971), and this drives the $O(c)$ polymer stress. As both the elastic torque and the PIST can be expressed as a function of the polymer stress, the $O(c)$ rotation rate of the particle is obtained from the analytical velocity field from the slender body theory.

The polymer relaxation time ($\lambda$) is non-dimensionalized with the shear rate ($\dot {\gamma }$) to yield Deborah number, $De(=\lambda \dot {\gamma })$ ($c = 0$ or $De = 0$ implies a Newtonian fluid). We find $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$, to be independent of $De$ and initially consider $\dot {\boldsymbol {p}}^{(1)}_{PIST}$, separately in the low and high $De$ limits. In the low $De$ limit, $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$ and $O(De^0)$ $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ are the rotation rates due to the additional torque from the rate of strain and fluid pressure, respectively, of the Newtonian fluid with an additional viscosity $c$. Therefore, in the low $De$ limit, the sum of $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$ and the $O(De^0)$ contribution to $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ is identical to the rotation due to the Newtonian torque. Since the Newtonian torque is already balanced to be zero at the leading order in $c$, the first viscoelastic effect on the rotation rate for $De \ll 1$ arises at $O(c \,De)$ and is entirely from the polymer-induced solvent stress. In the high $De$ limit, $\dot {\boldsymbol {p}}^{(1)}_{PIST}\sim O(De)$ and is hence $O(De)$ larger than $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$. Therefore, the particle-induced polymer stress is the primary agent changing the particle dynamics in both the low and high $De$ limits. By interpolating the rotation rate in $De$ between the large and small (up to $O(De)$) $De$ limits, we obtain a uniformly valid equation for the particle orientation. We analyse the effect of particle aspect ratio, $\kappa$ and viscoelastic fluid parameters, polymer concentration ($c$), and relaxation rate ($De$) on this equation valid for any $De$.

Depending upon $c\, De$ and particle aspect ratio ($\kappa$), several qualitatively different particle orientation behaviours are possible. At constant $c\, De$, $c$ affects the behaviour only quantitatively. In a Newtonian fluid ($c\, De=0$), a particle undergoes periodic motion on an initial-condition-dependent Jeffery orbit around the vorticity axis. As expected from symmetry, the vorticity axis and the FGP remain invariant within the orientational space (unit sphere). Viscoelasticity (for all $c\, De$ and $\kappa$) does not affect the rotation rate of the log-rolling state on the vorticity axis. In the FGP, even very small $c\, De$ creates a bottleneck near the flow direction where the particle slows down in proportion to $c\, De$. Thus, viscoelasticity increases the period of tumbling motion within the FGP. At $c\, De=c\,De_{crit}=(4\log (2\kappa )-6)/\kappa$, an infinite period bifurcation occurs, and two fixed points arise in the FGP near the flow direction. Therefore, for $c\, De>c\,De_{crit}$, the tumbling motion within the FGP does not exist, and a particle initially oriented within FGP ultimately aligns near the flow direction.

In a Newtonian fluid, the phase flow of the dynamical system of the particle's orientational drift is of equal and opposite magnitude in the vorticity direction about the flow–vorticity plane. Thus the Jeffery orbits, in that case, are symmetric about the flow–vorticity plane. For $0< c\, De< c\,De_{crit}$ and $\kappa >\kappa _{FV1}$, the phase flow remains qualitatively similar, but the magnitude of the phase velocity in the vorticity direction is larger on the downstream side of the phase-flow. Thus, a particle starting close to (but not in) the FGP spirals towards the vorticity direction. $\kappa _{FV1}\approxeq 17$ for $c=0.005$ and slightly reduces with $c$. For $\kappa <\kappa _{FV1}$, spiralling away from the FGP occurs for a $c\, De$ up to a value slightly less than $c\,De_{crit}$ indicated by the $\zeta =0$ curve in figure 2. $R_{vort}^{(1)}$ in figure 2 is the region where the particle spirals away from the FGP. It continues to spiral towards the vorticity axis as it moves away from the FGP until it comes close to the axis. A particle with an initial orientation arbitrarily close to, but not on, the vorticity axis spirals away from the vorticity axis. Thus, in $R_{vort}^{(1)}$, a particle ultimately undergoes a periodic motion in a limit cycle near the vorticity axis. The spiralling towards the vorticity axis due to viscoelasticity is observed in the experiments of Gauthier et al. (1971) and Iso et al. (1996a,b).

For $\kappa <\kappa _{FV1}$ and between $c\, De$ corresponding to $\zeta =0$ and $c\, De=c\,De_{crit}$, marked as $R_{flow-vort}^{(1)}$ in figure 2, the phase velocity in the vorticity direction downstream of the flow–vorticity plane and close to the FGP changes direction as compared with the $R_{vort}^{(1)}$ region discussed above or the case for a Newtonian fluid. Thus, a particle with an initial orientation close to the FGP spirals into the FGP, where it ultimately undergoes tumbling motion (albeit with a larger period than in Newtonian fluid). A particle starting farther away from the FGP or near the vorticity axis spirals towards the periodic orbit (similar to $R_{vort}^{(1)}$) near the vorticity axis.

The phase space for the dynamical system of a particle in a Newtonian fluid consists of a concatenation of neutral centres. The periodic orbit at small angles away from the vorticity axis for the case of a viscoelastic fluid ($c\, De>0$) is a stable limit cycle. It is the only stable attractor in the phase space in $R_{vort}^{(1)}$. As $c\, De$ increases at large $\kappa \gtrapprox 50$ in the small $c\, De$ regime, the limit cycle shrinks as it goes towards the vorticity axis. The size of the limit cycle is not affected by $c\, De$ in the small $c\, De$ regime for $\kappa \lessapprox 20$. The previous studies by Leal (1975) using the slender body theory at low $De$ have considered only the $O(1/(2\log (2\kappa )-3))$ velocity disturbance and found viscoelasticity to cause a slender particle to spiral towards the vorticity axis. However, the polymeric torque due to the force doublet and $O(1/\kappa ^2)$ Stokeslet, $O(1/\kappa ^2)$ velocity disturbance from slender body theory of Cox (1971) that accounts for the fluid velocity disturbance the particle makes when in the flow–vorticity plane ($p_2=0$), makes the vorticity axis an unstable spiral. Thus, a stable limit cycle occurs from the competition of polymeric torques from the disturbance at $O(1/\kappa ^2)$ and the Stokeslet flow at order $O(\,p_2/(2\log (2\kappa )-3))$. Unlike Leal (1975) who attributes it to second normal stress difference, our theory predicts the effect of viscoelasticity in a second-order fluid to arise from the first normal stress difference of the fluid. In $R_{vort}^{(1)}$, the FGP is an unstable limit cycle and, in $R_{flow-vort}^{(1)}$, it is a stable limit cycle. In both cases, the stable limit cycle exists near the vorticity axis. Thus, in $R_{flow-vort}^{(1)}$, there is an unstable limit cycle between the two stable ones, and there are two possible final orientation behaviours, tumbling within the FGP and periodic motion close to the vorticity axis.

For $c\, De>c\,De_{crit}=(4\log (2\kappa )-6)/\kappa$, the fixed points in the FGP leads to monotonic particle drift towards or away from the FGP instead of spiralling motion. For $\kappa >\kappa _{FV1}$ when the fixed points first appear at $c\,De_{crit}$ the one closer to the flow direction is a saddle node, and the other is an unstable node. Thus, a particle starting close to the FGP migrates along the stable manifold of the saddle point up towards the flow direction and then monotonically drifts along the flow–vorticity plane towards the periodic orbit close to the vorticity axis. The region where this qualitative behaviour occurs is marked as $R_{vort}^{(2)}$ in figure 2 and the orientation trajectories are similar to the experimental observations of Bartram et al. (1975). Also, for $\kappa >\kappa _{FV1}$, upon increasing $c\, De$ beyond a certain value indicated by the $\lambda _{2,{flow}}^{0-}=0$ boundary in figure 2, the saddle point changes into a stable node and hence a particle with an initial orientation close to FGP ultimately becomes stably flow aligned. This is marked as $R_{flow-vort}^{(2)}$ in figure 2. In this region, depending upon the initial orientation, a particle finally obtains either a stable orientation close to the flow direction or undergoes periodic motion close to the vorticity axis. For $c\, De$ beyond the $\lambda _{2,{flow}}^{0--}=0$ boundary of figure 2 for $\kappa >\kappa _{FV1}$, or for $c\, De>c\,De_{crit}$ when $\kappa <\kappa _{FV1}$, the fixed point in the FGP away from the flow direction is a saddle node with its stable manifold perpendicular to the FGP. In the part of this region marked as $R_{flow-vort}^{(3)}$ in figure 2, the final particle behaviour is similar to that in $R_{flow-vort}^{(2)}$ discussed above. But, due to a saddle node instead of an unstable fixed point farther from the flow direction in the FGP, a greater proportion of the initial orientations in $R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$ lead the particle towards a flow aligned state. Upon increasing $c\, De$ beyond $c\, De|_{cut-off}^{Limit}$, in the regions marked as $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$, the particle at all initial orientations apart from the log-rolling state at the vorticity axis leads to a flow-aligned state. This is because the limit cycle close to the vorticity axis does not exist and the vorticity axis is an unstable spiral ($R_{flow}^{(1)}$) or node ($R_{flow}^{(2)}$). The flow-alignment of the particle was observed in the experiments of Iso et al. (1996a) at larger elasticity ($c$).

We are able to obtain qualitative agreement with the low shear rate experiments of Gauthier et al. (1971) and Iso et al. (1996b) which correspond to low $c\, De$ regime ($R_{vort}^{(1)}$ in figure 2). In these experiments, the particle spirals towards the vorticity axis, and our theory predicts spiralling towards the stable limit cycle near the vorticity axis. While the authors Gauthier et al. (1971) and Iso et al. (1996b) state that the particle approaches the vorticity axis, this behaviour is not shown. Moreover, in the time series plot of angle with the vorticity axis in the results of Iso et al. (1996b) a deviation from zero and a small oscillatory behaviour can be observed, indicating the possibility of the existence of a stable limit cycle near the vorticity axis in their experiments. The high shear rate experiments of Bartram et al. (1975) were performed by releasing a $\kappa =9.1$ rod close to the gradient axis. The particle initially travels nearly parallel to the FGP and approaches the flow direction. Upon perturbing, it moves out of the FGP along the flow–vorticity plane in a monotonic fashion before it starts spiralling near the vorticity axis. This is qualitatively similar to the orientation dynamics in the $R_{vort}^{(2)}$ region shown in figure 2 where a saddle point exists near the flow direction (stable manifold in the FGP) and a stable limit cycle is near vorticity axis. Our theory is, however, unable to fully explain the high shear rate experiments of Iso et al. (1996b) performed for $\kappa =34.4$, $De=3$ and $c=0.39$. Iso et al. (1996b) observe the particle to remain at an angle between $5^\circ$ to $60^\circ$ from the vorticity axis near the flow–vorticity plane after spiralling or moving monotonically away from the FGP. The regions $R_{flow-vort}^{(2)}$, $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ (shown in figure 2) do have an attractor near the flow–vorticity plane, but it is a saddle point. From (5.24) ${p}_{3,saddle}=0.7$ for $c\, De=1.17$ and $\kappa =34.4$, i.e. the saddle point is at an angle of $25.2^\circ$ from vorticity axis. In our theory, the absence of a stable fixed point in the flow–vorticity plane between the vorticity and flow directions may be due to neglecting shear thinning, finite $c$ effects, finite polymer length and polymer entanglement. We use an Oldroyd-B equation to model the polymer stress. Finite polymer length may be captured by the FENE-P model and polymer entanglement likely at larger $c$ by the Giesekus model Bird et al. (1987). Future numerical investigations may be used to test our theory and further clarify the previous experimental findings.

Our findings have a major implication in using dilute (low $c$) polymeric liquids to achieve desired properties such as strength and anisotropy of products manufactured from dilute fibre-filled suspensions mentioned in the introduction. At very low $c\, De$, all the fibres will eventually have orientations close to the vorticity axis. Therefore, adding a small polymer concentration to the fluid in roll-to-roll manufacturing can lead to low resistance films with higher anisotropy and hence better quality. Since the limit cycle in the low $c\, De$ regime becomes closer to the vorticity axis as $c\, De$ increases, the anisotropy can be tuned by changing the shear rate or polymer relaxation time ($De$ is the product of shear rate and polymer relaxation time). Even higher anisotropy can be obtained if very large $De$ can be achieved since at a large $c\, De$, the fibres with all initial orientations ultimately align near the flow direction. For moderate values of $c\, De$, the flow field that precedes a period of simple shear may preorient the fibres somewhat closer to the FGP or the vorticity axis. This will determine whether the fibres lie in the basin of attraction of attractor near the flow direction or the attractor near the vorticity axis. The polymer stresses in the simple shear flow can then drive the particles to a single final orientation leading to highly aligned fibres.