3.1. Heavy particle in creeping shear flow

Jeffery (Reference Jeffery1922) derived analytical expressions for the torque on an ellipsoidal particle in a linear Stokes flow ( $\mathit{Re}_{p}=0$ ) for any angular velocity vector. By coupling this torque to the equation of motion for the particle, a differential equation could be obtained and integrated, which was done numerically by Lundell & Carlsson (Reference Lundell and Carlsson2010). This differential equation in non-dimensional form has a dependency on the Stokes number $\mathit{St}={\it\alpha}\,\mathit{Re}_{p}$ , which determines the effect of particle inertia relative to viscous forces. Even though the torque was derived for $\mathit{Re}_{p}=0$ , a solution could be found for finite $\mathit{St}$ . The Stokes flow solution can, however, be considered valid even in a low- $\mathit{Re}_{p}$ regime. As mentioned in the introduction, Lundell & Carlsson (Reference Lundell and Carlsson2010) and Nilsen & Andersson (Reference Nilsen and Andersson2013) reported a transition for elongated heavy spheroids in linear shear flow from tumbling to rotating, due to particle inertial forces dominating viscous damping forces. The period is found to be a decreasing function of the Stokes number, with the period $GT_{J}=2{\rm\pi}(r_{p}^{-1}+r_{p})$ obtained at $\mathit{St}=0$ and the period $GT_{H}=4{\rm\pi}$ as $\mathit{St}\rightarrow \infty$ . The transition is characterized by $\mathit{St}_{0.5}$ , which is defined as the Stokes number at which the spheroid rotates with a period of $(GT_{J}+GT_{H})/2$ . For a particle of $r_{p}=4$ , this number is $\mathit{St}_{0.5}=307$ , and was observed by Lundell & Carlsson (Reference Lundell and Carlsson2010) to be an increasing function of $r_{p}$ . It should be mentioned that for a particle with this quite low value of $r_{p}$ , the transition from tumbling to rotating will be rather smooth since $GT_{J}=8.5{\rm\pi}$ , which is not much larger than $GT_{H}=4{\rm\pi}$ . For a more slender particle the transition is sharper, as was also observed by Lundell & Carlsson (Reference Lundell and Carlsson2010) and Nilsen & Andersson (Reference Nilsen and Andersson2013).

The transition is illustrated in figure 3, where the angular velocity in the final rotational state is plotted against the angle at different $\mathit{St}$ values. In the figure, a comparison is also made with the present LB-EBF model, where a prolate spheroid with $r_{p}=4~(a=12,b=3)$ is simulated for $\mathit{Re}_{p}=0.5$ (with shear rate lowered to $G=1/34\,560$ to ensure low $\mathit{Kn}$ ) in a computational domain of $N=120$ . The density ratio was chosen to be ${\it\alpha}=2\times 10^{0},2\times 10^{1},\dots ,2\times 10^{5}$ , leading to the same range of Stokes numbers as in the analytical case. The simulated results are compared with the analytical results at the same Stokes number. The simulated results agree well with the analytical ones, with some small discrepancies which can be explained by the slight non-zero $\mathit{Re}_{p}$ in the simulations.

Figure 3. Phase diagrams showing the motion of a spheroid with $r_{p}=4$ in a flow with negligible fluid inertia: our numerical model at $\mathit{Re}_{p}=0.5$ (circles) is compared with analytical results from Lundell & Carlsson (Reference Lundell and Carlsson2010) (solid lines) for: (a)  $\mathit{St}=1$ ; (b)  $\mathit{St}=10$ ; (c)  $\mathit{St}=100$ ; (d)  $\mathit{St}=1000$ ; (e)  $\mathit{St}=10\,000$ ; (f)  $\mathit{St}=100\,000$ . Panels (a)–(f) show the transitions from tumbling to rotating; in each panel the dashed line indicates the location of $\dot{{\it\phi}}_{min}$ during a period.

3.2. Neutrally buoyant particle ( ${\it\alpha}=1,r_{p}=4$ ) at $\mathit{Re}_{p}>0$

Rosén et al. (Reference Rosén, Lundell and Aidun2014) investigated a prolate spheroid with $r_{p}=4$ and reported that the particle inertia and fluid inertia are always in competition, even for a neutrally buoyant particle. As fluid inertia, associated with $\mathit{Re}_{p}$ , increases, particle inertia also increases since it is associated with the Stokes number $\mathit{St}=1\cdot \mathit{Re}_{p}$ . Therefore, if the particle is released at an orientation where the symmetry axis gains enough angular momentum, centrifugal forces will take it to a planar tumbling motion. Otherwise, fluid inertia effects will dominate, leading to non-planar motion (i.e. log-rolling, inclined rolling, inclined kayaking or kayaking). At $\mathit{Re}_{p}>\mathit{Re}_{c}$ , fluid inertia can be strong enough to completely stop the motion, and the particle will end up in a steady state. The transitions for a neutrally buoyant particle with $r_{p}=4$ , as found by Rosén et al. (Reference Rosén, Lundell and Aidun2014), are summarized in table 2, and the dynamical behaviour at $\mathit{Re}_{p}=60$ (T/LR), 70 (T/IR), 74 (T/IK) and 80 (T) is shown in figure 4.

Table 2. Summary of all dynamical transitions of a neutrally buoyant spheroidal particle with $r_{p}=4$ in a linear shear flow, as found by Rosén et al. (Reference Rosén, Lundell and Aidun2014).

Figure 4. Trajectories of a neutrally buoyant particle with $r_{p}=4$ initialized at rest from five different initial orientations: (a) particle motion projected onto the $xy$ -plane by $(p_{x},p_{y})=(\sin {\it\theta}\cos {\it\phi},\sin {\it\theta}\sin {\it\phi})$ ; (b) the five initial orientations according to (3.1) projected onto the $xy$ -plane, where the grey circles indicate the initial conditions, which due to symmetry in the flow problem are equivalent to the initial orientations 2–5, and the large circle represents orientations in the flow-gradient plane (when the particle is in planar rotation at ${\it\theta}={\rm\pi}/2$ ); projected trajectories from simulations at (c)  $\mathit{Re}_{p}=60$ , (d)  $\mathit{Re}_{p}=70$ , (e)  $\mathit{Re}_{p}=74$ and (f)  $\mathit{Re}_{p}=80$ , where the empty circles represent the initial orientations and the filled circles the orientations at $G\cdot t=1000$ .

In the subcritical Hopf bifurcation at $\mathit{Re}_{\mathit{LR}}$ , the log-rolling solution is stabilized, and there are two coexisting stable rotational states, one of which is associated with fluid inertial effects (log-rolling) and the other with particle inertial effects (tumbling). An unstable limit cycle (sketched as a dashed curve in figure 4 b,c) is created in this transition. This limit cycle is located between the two solutions and distinguishes the initial orientations that lead to each rotational state. To illustrate this, a particle is simulated at $\mathit{Re}_{p}=60,70,74$ and 80 with five different initial orientations according to

(3.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}({\it\phi}_{0},{\it\theta}_{0})_{1}=(0,0),\\ \displaystyle ({\it\phi}_{0},{\it\theta}_{0})_{2{-}5}=\left(\frac{{\rm\pi}n}{2},\frac{{\rm\pi}m}{8}\right),\quad n=0,1~\text{and}~m=1,2.\end{array}\right\}\end{eqnarray}$$

In figure 4(a) the projection $(p_{x},p_{y})=(\sin {\it\theta}\cos {\it\phi},\sin {\it\theta}\sin {\it\phi})$ of a simulated particle trajectory is displayed, and figure 4(b) shows the projection of the five initial orientations from (3.1). Figure 4(c) shows that at $\mathit{Re}_{p}=60$ , the initial orientations 1, 2 and 4 are inside the unstable limit cycle and lead to log-rolling, whereas the initial orientations 3 and 5 are outside the limit cycle and end up in the tumbling state. The size of the limit cycle thus determines how probable it is that a random orientational configuration of particles, initialized at rest, will end up in the log-rolling state. After the supercritical pitchfork bifurcation at $\mathit{Re}_{\mathit{PF}}$ , the unstable limit cycle shrinks and reshapes, as can be seen in figure 4(d,e); finally, the limit cycle is annihilated in a saddle-node bifurcation of limit cycles at $\mathit{Re}_{T}$ , as seen in figure 4(f).

The transition at $\mathit{Re}_{p}=\mathit{Re}_{c}$ , for a particle of $r_{p}=4$ , is through an infinite-period saddle-node bifurcation with a diverging tumbling period close to $\mathit{Re}_{p}=\mathit{Re}_{c}$ (Ding & Aidun Reference Ding and Aidun2000; Rosén et al. Reference Rosén, Lundell and Aidun2014). The bifurcation gives rise to two equilibrium orientations, where there is zero torque on the particle: one stable orientation ${\it\phi}_{c}$ and one unstable orientation ${\it\phi}_{us}$ ; see figure 5. This means that for a small deviation ${\it\phi}_{c}\pm {\it\delta}{\it\phi}$ or ${\it\phi}_{us}+{\it\delta}{\it\phi}$ , the particle will go to ${\it\phi}_{c}$ , but a small negative deviation ${\it\phi}_{us}-{\it\delta}{\it\phi}$ will make the particle rotate around the vorticity direction until it is almost aligned with the flow again at ${\it\phi}_{c}-{\rm\pi}$ . The transition at $\mathit{Re}_{c}$ is based on the balance between the torque on the particle from the primary flow and the opposing torque from the secondary flow, as explained in Ding & Aidun (Reference Ding and Aidun2000). The magnitude of the torque from the secondary flow increases with the Reynolds number, and at the critical value it balances the torque from the primary flow. The quantitative value of $\mathit{Re}_{c}$ predicted numerically depends on the exact point of balance reached. This makes the value of $\mathit{Re}_{c}$ strongly dependent on precise prediction of the flow around the particle. This is true also for determination of the exact values of $\mathit{Re}_{PF},\mathit{Re}_{\mathit{Hopf}}$ and $\mathit{Re}_{T}$ .

Figure 5. Illustration of equilibrium angles at $\mathit{Re}_{p}>\mathit{Re}_{c}$ : (a) stable equilibrium angle at ${\it\phi}_{c}$ ; (b) unstable equilibrium angle at ${\it\phi}_{us}$ .

3.3. Bifurcation diagram

In practical applications, the solid-to-fluid density ratio of the particle and fluid rarely varies; instead, it is of greater significance what the rotational states of the particles are with constant ${\it\alpha}$ and varying $\mathit{Re}_{p}$ . This is best illustrated with a bifurcation diagram, shown in figure 6, following Rosén et al. (Reference Rosén, Lundell and Aidun2014).

To draw bifurcation diagrams, a norm is needed to distinguish between the different rotational states. Rosén et al. (Reference Rosén, Lundell and Aidun2014) used the measure $N_{C}=C_{norm}(1+|{\bf\omega}/G|)$ , where $C_{norm}$ is an orbit parameter defined as

(3.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle C_{norm}=\frac{C}{C+1},\\ C=r_{p}^{-1}\tan {\it\theta}(r_{p}^{2}\sin ^{2}{\it\phi}+\cos ^{2}{\it\phi})^{1/2}.\end{array}\right\}\end{eqnarray}$$

This yields:

1. (i) $N_{C}=\text{const.}$

   1. (a) $N_{C}=0$ , log-rolling;

   2. (b) $0<N_{C}<1$ , inclined rolling;

   3. (c) $N_{C}=1$ , steady state;

2. (ii) $N_{C}=f(t)$

   1. (a) $N_{C,mean}\geqslant 0.5$ , tumbling/rotating;

   2. (b) $N_{C,mean}<0.5$ , inclined kayaking, kayaking.

The definition above also makes it easy to distinguish planar rotational states ( $N_{C}\geqslant 1$ ) from non-planar rotational states ( $N_{C}<1$ ). The oscillating states for which $N_{C}=f(t)$ are shown by plotting $N_{C,max}$ and $N_{C,min}$ during a period, represented by the grey area in figure 6. The dashed lines illustrate the approximate location of unstable fixed points, while the open circles illustrate the approximate location of $N_{C,max}$ of the unstable cycles. The bifurcation diagram reveals the following hysteresis behaviour: if one starts on the non-planar branch (where $N_{C}<1$ ) and increases or decreases $\mathit{Re}_{p}$ , the particle will eventually end up in a tumbling motion for $\mathit{Re}_{p}>\mathit{Re}_{T}$ or $\mathit{Re}_{p}<\mathit{Re}_{\mathit{LR}}$ ; however, since the tumbling solution is also available at $\mathit{Re}_{\mathit{LR}}<\mathit{Re}_{p}<\mathit{Re}_{T}$ , if $\mathit{Re}_{p}$ is now varied the particle will always stay in this planar motion.

Figure 6. Bifurcation diagram for a spheroid with $r_{p}=4$ ; the zoomed-in region with $\mathit{Re}_{p}$ between 60 and 100 is displayed on the right. The stable rotational states have either a fixed orientation (solid lines) or an oscillating orientation (grey area), where in the latter case the norm is plotted using $N_{C,max}$ and $N_{C,min}$ during an oscillation period; the unstable fixed points are schematically drawn with dashed lines, and the approximate location of $N_{C,max}$ for the unstable limit cycle is indicated by the open circles.

3.4. Motivation for the present work

The knowledge from previous work can be summarized in $\mathit{Re}_{p}$ – $\mathit{St}$ parameter space, as is done in figure 7 for a prolate spheroid with $r_{p}=4$ . The entire regime of $\mathit{St}>0$ and $\mathit{Re}_{p}=0$ is well understood from the work of Lundell & Carlsson (Reference Lundell and Carlsson2010), and the line where $\mathit{St}=\mathit{Re}_{p}$ , i.e.  ${\it\alpha}=1$ , is well understood from the work of Rosén et al. (Reference Rosén, Lundell and Aidun2014) and others.

Figure 7. Dynamical transitions of a prolate spheroidal particle with $r_{p}=4$ in the $\mathit{Re}_{p}$ – $\mathit{St}$ parameter space, as known from previous work by Lundell & Carlsson (Reference Lundell and Carlsson2010) and Rosén et al. (Reference Rosén, Lundell and Aidun2014). As in table 1, $\text{T}=\text{tumbling}$ , $\text{R}=\text{rotating}$ , $\text{LR}=\text{log-rolling}$ , $\text{IR}=\text{inclined rolling}$ , $\text{IK}=\text{inclined kayaking}$ , $\text{K}=\text{kayaking}$ and $\text{S}=\text{steady state}$ .

The region where $\mathit{St}>\mathit{Re}_{p}$ is not well documented. The critical Reynolds number $\mathit{Re}_{c}$ is defined for a stationary particle and hence is dependent only on the geometry of the problem; this number therefore cannot depend on $\mathit{St}$ . The dependence of the other critical Reynolds numbers ( $\mathit{Re}_{\mathit{LR}}$ , $\mathit{Re}_{PF}$ , $\mathit{Re}_{\mathit{Hopf}}$ and $\mathit{Re}_{T}$ ) on the Stokes number was not known prior to this work. Furthermore, the dependence of $\mathit{St}_{0.5}$ on $\mathit{Re}_{p}$ and how all these transitions depend on the aspect ratio $r_{p}$ are also not clear. All of these questions will be addressed in this work, where a complete picture of the dynamics of a single prolate spheroid will be provided. The results help to explain contradictions in the literature and provide a complete framework for analysing rotational motion of a particle in shear flow.

4.1. Comments about the sensitivity of $\mathit{Re}_{c}$ to the flow field near the particle

The transition at $\mathit{Re}_{c}$ is very sensitive to small changes in the fluid field around the particle, as it is caused by a very delicate balance of torques on the particle. At $\mathit{Re}_{p}=\mathit{Re}_{c}$ , the negative torque $T^{-}$ from the primary flow is exactly equal to the positive torque $T^{+}$ from the secondary flow (illustrated in figure 8). Therefore, a small error of $\pm 1\,\%$ in determining $T^{+}$ or $T^{-}$ can lead to a large error in the value of $\mathit{Re}_{c}$ . This stringent requirement of accuracy in determining the torque results in difficulty pinpointing the critical Reynolds number using any numerical method without extreme grid refinement around the particle.

Figure 8. Qualitative flow field close to steady state at $\mathit{Re}_{p}\approx \mathit{Re}_{c}$ ; the primary flow (1) gives a negative torque (clockwise direction in left panel) on the particle, and the secondary flow (2) gives a positive torque (anticlockwise direction in left panel) on the particle (Ding & Aidun Reference Ding and Aidun2000).

The aim of the present work is to describe dynamical states of the particle in shear flow, and a moving grid is therefore crucial. The qualitative changes in the fluid field as $\mathit{Re}_{p}$ changes are not sensitive to resolution, so the order in which the dynamical states occur and the particle behaviour in each dynamical region are resolution-independent. This is probably why the sensitivity described above has not been reported before. For example, the present model matches the transitional Reynolds numbers of Huang et al. (Reference Huang, Yang, Krafczyk and Lu2012b ) when the same resolution is used (Rosén et al. Reference Rosén, Lundell and Aidun2014). The critical Reynolds number $\mathit{Re}_{c}$ was estimated by Huang et al. (Reference Huang, Yang, Krafczyk and Lu2012b ) to be $\mathit{Re}_{c}=445\pm 5$ for the case of $r_{p}=2$ and ${\it\kappa}=0.5$ .

Special consideration must be given to two parameters in lattice Boltzmann methods, due to the kinetic basis of the approach. One consideration is that the maximum Knudsen number $\mathit{Kn}_{max}$ must be small enough in order to correctly describe a continuum fluid. This number is defined close to $\mathit{Re}_{c}$ as $\mathit{Kn}_{max}=2\sqrt{3}\,a\,Gr_{p}/\mathit{Re}_{c}$ , and can be seen in figure 9(b) to be small (approximately $O(10^{-3})$ ) and decreasing with $N$ . These effects are therefore considered negligible. The second consideration is that the maximum lattice Boltzmann Mach number in the simulation, $Ma_{max}$ , should be small in order to eliminate any compressible effects in the fluid. This number is linearly increasing with $N$ via the relation $Ma_{max}=\sqrt{3}GN/2$ . Therefore, to keep $Ma_{max}$ small, as $N$ increases $G$ must become smaller to keep $GN$ small. The smaller the value of $G$ , the larger the number of time steps required to reach the equilibrium state. To pinpoint the quantitative value of $\mathit{Re}_{c}$ requires very large $N$ , very small $G$ and very long computational time. However, in order to capture the dynamics of the particle motion and the equilibrium states, it is not necessary to pinpoint the value of $\mathit{Re}_{c}$ .

Figure 9. Resolution dependence of the critical Reynolds number $\mathit{Re}_{c}$ for a prolate spheroid with $r_{p}=4~(a=N/10)$ and ${\it\kappa}=0.2$ : (a)  $\mathit{Re}_{c}$ versus the spatial resolution $N$ with constant shear rate $G=1/2400$ or $G=1/1800$ ; (b)  $Ma_{max}$ and $\mathit{Kn}_{max}$ versus the spatial resolution  $N$ .

To investigate the desired range of $\mathit{Re}_{p}$ , the shear rate was set to $G=1/600$ , which was seen by Rosén et al. (Reference Rosén, Lundell and Aidun2014) to give correct behaviour of particle motion. However, the physical transitional $\mathit{Re}_{p}$ values, for example $\mathit{Re}_{c}$ , could differ from the values presented here, since $G$ is not small enough (and $N$ not large enough) for the solution to have converged to the precise value of $\mathit{Re}_{c}$ (Rosén et al. Reference Rosén, Lundell and Aidun2014). Figure 9 shows that the value of $\mathit{Re}_{p}$ converges to the accurate $\mathit{Re}_{c}$ as the resolution $N$ increases (with particle size determined by $a=N/10$ ). It is found that, due to the slow convergence with respect to the spatial resolution, getting a converged value of $\mathit{Re}_{c}$ will be extremely computationally expensive; the values obtained are in the range of 69–132.

Consequently, to obtain a converged solution of $\mathit{Re}_{c}$ not only requires high spatial resolution through a high value of $N$ but also a low value of the shear rate $G$ in order to be able to neglect compressible effects. Because the lattice Boltzmann equation is explicit in time integration, the shear rate $G$ is simultaneously connected to the time resolution of the particle rotation, and more time steps are required for lower values of $G$ . Furthermore, as $\mathit{Re}_{p}$ approaches the critical point, i.e.  $\mathit{Re}_{c}$ , the period of oscillation increases to infinity asymptotically. In essence, the computational cost of achieving a simulation that is fully converged in space with a low $Ma_{max}$ is huge.

However, a ‘steady state’ approach can be employed to accurately pinpoint $\mathit{Re}_{c}$ , as shown by Ding & Aidun (Reference Ding and Aidun2000). Determining the value of $\mathit{Re}_{c}$ does not require performing a transient simulation, since it is defined as the minimum $\mathit{Re}_{p}$ at which a stationary particle with a certain orientation ${\it\phi}~({\it\theta}={\rm\pi}/2)$ is experiencing zero torque. Therefore, one can fix the orientation of the particle, compute steady-state flow around the particle, and then change the orientation until the torque is balanced at the minimum value of $\mathit{Re}_{p}$ , which would be $\mathit{Re}_{c}$ . This approach was demonstrated by Ding & Aidun (Reference Ding and Aidun2000) with the lattice Boltzmann method. However, because in this case the flow is steady state, one can solve the steady-state Navier–Stokes equation to pinpoint the value of $\mathit{Re}_{c}$ without any need for time integration. This was done by setting up the problem in Comsol Multiphysics 4.3a. The same flow problem was considered, using a fixed particle and evaluating the torque for different particle orientations ${\it\phi}$ and $\mathit{Re}_{p}$ . The minimum $\mathit{Re}_{p}$ at which a zero torque exists is found to be $\mathit{Re}_{c}=131\pm 5$ , and the corresponding orientation is ${\it\phi}=6.6^{\circ }\pm 0.5^{\circ }$ . This value seems reasonable, considering the LB-EBF data in figure 9. Using the same Comsol Multiphysics model, an estimate of $\mathit{Re}_{c}=476\pm 10$ is obtained for the case of $r_{p}=2$ and ${\it\kappa}=0.5$ , as opposed to the lattice Boltzmann simulations which gave $\mathit{Re}_{c}=445\pm 5$ .

To conclude this section, we observe that the same sequence of transitions between dynamical states as $\mathit{Re}_{p}$ and $\mathit{St}$ are changed is found regardless of the resolution (as long as the particle is large enough, as discussed in § 2.1). By fixing the spatial resolution and shear rate to constant values, as is done here, the resolution-dependence can be neglected and the transitional $\mathit{Re}_{p}$ can be compared between the simulations. However, in order to pinpoint the physical transitional Reynolds numbers, other methods must be employed. At present this is feasible in cases where the particle is steady (such as $\mathit{Re}_{c}$ ), but not for transitions from one moving state to another.

4.2. Prediction of $\mathit{St}_{c}$

In the $\mathit{Re}_{p}=0$ case, the spheroid in the tumbling motion has a minimum angular velocity $\dot{{\it\phi}}_{min}$ , which increases with $\mathit{St}$ as can be deduced from figure 3. On the other hand, for a neutrally buoyant particle, $\dot{{\it\phi}}_{min}$ decreases as $\mathit{Re}_{p}$ increases from 10 to 50 (figure 10) and is the source of the infinite-period saddle-node bifurcation, which occurs at $\mathit{Re}_{p}=\mathit{Re}_{c}\approx 89$ , as $\dot{{\it\phi}}_{min}$ becomes zero somewhere on the orbit. When $\mathit{St}$ is increased to 1000 at $\mathit{Re}_{p}=100$ , an orbit with $\dot{{\it\phi}}_{min}>0$ can again be found. Since the steady state must still exist, as it is dependent only on geometry, two solutions (tumbling and steady state) will exist above some critical Stokes number $\mathit{St}_{c}$ . Hence, a heavy particle can maintain a tumbling motion due to particle inertia given sufficient initial angular momentum. If the particle is not initialized with enough angular momentum, the fluid inertia will dominate and take the particle to a steady state.

Figure 10. Phase diagrams of a tumbling particle with $r_{p}=4$ , for: (a)  $\mathit{St}=10$ ; (b)  $\mathit{St}=50$ ; (c)  $\mathit{St}=100$ ; (d)  $\mathit{St}=1000$ . In each panel the dashed line indicates the location of $\dot{{\it\phi}}_{min}$ , which does not exist in panel (c) (where $\mathit{Re}_{p}=\mathit{St}=100$ ), since there is no stable tumbling solution; the steady state and the tumbling solution can coexist at $\mathit{Re}_{p}=100$ if the Stokes number is increased to $\mathit{St}=1000$ , as in panel (d).

4.3. Rotational states

The results from the present numerical study for a prolate spheroid of $r_{p}=4$ can be summarized in a single state plot as shown in figure 11. The figure shows the different rotational states that can be present for a certain pair of $\mathit{Re}_{p}$ and $\mathit{St}$ values. At every $\mathit{Re}_{p}$ , there is a periodic solution at sufficiently high $\mathit{St}$ where the particle rotates with axis of symmetry in the flow-gradient plane. This solution is tumbling in the moderate $\mathit{St}$ regime, but at high $\mathit{St}$ there is a transition to rotating motion. Since this solution is associated with high particle inertia, it will most likely be the preferred rotational state if the particle is initialized with sufficient angular momentum. Figure 11 shows thirteen distinct regimes depending on $\mathit{Re}_{p}$ , which are listed below (an impressive number indeed for such a seemingly simple case). The details on how the transitions are determined will be given in the next subsection.

Figure 11. State plot showing the different rotational states for a spheroid with $r_{p}=4$ , depending on the values of $\mathit{Re}_{p}$ and $\mathit{St}$ ( $\text{T}=\text{tumbling}$ , $\text{R}=\text{rotating}$ , $\text{LR}=\text{log{-}rolling}$ , $\text{IR}=\text{inclined rolling}$ , $\text{IK}=\text{inclined kayaking}$ , $\text{K}=\text{kayaking}$ , $\text{S}=\text{steady state}$ ); the error bar under the $\mathit{Re}_{p}$ -axis indicates the range in which $\mathit{Re}_{c}$ has been found due to a dependence on resolution, and the black cross indicates the value of $\mathit{Re}_{c}$ obtained from Comsol (see § 4.1).

4.4. Determining $\mathit{St}_{c}$ as function of $\mathit{Re}_{p}$

4.4.2. Equilibrium angles

In order to find the critical density ratio ${\it\alpha}_{c}$ , the equilibrium angle ${\it\phi}_{us}$ must be found. Both ${\it\phi}_{us}$ and ${\it\phi}_{c}$ are found by using the fact that there should be zero torque on the particle, which means that $\dot{{\it\phi}}\approx 0$ should be achieved shortly after the particle is initialized in an equilibrium angle. This fact is used to systematically find the equilibrium angles for a given $\mathit{Re}_{p}$ with a half-interval search (the exact algorithm is described in appendix B).

In figure 12, ${\it\phi}_{c}$ and ${\it\phi}_{us}$ are plotted as functions of $\mathit{Re}_{p}$ . It can be seen that ${\it\phi}_{us}$ and ${\it\phi}_{c}$ are decreasing and increasing functions of $\mathit{Re}_{p}$ , respectively. A seventh-order polynomial is fitted to the computed data points, and using the polynomial function the bifurcation is found to occur at $\mathit{Re}_{c}\approx 90$ , which is slightly higher than the previously obtained value of $\mathit{Re}_{c}\approx 89$ .

Figure 12. Equilibrium angles ${\it\phi}_{c}$ (circles) and ${\it\phi}_{us}$ (crosses) plotted as functions of $\mathit{Re}_{p}$ for a spheroid with $r_{p}=4$ ; the data are fitted with a polynomial function (solid line), for which the minimum is found at $\mathit{Re}_{p}=\mathit{Re}_{c}=90$ .

4.4.3. Critical density ratio as a function of $\mathit{Re}_{p}$

The critical density ratio, ${\it\alpha}_{c}$ , is obtained as described in § 4.4.1, with ${\it\delta}{\it\phi}$ chosen sufficiently large so that the particle is definitely moving with $\dot{{\it\phi}}<0$ . In this study, ${\it\delta}{\it\phi}$ was chosen to be $1^{\circ }$ . The critical density ratio ${\it\alpha}_{c}$ is found by using the fact that the particle with ${\it\alpha}={\it\alpha}_{c}$ initialized at rest from ${\it\phi}_{0}={\it\phi}_{us}-{\it\delta}{\it\phi}$ should reach ${\it\phi}={\it\phi}_{0}-{\rm\pi}$ with $\dot{{\it\phi}}=0$ . The critical density ratio is found by using half-interval search (the exact algorithm is provided in appendix C).

Figure 13 shows ${\it\alpha}_{c}$ as a function of $\mathit{Re}_{p}$ , where the data points are fitted with a rational function of the form ${\it\alpha}_{c}=P(\mathit{Re}_{p})/Q(\mathit{Re}_{p})$ , with $P$ and $Q$ being fifth- and fourth-degree polynomials, respectively. The error bars indicate the minimum intervals ${\rm\Delta}\mathit{St}=10$ used in the half-interval algorithm. The maximum critical density ratio of ${\it\alpha}_{c,max}\approx 2.3$ is found close to $\mathit{Re}_{p}\approx 125$ , while at higher $\mathit{Re}_{p}$ particle inertia plays a greater role and ${\it\alpha}$ does not need to be much larger than unity in order to find a periodic state. It can also be seen from the figure that ${\it\alpha}_{c}>1$ as $\mathit{Re}_{p}=\mathit{Re}_{c}$ . This codimension-two point ${\it\alpha}_{cd2}={\it\alpha}_{c}(\mathit{Re}_{c})\approx 1.6$ , where two bifurcations coexist (the saddle-node bifurcaction at $\mathit{Re}_{c}$ and the homoclinic bifurcation at ${\it\alpha}_{c}$ ), is found by evaluating the fitted rational function at $\mathit{Re}_{c}$ . A more detailed discussion of this bifurcation will be given in § 5. The line of $\mathit{St}_{c}$ in figure 11 is drawn by the relation $\mathit{St}_{c}={\it\alpha}_{c}\,\mathit{Re}_{p}$ , which connects to the codimension-two point at $\mathit{St}_{cd2}={\it\alpha}_{cd2}\mathit{Re}_{c}$ .

Figure 13. The critical density ratio ${\it\alpha}_{c}$ plotted as a function of $\mathit{Re}_{p}$ for a spheroid with $r_{p}=4$ (solid blue line), fitted with a rational function (solid black line). The critical Reynolds number $\mathit{Re}_{c}=90$ (dashed line) is found through polynomial fitting of equilibrium angles, and the codimension-two point, where ${\it\alpha}_{c}$ connects to the $\mathit{Re}_{c}$ line, is found at ${\it\alpha}={\it\alpha}_{cd2}\approx 1.6$ . Black crosses indicate ${\it\alpha}=2,2.5$ and 3 at $\mathit{Re}_{p}=150$ , for which the phase diagrams in figure 14 are drawn. Only planar rotational states, namely tumbling (T) and steady state (S), are shown in the figure.

4.5. Determining $\mathit{St}_{0.5}$ as function of $\mathit{Re}_{p}$

As already discussed, it was observed by Lundell & Carlsson (Reference Lundell and Carlsson2010) and Nilsen & Andersson (Reference Nilsen and Andersson2013) that the period of the tumbling motion is a steadily decreasing function of $\mathit{St}$ in the absence of fluid inertia ( $\mathit{Re}_{p}=0$ ). At $\mathit{St}=0$ , the particle is tumbling with period $GT_{J}=2{\rm\pi}(r_{p}^{-1}+r_{p})$ ; as $\mathit{St}\rightarrow \infty$ , the particle rotates at a constant angular velocity with period $GT_{H}=4{\rm\pi}$ . The transition between the two regimes is characterized by $\mathit{St}_{0.5}$ , which is defined as the Stokes number at which the spheroid rotates with a period of $(GT_{J}+GT_{H})/2$ . However, $\mathit{St}_{0.5}$ was defined only for $\mathit{Re}_{p}=0$ , since it refers to the period without any fluid or particle inertia. For finite $\mathit{Re}_{p}>0$ , the period is also seen to decrease as $\mathit{St}$ increases; therefore in this work we use the same definition of $\mathit{St}_{0.5}$ , i.e. the Stokes number at which the period is $GT=(GT_{J}+GT_{H})/2={\rm\pi}(r_{p}^{-1}+r_{p}+2)$ . The value of $\mathit{St}_{0.5}$ is found by using a half-interval search (see appendix D for the exact algorithm). Since the value at $\mathit{Re}_{p}=0$ can be analytically found to be $\mathit{St}_{0.5}=307$ , the line with the simulated values has been connected to this analytical value in figure 11.

4.6. Determining $\mathit{Re}_{\mathit{LR}},\mathit{Re}_{PF},\mathit{Re}_{\mathit{Hopf}}$ and $\mathit{Re}_{T}$ as functions of $\mathit{St}$

5.1. Bifurcation at ${\it\alpha}={\it\alpha}_{c}~(St=\mathit{St}_{c})$

From figure 13 it can be seen that for a particle of ${\it\alpha}>{\it\alpha}_{c,max}$ , fluid inertia will never be enough to fully stop a tumbling motion, and there will be no diverging period when $\mathit{Re}_{p}$ increases or decreases (a steady state will still always coexist at $\mathit{Re}_{p}>\mathit{Re}_{c}$ ). However, by fixing $\mathit{Re}_{p}>\mathit{Re}_{c}$ and decreasing ${\it\alpha}$ , the period diverges at ${\it\alpha}={\it\alpha}_{c}$ . This transition is not the infinite-period saddle-node bifurcation seen at $\mathit{Re}_{p}=\mathit{Re}_{c}$ , since no new fixed points are created. To understand this process, phase diagrams of the motion are drawn in figure 14(b–d), where a particle is initialized at rest from eight different initial orientations, $({\it\phi},{\it\theta})_{1{-}8}=({\rm\pi}n/8,{\rm\pi}/2)$ for $n=0,1,\dots ,7$ (see figure 14 a), at $\mathit{Re}_{p}=150$ and three different density ratios ${\it\alpha}=2,2.5$ and 3 (crosses in figure 13). As shown in figure 12, the equilibrium angles at $\mathit{Re}_{p}=150$ are ${\it\phi}_{us}\approx 0^{\circ }$ and ${\it\phi}_{c}\approx 15.5^{\circ }$ .

Figure 14. Phase diagrams showing homoclinic bifurcation for a spheroid with $r_{p}=4$ at $\mathit{Re}_{p}=150$ , initialized at rest in the flow-gradient plane ( ${\it\theta}=90^{\circ }$ ) from eight different initial orientations, shown in (a), and with three different density ratios: (b)  ${\it\alpha}=3$ ; (c)  ${\it\alpha}=2.5$ ; (d)  ${\it\alpha}=2$ . A limit cycle (tumbling) is seen to collapse with the saddle ${\it\phi}_{us}$ (open circle) at ${\it\alpha}={\it\alpha}_{c}=2.2$ ; at ${\it\alpha}<{\it\alpha}_{c}$ , all initial conditions lead to the sink ${\it\phi}_{c}$ (solid circle) and thus stay in steady state.

Starting with the case where ${\it\alpha}=3$ in figure 14(b), the point at $(\dot{{\it\phi}},{\it\phi})=(0,{\it\phi}_{us})$ is a saddle (open circle, with stable dynamics in one eigendirection and unstable dynamics in another eigendirection) while $(\dot{{\it\phi}},{\it\phi})=(0,{\it\phi}_{c})$ is a sink (filled circle, with stable dynamics, i.e. nearby trajectories go to this point). The trajectory of the initial condition ${\it\phi}_{1}$ almost coincides with the saddle (i.e.  ${\it\phi}_{1}={\it\phi}_{us}$ ) and can be considered a branch of the unstable manifold of the saddle. This manifold asymptotically approaches the stable limit cycle, i.e. a tumbling motion. At some trajectory between the trajectories of initial conditions ${\it\phi}_{6}$ and ${\it\phi}_{7}$ , there is a branch of the stable manifold of the saddle; this means that initial conditions above it (e.g. initial conditions 2–6) lead to the sink, i.e. steady state. As ${\it\alpha}$ is decreased close to ${\it\alpha}_{c}$ , the stable limit cycle goes very near the saddle, as seen in figure 14(c). At ${\it\alpha}={\it\alpha}_{c}\approx 2.2$ , the stable limit cycle is destroyed as the limit cycle merges with the unstable manifold of the saddle. The manifold becomes a homoclinic orbit which connects with itself at the saddle. At ${\it\alpha}<{\it\alpha}_{c}$ , even the unstable manifold of the saddle leads to the sink and therefore the steady state is the only solution (figure 14 d).

The behaviour is very similar to that of a damped pendulum driven by a constant torque (Strogatz Reference Strogatz1994), which undergoes a homoclinic bifurcation when the driving torque is varied under low damping conditions. The period of the limit cycle is divergent, but unlike the case of an infinite-period saddle-node bifurcation, where the period scales as $GT\propto |{\it\beta}-{\it\beta}_{c}|^{-0.5}$ , the period close to the homoclinic bifurcation scales as $GT\propto \ln |{\it\beta}-{\it\beta}_{c}|$ , where ${\it\beta}$ is any parameter leading to the transition at ${\it\beta}_{c}$ . In figure 15(a,b), the period is plotted as a function of ${\it\alpha}-{\it\alpha}_{c}$ and is compared with the two types of scaling at $\mathit{Re}_{p}=150$ . It is clear that the logarithmic scaling gives a much better fit, and the existence of a homoclinic bifurcation at ${\it\alpha}={\it\alpha}_{c}$ is confirmed. The reason that the power law fails to describe the diverging period in this case is that the angular velocity becomes a non-differentiable function of the angle close to the saddle (seen also in figure 13 b), which is an important criterion for this scaling (Strogatz Reference Strogatz1994; Rosén et al. Reference Rosén, Lundell and Aidun2014). The shape of $\dot{{\it\phi}}({\it\phi})$ close to ${\it\phi}_{us}$ is determined by the two eigendirections of the saddle when ${\it\alpha}\gtrsim {\it\alpha}_{c}$ and therefore has a non-parabolic shape.

Figure 15. Scaling of the diverging period at $\mathit{Re}_{p}=150$ for a spheroid with $r_{p}=4$ : (a) data fitted using the assumption of a homoclinic bifurcation with $GT\propto \ln |{\it\beta}-{\it\beta}_{c}|$ ; (b) data fitted using the assumption of an infinite-period saddle-node bifurcation with $GT\propto |{\it\beta}-{\it\beta}_{c}|^{-0.5}$ .

Combining this knowledge with what was found regarding the coexistence of the tumbling solution and a steady state (figure 13), the following conclusions can be drawn for a prolate spheroid of aspect ratio $r_{p}=4$ with increasing $\mathit{Re}_{p}$ .

1. (i) If ${\it\alpha}<{\it\alpha}_{cd2}\approx 1.6$ , the transition from tumbling to steady state will be through an infinite-period saddle-node bifurcation at $\mathit{Re}_{p}=\mathit{Re}_{c}$ , with a diverging period scaling as $GT\propto |\mathit{Re}_{p}-\mathit{Re}_{c}|^{-0.5}$ .

2. (ii) If ${\it\alpha}_{cd2}<{\it\alpha}<{\it\alpha}_{c,max}$ (i.e.  $1.6\lesssim {\it\alpha}\lesssim 2.3$ ), the transition will be through a homoclinic bifurcation at $\mathit{Re}_{p}=\mathit{Re}_{c1}$ , where $\mathit{Re}_{c1}$ is defined as the $\mathit{Re}_{p}$ value at which ${\it\alpha}={\it\alpha}_{c}$ . The tumbling period will diverge according to $GT\propto \ln |\mathit{Re}_{p}-\mathit{Re}_{c1}|$ .

3. (iii) If ${\it\alpha}={\it\alpha}_{cd2}\approx 1.6$ , a transition from tumbling to steady state will occur at $\mathit{Re}_{p}=\mathit{Re}_{c}$ , but no distinction can be made between the types of bifurcations and both types of scalings will apply.

4. (iv) If ${\it\alpha}>{\it\alpha}_{c,max}\approx 2.3$ , there is no diverging period at any $\mathit{Re}_{p}$ , and a tumbling solution will always exist.

5.2. Bifurcation diagrams

The bifurcation diagrams for ${\it\alpha}=2$ and 3 are shown in figure 16(a,b) in the same way as in § 3.3. Looking at figure 16(a), for ${\it\alpha}=2$ , which lies in the interval ${\it\alpha}_{cd2}<{\it\alpha}<{\it\alpha}_{c,max}$ , there are two critical particle Reynolds numbers which correspond to intersections of the ${\it\alpha}_{c}$ line (see figure 13). Due to higher particle inertia, the fluid inertia is not sufficient to stop the tumbling motion at $\mathit{Re}_{p}=\mathit{Re}_{c}$ (note that $\mathit{Re}_{c}$ is dependent only on geometry and fluid properties and is not changed due to high particle inertia). The tumbling period instead diverges to infinity at $\mathit{Re}_{p}=\mathit{Re}_{c1}>\mathit{Re}_{c}$ through a homoclinic bifurcation. Similar behaviour is seen when $\mathit{Re}_{p}$ is decreased from high values, with the tumbling period diverging to infinity at $\mathit{Re}_{c2}$ . In the region $\mathit{Re}_{c1}<\mathit{Re}_{p}<\mathit{Re}_{c2}$ , the steady state is thus the only solution, while a tumbling solution coexists in the regions $\mathit{Re}_{c}<\mathit{Re}_{p}<\mathit{Re}_{c1}$ and $\mathit{Re}_{p}>\mathit{Re}_{c2}$ . For the particle with ${\it\alpha}=3$ (figure 16 b), the density ratio is higher than ${\it\alpha}_{c,max}$ , and hence there is no diverging period and a tumbling solution will always coexist with the steady state at $\mathit{Re}_{p}>\mathit{Re}_{c}$ . The bifurcation diagrams indicate the following hysteresis behaviour as $\mathit{Re}_{p}$ varies.

1. (i) For a particle with ${\it\alpha}_{cd2}<{\it\alpha}<{\it\alpha}_{c,max}$ (figure 16 a) in a tumbling state at $\mathit{Re}_{p}<\mathit{Re}_{c}$ , increasing $\mathit{Re}_{p}$ will lead to a diverging period and a transition to steady state at $\mathit{Re}_{p}=\mathit{Re}_{c1}$ . However, to reach a tumbling state again, $\mathit{Re}_{p}$ must be decreased to $\mathit{Re}_{p}=\mathit{Re}_{c}$ , since the steady-state solution (bold line at $N_{C}=1$ ) exists also for $\mathit{Re}_{c}<\mathit{Re}_{p}<\mathit{Re}_{c1}$ .

2. (ii) For a particle with ${\it\alpha}<{\it\alpha}_{c,max}$ (figure 16 a) in a tumbling state at $\mathit{Re}_{p}>\mathit{Re}_{c2}$ , decreasing $\mathit{Re}_{p}$ will lead to a diverging period and a transition to steady state at $\mathit{Re}_{p}=\mathit{Re}_{c2}$ . Upon increasing $\mathit{Re}_{p}>\mathit{Re}_{c2}$ , the particle will stay in steady state.

3. (iii) For a particle with ${\it\alpha}>{\it\alpha}_{c,max}$ (figure 16 b) in a steady state, decreasing $\mathit{Re}_{p}$ past $\mathit{Re}_{c}$ will make the particle enter a tumbling motion. If $\mathit{Re}_{p}$ is increased again, the particle will stay in a tumbling rotational state.

Figure 16. Bifurcation diagrams for a spheroid of $r_{p}=4$ with: (a)  ${\it\alpha}=2$ ; (b)  ${\it\alpha}=3$ . In each panel, the zoomed-in region for $\mathit{Re}_{p}$ between 60 and 100 is displayed on the right; the stable rotational states have either a fixed orientation (solid lines) or an oscillating orientation (grey areas), where in the latter case the norm is plotted using $N_{C,max}$ and $N_{C,min}$ during an oscillation period. The unstable fixed points are schematically drawn with dashed lines, and the approximate location of $N_{max}$ for the unstable limit cycle is indicated by open circles.

8.1. Dependence of aspect ratio on fluid inertial transitions

In § 4.6, it was seen that the transitions at $\mathit{Re}_{PF}$ and $\mathit{Re}_{c}$ are independent of the density ratio and are caused only by fluid inertia. Therefore the $r_{p}$ -dependence is analysed for these transitions using a neutrally buoyant particle ( ${\it\alpha}=1$ ).

In order to find $\mathit{Re}_{PF}$ , the same procedure as described in § 4.6.2 was performed for aspect ratios $r_{p}=2,2.1,2.2,2.3,2.4,2.6,2.8,3,4,5$ and 6. The critical Reynolds number $\mathit{Re}_{c}$ was found by first finding the equilibrium angles, using the procedure described in § 4.4.2, and then fitting a polynomial function to the data. The fitted polynomials are shown in figure 19. The value of $\mathit{Re}_{c}$ is found by seeking the minimum of the polynomial function. Apart from the aspect ratios used for $\mathit{Re}_{PF}$ , aspect ratios of $r_{p}=3.1,3.2,\dots ,3.9$ were also investigated. When $r_{p}$ is increased, the difference between ${\it\phi}_{us}$ and ${\it\phi}_{c}$ decreases, and at the same time $\mathit{Re}_{c}$ decreases.

Figure 19. Fitted polynomials of the equilibrium angles ${\it\phi}_{c}$ (red solid line) and ${\it\phi}_{us}$ (blue dashed line) as functions of $\mathit{Re}_{p}$ for spheroids with $r_{p}=3.0,3.1,\dots ,3.9,4.0,5.0$ and 6.0; the arrows indicate the effects of increasing $r_{p}$ .

Figure 20 illustrates how $\mathit{Re}_{PF}$ and $\mathit{Re}_{c}$ vary as functions of $r_{p}$ . The transitional Reynolds numbers are found to be decreasing functions of $r_{p}$ , and since the transitions at $\mathit{Re}_{PF}$ and $\mathit{Re}_{c}$ depend only on fluid inertia (i.e. the emerging solutions correspond to steady states of the fluid field), the two critical numbers decay in a similar manner as $r_{p}$ is increased. However, due to the resolution issues discussed in § 4.1, it may not be easy to compare the values because the spatial resolution cannot be held exactly constant when the geometry of the particle is changed. Therefore, in figure 20 the trend of $\mathit{Re}_{c}$ is analysed at higher resolution ( $N=240$ ) and by using Comsol, following the procedure described in § 4.1. It is clear that the critical Reynolds numbers are indeed decreasing with increasing $r_{p}$ . The sensitivity of these transitions indicates that the force balance which causes the transitions at $\mathit{Re}_{PF}$ and $\mathit{Re}_{c}$ becomes even more sensitive to the resolution as $r_{p}$ gets larger.

Figure 20. Plots of $\mathit{Re}_{c}$ and $\mathit{Re}_{PF}$ as functions of $r_{p}$ .

8.2. Dependence of aspect ratio on particle inertial transitions

To study the effect of aspect ratio on ${\it\alpha}_{c}$ , the simulation procedure from § 4.4 was performed for $r_{p}=3{-}6$ . The critical density ratio ${\it\alpha}_{c}$ is plotted as a function of $\mathit{Re}_{p}$ and $r_{p}$ in figure 21. From this figure it is observed that both ${\it\alpha}_{c,max}$ and ${\it\alpha}_{cd2}$ increase rapidly with increasing $r_{p}$ . These dependencies are better illustrated in figure 22, and it seems that ${\it\alpha}_{c,max}$ and ${\it\alpha}_{cd2}$ are linearly dependent on $r_{p}$ . It can also be seen that below a certain aspect ratio, $r_{p}=r_{p,c1}\approx 3.2$ , a critical density ratio of ${\it\alpha}_{c}>1$ does not exist. This means that neutrally buoyant particles with low $r_{p}$ will never experience a diverging period in the tumbling motion, a fact that could also be seen from the results of Huang et al. (Reference Huang, Wu and Lu2012a ,Reference Huang, Yang, Krafczyk and Lu b ). In a range $r_{p,c1}<r_{p}<r_{p,c2}~(3.2\lesssim r_{p}\lesssim 3.7)$ , it is found that ${\it\alpha}_{c,max}>1$ but ${\it\alpha}_{cd2}<1$ . Thus, in this range of aspect ratios, the transition to steady state is through a homoclinic bifurcation at $\mathit{Re}_{p}=\mathit{St}_{c}$ for neutrally buoyant particles.

Figure 21. The critical density ratio ${\it\alpha}_{c}$ plotted as a function of $\mathit{Re}_{p}$ (solid curve) and $\mathit{Re}_{c}$ (dashed line) for a spheroid with $r_{p}=3.5,4,5$ and 6; the arrows show the effect of increasing  $r_{p}$ .

At larger aspect ratios, $r_{p}>r_{p,c2}$ , both ${\it\alpha}_{c,max}$ and ${\it\alpha}_{cd2}$ are greater than $1$ and the transition to steady state will instead be through an infinite-period saddle-node bifurcation at $\mathit{Re}_{p}=\mathit{Re}_{c}$ for neutrally buoyant particles.

The fact that fluid inertia effects dominate over particle inertia effects when $r_{p}\rightarrow \infty$ is also supported by physical reasoning. Particle inertia torques scale with the magnitude of the inertial tensor, which behaves as $1/r_{p}^{2}$ for large $r_{p}$ , whereas torques from fluid stresses on the particle probably scale with the surface area of the particle, behaving as $1/r_{p}$ for large  $r_{p}$ .

Furthermore, the $\mathit{Re}_{p}$ value at which ${\it\alpha}_{c,max}$ occurs seems to get closer to $\mathit{Re}_{c}$ , i.e.  ${\it\alpha}_{cd2}$ approaches ${\it\alpha}_{c,max}$ as $r_{p}$ increases. At $r_{p}=r_{p,c3}\approx 8.0$ they coincide, and this could mean that for more slender particles, no $\mathit{Re}_{c1}$ will exist and a transition from tumbling to steady state for particles with ${\it\alpha}<{\it\alpha}_{cd2}$ will always be through an infinite-period saddle-node bifurcation with increasing $\mathit{Re}_{p}$ . However, the increased sensitivity of determining $\mathit{Re}_{c}$ for more slender particles will result in increased sensitivity for the determination of ${\it\alpha}_{cd2}$ as well. This will probably decrease the slope of ${\it\alpha}_{cd2}$ in figure 22, and thus it cannot be concluded that $r_{p,c3}$ exists.

Figure 22. Plots of the codimension-two point ${\it\alpha}_{cd2}$ (crosses with fitted solid line) and the maximum critical density ratio ${\it\alpha}_{c,max}$ (circles with fitted dashed line) as functions of the particle aspect ratio $r_{p}$ . At a certain $r_{p}=r_{p,c3}$ , the fitted lines intersect and ${\it\alpha}_{c,max}={\it\alpha}_{cd2}$ for $r_{p}>r_{p,c3}$ .

Lundell & Carlsson (Reference Lundell and Carlsson2010) and Nilsen & Andersson (Reference Nilsen and Andersson2013) already studied the aspect ratio dependence on $\mathit{St}_{0.5}$ with negligible fluid inertia and found it to be an increasing function of $r_{p}$ . The same can be said of the problem with fluid inertia in figure 23, where $\mathit{St}_{0.5}$ is plotted as function of $\mathit{Re}_{p}$ and $r_{p}$ . The graphs are connected to the analytical value at $\mathit{Re}_{p}=0$ . From the dent in the graph at $\mathit{Re}_{p}=10$ , it can be seen that the error is larger at higher values of $r_{p}$ . This is because small numerical errors in determining $\dot{{\it\phi}}_{min}$ during the slow part of the tumbling motion for slender particles have a large influence on the error in the period, which is used for the determination of $\mathit{St}_{0.5}$ .

Figure 23. Plots of $\mathit{St}_{0.5}$ as a function of $\mathit{Re}_{p}$ for a spheroid with $r_{p}=3,4$ and 5; the black arrow shows the effect of increasing  $r_{p}$ .

8.3. Transient behaviour and orbit drift

In the region $\mathit{Re}_{\mathit{LR}}<\mathit{Re}_{p}<\mathit{Re}_{T}$ , the particle can end up in one of two final rotational states depending on the initial conditions. This splitting of solutions is described by a subcritical Hopf bifurcation and the emergence of an unstable limit cycle, outside of which initial orientations go to tumbling (associated with particle inertia), while initial orientations inside go to the coexisting solution (associated with fluid inertia). The size of the unstable limit cycle is therefore a measure of the probability of a randomly oriented particle initialized at rest ending up in the respective states. This was described in § 3.2. To determine how the unstable limit cycle changes shape due to the aspect ratio, a comparison is made of the initial orientation dependence for $r_{p}=3$ and 6 around $\mathit{Re}_{PF}$ . The simulation was initialized with the particle at rest in five different initial orientations given by (3.1), and was run up to $Gt=400$ .

Figure 24(a,b) show the projected trajectories for $r_{p}=3(\mathit{Re}_{p}=79)$ and $r_{p}=6$ $(\mathit{Re}_{p}=60)$ , respectively. It is clear that the unstable limit cycle shrinks with increasing aspect ratio, which leads to the conclusion that as $r_{p}$ increases, more initial orientations will lead to a tumbling motion. This is also supported by physical reasoning, as the inertial tensor element for minor axis rotation, $\unicode[STIX]{x1D610}_{minor}\propto r_{p}^{-2}(1+r_{p}^{-2})$ , is orders of magnitude greater than the element for major axis rotation, $\unicode[STIX]{x1D610}_{major}\propto 2r_{p}^{-4}$ , when $r_{p}\rightarrow \infty$ . It therefore seems reasonable that the tumbling state would be favoured over log-rolling as the particle becomes more slender.

Figure 24. Projected trajectories close to $\mathit{Re}_{PF}$ , for: (a)  $r_{p}=3$ and $\mathit{Re}_{p}=79$ ; (b)  $r_{p}=6$ and $\mathit{Re}_{p}=60$ . In each panel, the open circles indicate initial orientations and the filled circles indicate final orientations at $Gt=400$ ; the dashed line indicates the approximate location of an unstable limit cycle, such that initial orientations inside lead to log-rolling and initial orientations outside lead to tumbling.

8.4. Consequences for suspension rheology

In § 7, it was discussed how the intrinsic viscosity ${\it\eta}$ is modified by the rotational states at different $\mathit{Re}_{p}$ and $\mathit{St}$ . To investigate the $r_{p}$ -dependence of ${\it\eta}$ , three neutrally buoyant particles of aspect ratios $r_{p}=4,5$ and 6 are simulated and compared in figure 25. The analytical value of ${\it\eta}_{Jeffery}(\text{LR})$ at $\mathit{Re}_{p}=0$ is known to go as ${\it\eta}_{Jeffery}(\text{LR})\rightarrow 2$ when $r_{p}\rightarrow \infty$ , whereas ${\it\eta}_{Jeffery}(\text{T})\rightarrow \infty$ as $r_{p}\rightarrow \infty$ (Jeffery Reference Jeffery1922). At the same time ${\it\eta}(\mathit{Re}_{c})$ seems to decrease. This observation is also supported by the fact that ${\it\phi}_{c}(\mathit{Re}_{c})$ gets closer to zero as $r_{p}$ is increased (figure 19), which would make a smaller contribution to the suspension viscosity (Jeffery Reference Jeffery1922; Mueller, Llewellin & Mader Reference Mueller, Llewellin and Mader2009). The consequence is that the suspension goes from shear-thickening to shear-thinning behaviour on the tumbling branch as $r_{p}$ increases. This is in line with the results of Huang et al. (Reference Huang, Wu and Lu2012a ). It should be noted that, even though the shear-thickening behaviour of the non-planar branch still holds for high-aspect-ratio particles, the results from the previous section show that slender particles have a higher probability of ending up in a tumbling motion. Therefore the rheological properties of dilute fibre suspensions will rather be determined by the shear-thinning of the tumbling branch. The largest contribution to the suspension viscosity without inertia occurs when the particle is aligned with ${\it\phi}=\pm 45^{\circ }$ (Jeffery Reference Jeffery1922; Mueller et al. Reference Mueller, Llewellin and Mader2009). Knowing that the stable equilibrium angle ${\it\phi}_{c}$ is an increasing function of $\mathit{Re}_{p}$ (figure 12) and that the viscosity increases as this angle comes closer to ${\it\phi}=45^{\circ }$ , the shear-thickening behaviour on the steady-state branch when $\mathit{Re}_{p}>\mathit{Re}_{c}$ is believed to remain as the particles become more slender. However, the trend is that $\text{d}{\it\eta}/\text{d}\mathit{Re}_{p}$ decreases as $r_{p}$ increases for $\mathit{Re}_{p}>\mathit{Re}_{c}$ .

Figure 25. Intrinsic viscosity of a one-particle suspension plotted as a function of $\mathit{Re}_{p}$ with $r_{p}=4,5$ and 6 for a neutrally buoyant particle ( ${\it\alpha}=1$ ); the arrows show the effect of increasing  $r_{p}$ .

8.5. Low-aspect-ratio particles

In this section, we report our findings on the rotational behaviour of neutrally buoyant, low-aspect-ratio particles, which deviates from the behaviour described by Rosén et al. (Reference Rosén, Lundell and Aidun2014). A spheroid with $r_{p}=2$ is typically considered in this section.

Up to $\mathit{Re}_{p}=\mathit{Re}_{PF}$ , the behaviour is similar to that of the higher-aspect-ratio particles. The big difference occurs at $\mathit{Re}_{p}>\mathit{Re}_{PF}$ and is illustrated in figure 26 for $\mathit{Re}_{p}=210$ and three different initial orientations: $({\it\phi},{\it\theta})_{1}=(0,{\rm\pi}/8)$ , $({\it\phi},{\it\theta})_{2}=(0,0)$ and $({\it\phi},{\it\theta})_{3}=({\rm\pi}/2,{\rm\pi}/8)$ .

Figure 26. Projected trajectories of the $r_{p}=2$ particle’s symmetry axis $(p_{x},p_{y})=(\sin {\it\theta}\cos {\it\phi},\sin {\it\theta}\sin {\it\phi})$ at $\mathit{Re}_{p}=210$ , for different initial orientations that lead to: (a) inclined rolling; (b) ‘chaotic’ kayaking; (c) tumbling. The open circles indicate the initial orientations and the filled circles indicate the final orientations at $Gt=1000$ .

Here a third rotational state is seen to coexist with the tumbling and inclined rolling solutions. The unstable manifold of the log-rolling saddle is attracted by the other unstable manifold of the same saddle. The result is a rotational state which is seemingly chaotic. While all the other rotational states discussed previously could be attributed to either fluid or particle inertial effects, this ‘chaotic’ kayaking state seems to involve both effects. Fluid inertia dominates close to the log-rolling saddle, but the angular acceleration of the particle leads to increased particle inertial effects, bringing it close to the saddle again. The simulations were not run for long enough (simulation time $Gt=1000$ ) to exclude the possibility that the particle actually does not settle to a stable limit cycle, and thus it cannot be concluded that this state is truly chaotic.

The reason for both these effects being present can be understood by again looking at the terms in the inertial tensor of the spheroid. Using the same reasoning as in § 8.3, we find that the inertial tensor element for minor axis rotation, $\unicode[STIX]{x1D610}_{minor}\propto r_{p}^{-2}(1+r_{p}^{-2})$ , is of the same order of magnitude as the inertial tensor element for major axis rotation, $\unicode[STIX]{x1D610}_{major}\propto 2r_{p}^{-4}$ , when $r_{p}\gtrsim 1$ . Particle inertia therefore cannot be neglected when the particle is in non-planar motion and will play a role in any type of rotation.

Above $\mathit{Re}_{p}=215$ , no stable fixed point was found, and due to the existence of the ‘chaotic’ kayaking state, it is impossible to see by simple observation whether a supercritical Hopf bifurcation has occurred or if the fixed point becomes unstable because it has moved outside the unstable limit cycle. A value of $\mathit{Re}_{\mathit{Hopf}}$ thus becomes very hard to determine.

As $\mathit{Re}_{p}$ is increased further, the oscillations become so large that the particle eventually spirals out to the stable limit cycle corresponding to tumbling.

Another difference for the $r_{p}=2$ spheroid is observed when analysing the tumbling motion. Above $\mathit{Re}_{p}\approx 210$ , the particle in tumbling motion does not necessarily rotate with axis of rotation parallel to the vorticity axis. Depending on the initial conditions, it may end up in a tumbling motion for which the angle between the vorticity axis and the axis of rotation is somewhere between 0 and ${\rm\pi}/6$ . For one value of $\mathit{Re}_{p}$ , the particle can thus end up in one of an infinite number of these inclined tumbling states. The explanation might be that the particle is interacting with its own wake, which creates forces that counteract the centrifugal forces that want to take the particle to a rotation in the flow-gradient plane.

It is clear that the rotational behaviour of low-aspect-ratio particles, i.e. ones with $r_{p}\leqslant 2$ , is very complex, and further investigation is needed to make a full characterization. This is beyond the scope of the present study.