2.1. Mechanistic origins of stratification-induced force and torque

The derivation of the expressions in (2.6) is in appendix A and relies on an exact decomposition of the velocity and pressure fields into two components such that

(2.7) \begin{equation} \mathbf {u}=\mathbf {u}^{\textit{Stokes}}+\mathbf {u}^{\textit{Stratified}},\text { and } p=(1+{\beta }\eta '/\eta _0)p^{\textit{Stokes}}+p^{\textit{Stratified}}. \end{equation}

The variables $\mathbf {u}^{\textit{Stokes}}$ and $p^{\textit{Stokes}}$ represent the velocity and pressure in a uniform viscosity fluid (Stokes flow) with viscosity equal to that at the particle’s center, whereas $\mathbf {u}^{\textit{Stratified}}$ and $p^{\textit{Stokes}}{\beta }\eta '/\eta _0+p^{\textit{Stratified}}$ are velocity and pressure induced by viscosity variation. The pressure component $p^{\textit{Stokes}}{\beta }\eta '/\eta _0$ is the Stokes pressure at any position incremented by the local viscosity. The fluid stress, $\boldsymbol {\sigma }$ , consists of three parts,

(2.8) \begin{equation} \boldsymbol {\sigma }=\boldsymbol {\sigma }^{\textit{Stokes}}+({\beta }\eta '/\eta _0)\boldsymbol {\sigma }^{\textit{Stokes}}+\boldsymbol {\sigma }^{\textit{Stratified}}, \end{equation}

with its constituent stresses,

(2.9) \begin{align} \boldsymbol {\sigma }^{\textit{Stokes}}&=-p^{\textit{Stokes}}\mathbf {I}+\eta _0\big(\nabla \mathbf {u}^{\textit{Stokes}}+\big(\nabla \mathbf {u}^{\textit{Stokes}}\big)^{\rm T}\big),\text { and,} \end{align}

(2.10) \begin{align} \boldsymbol {\sigma }^{\textit{Stratified}}&=-p^{\textit{Stratified}}\mathbf {I}+(\eta _0+{\beta }\eta ')\big[\nabla \mathbf {u}^{\textit{Stratified}}+\big(\nabla \mathbf {u}^{\textit{Stratified}}\big)^{\rm T}\big]. \\[6pt] \nonumber \end{align}

After applying the decomposition, the governing equations are,

(2.11) \begin{eqnarray} \nabla \cdot \mathbf {u}^{\textit{Stokes}}=0, \hspace {0.2in}\nabla \cdot \boldsymbol {\sigma }^{\textit{Stokes}}=0, \end{eqnarray}

(2.12) \begin{eqnarray} \nabla \cdot \mathbf {u}^{\textit{Stratified}}=0, \hspace {0.2in}\nabla \cdot \boldsymbol {\sigma }^{\textit{Stratified}}+\frac {\nabla \eta }{\eta _0}.\boldsymbol {\sigma }^{\textit{Stokes}}=0, \end{eqnarray}

subject to the boundary conditions

(2.13) \begin{eqnarray} \mathbf {u}^{\textit{Stokes}}=\mathbf {u}_{\textit{particle}}, \text {on the particle surface, and, }\mathbf {u}^{\textit{Stokes}}\rightarrow \mathbf {u}_\infty \text { as }\mathbf {x}\rightarrow \mathbf {x}_{{out}},\text { and,}\,\,\,\,\,\,\,\,\,\,\,\, \end{eqnarray}

(2.14) \begin{eqnarray} \mathbf {u}^{\textit{Stratified}}=0, \text {on the particle surface, and, as }\mathbf {x}\rightarrow \mathbf {x}_{{out}}, \end{eqnarray}

where as in (2.3), the outer boundary, $\mathbf {x}_{{out}}$ , may be in the far-field, at a solid wall or a periodic boundary. Equations (2.11) and (2.13) governing the evolution of $\boldsymbol {\sigma }^{\textit{Stokes}}$ are the same as the original equations (2.1) and (2.3) but in a fluid with uniform viscosity, $\eta _0(t)$ . The velocity induced by the variable viscosity effect, $\mathbf {u}^{\textit{Stratified}}$ , has zero boundary conditions, equation (2.14), because the imposed flow and particle motion are already accounted for in the boundary conditions for Stokes velocity, $\mathbf {u}^{\textit{Stokes}}$ in (2.13). Summing (2.11) and (2.12) recovers the original formulation in (2.1). We have not made any assumption in decomposing the original system ((2.1) to (2.3)) into the two components given by (2.11) to (2.14). This exact decomposition is possible because the original mass and momentum equations are linear in velocity and pressure. Periodic boundary conditions (not in (2.13) and (2.14)) are also compatible with this decomposition.

To support the discussion and analysis of § 4, it is useful to describe the additional stratification-induced force and torque on a particle that are not present in classical Stokes flow. From (2.8) the extra stress in a variable viscosity fluid is $({\beta }\eta '(\mathbf {x})/\eta _0)\boldsymbol {\sigma }^{\textit{Stokes}}+\boldsymbol {\sigma }^{\textit{Stratified}}$ which leads to the $\mathcal {O}(\beta )$ force or torque $-2{\beta }\int _{\textit{Fluid}} {\rm d}V\hspace {0.1in}\eta '(\mathbf {e}^{\textit{Stokes}}-\mathbf {E}_\infty ):\nabla \mathbf {b}$ (with $\mathbf {b}=\mathbf {b}^{\mathbf {f}}$ or $\mathbf {b}^{\mathbf {q}}$ ). The net force and torque, shown in (2.6), acting on the particle can be represented as, $\mathbf {f}=\eta _0\mathbf {f}^{\textit{Stokes}}+\mathbf {f}^{\textit{Stratified}}_{{A}}+\mathbf {f}^{\textit{Stratified}}_{{B}}$ and $\mathbf {q}=\eta _0\mathbf {q}^{\textit{Stokes}}+\mathbf {q}^{\textit{Stratified}}_{{A}}+\mathbf {q}^{\textit{Stratified}}_{{B}}$ . Here, part of the force and torque arise from the Stokes stress acting on the particle surface immersed in varying viscosity fluid,

(2.15) \begin{align} \mathbf {f}^{\textit{Stratified}}_{{A}}&={\beta }\int _{r_p}{\rm d}S\hspace {0.05in}((\eta '/\eta _0)\mathbf {n}\cdot \boldsymbol {\sigma }^{\textit{Stokes}}), \nonumber \\ \mathbf {q}^{\textit{Stratified}}_{{A}}&={\beta }\int _{r_p}{\rm d}S\hspace {0.05in}(\eta '/\eta _0)\mathbf {x}\times (\mathbf {n}\cdot \boldsymbol {\sigma }^{\textit{Stokes}}), \end{align}

and the final contributions come from the velocity and pressure induced by variable viscosity,

(2.16) \begin{eqnarray} \begin{split} \mathbf {f}^{\textit{Stratified}}_{{B}}&=\int _{r_p}{\rm d}S\hspace {0.05in}(\mathbf {n}\cdot \boldsymbol {\sigma }^{\textit{Stratified}})\\&= -2\beta \int _{\textit{Fluid}} {\rm d}V\hspace {0.1in}\eta '(\mathbf {e}^{\textit{Stokes}}-\mathbf {E}_\infty ):\nabla \mathbf {b}^{\mathbf {f}}-\mathbf {f}^{\textit{Stratified}}_{{A}}{+\mathcal {O}(\beta ^2)}, \\ \mathbf {q}^{\textit{Stratified}}_{{B}}&=\int _{r_p}{\rm d}S\hspace {0.05in}\mathbf {x}\times (\mathbf {n}\cdot \boldsymbol {\sigma }^{\textit{Stratified}})\\&= -2\beta \int _{\textit{Fluid}} {\rm d}V\hspace {0.1in}\eta '(\mathbf {e}^{\textit{Stokes}}-\mathbf {E}_\infty ):\nabla \mathbf {b}^{\mathbf {q}}- \mathbf {f}^{\textit{Stratified}}_{{A}}{+\mathcal {O}(\beta ^2)}. \end{split} \end{eqnarray}

The force and torque, $\mathbf {f}^{\textit{Stratified}}_{{B}}$ and $\mathbf {q}^{\textit{Stratified}}_{{B}}$ , are due to the $\mathcal {O}(\beta )$ stratified stress governed by the $\mathcal {O}(\beta )$ stratified mass and momentum equations,

(2.17) \begin{equation} \nabla \cdot ({\mathbf {u}}^{\textit{Stratified}})^{(1)}=0, \hspace {0.2in}\nabla \cdot ({\boldsymbol {\sigma }}^{\textit{Stratified}})^{(1)}+\frac {\mathbf {d}}{\eta _0}\cdot \boldsymbol {\sigma }^{\textit{Stokes}}=0. \end{equation}

Here in (2.17) we have regularly expanded the flow variables in $\beta$ ,

(2.18) \begin{align} {\mathbf {u}}^{\textit{Stratified}}&=\beta ({\mathbf {u}}^{\textit{Stratified}})^{(1)}+\mathcal {O}(\beta ^2), \end{align}

(2.19) \begin{align} {p}^{\textit{Stratified}}&=\beta ({p}^{\textit{Stratified}})^{(1)}+\mathcal {O}(\beta ^2), \end{align}

(2.20) \begin{align} {\boldsymbol {\sigma }}^{\textit{Stratified}}&= \beta \left [-(\widetilde {p}^{\textit{Stratified}})^{(1)}\mathbf {I}+\eta _0[(\nabla {\mathbf {u}}^{\textit{Stratified}})^{(1)}+((\nabla {\mathbf {u}}^{\textit{Stratified}})^{(1)})^{\rm T}]\right ] +\mathcal {O}(\beta ^2). \end{align}

From here on, we focus on linear viscosity stratification such that the viscosity, $\eta$ , distribution is,

(2.21) \begin{equation} \eta =\eta _0+\beta \mathbf {d}\cdot \mathbf {x},\hspace {0.1in} ||\mathbf {d}||_2=1,\hspace {0.1in} ||\nabla \eta ||_2=\beta \hspace {0.1in} (\eta '(\mathbf {x},t)=\eta '(\mathbf {x})=\beta \mathbf {d}\cdot \mathbf {x}), \end{equation}

where unit vector $\mathbf {d}$ lies along the direction of viscosity stratification. For this case the parameter, $\beta$ in the regular perturbation expansion of the relevant flow variables is the magnitude of viscosity gradient.

Table 1.

Forces and torques (in Cartesian basis) on an axi- and fore-aft symmetric particle in body-fixed coordinates (where component 3 or the $z$ axis aligns with the particle orientation) for different flows within (3.1) for a linearly stratified fluid with viscosity given by (2.21), are $\eta _0\mathbf {f}^{\textit{Stokes}}+ \beta \mathbf {F}_{{strat}}^{\text {body}}\cdot \mathbf {d}^{\text {body}}{+\mathcal {O}(\beta ^2)}$ and $\eta _0\mathbf {q}^{\textit{Stokes}}+\beta \mathbf {Q}_{{strat}}^{\text {body}}\cdot \mathbf {d}^{\text {body}}{+\mathcal {O}(\beta ^2)}$ . The first part arises from fluid’s uniform viscosity, $\eta _0$ and the second from its viscosity gradient with magnitude $\beta =||\nabla \eta ||_2$ . The unit of $f_i$ is $l$ and that of $q_i$ , $F_{ij}^{\Gamma _k}$ and $Q_{ij}^{U_k}$ is $l^3$ , where $l$ is the chosen length scale.

3.1. Constant viscosity fluid

In a fluid with spatially constant viscosity, $\eta _0$ , a force, $\mathbf {g}$ and torque, $\mathbf {q}_{\boldsymbol {\Gamma }}$ leads to the following particle translation and rotation velocities,

(3.3) \begin{align} \mathbf {u}_{\textit{particle}}&=\frac {1}{\eta _0}\left [\mathbf {p}(\mathbf {p}\cdot {\mathbf {g}})\frac {f_1-f_3}{f_1f_3}+\frac {1}{f_1} {\mathbf {g}}\right ], \end{align}

(3.4) \begin{align} \boldsymbol {\omega }_{\textit{particle}}&=\frac {1}{\eta _0}\left [\mathbf {p}(\mathbf {p}\cdot {\mathbf {q}_{\boldsymbol {\Gamma }}})\frac {q_1-q_3}{q_1q_3}+\frac {1}{q_1} {\mathbf {q}_{\boldsymbol {\Gamma }}}\right ]. \\[6pt] \nonumber \end{align}

The significance of particle shape-dependent factors $f_1$ , $f_3$ , $q_1$ and $q_3$ can be ascertained from table 1. In particular, $f_1$ and $f_3$ respectively are the non-zero components of the hydrodynamic forces on the particle in a unit uniform flow of uniform viscosity fluid, aligned perpendicular and parallel to the axis. The factors $q_1$ and $q_3$ respectively are the non-zero components of hydrodynamic torques experienced by a fixed particle in the rotational flows $\mathbf {x}^{{body}}\cdot \boldsymbol {\Gamma }^{(7)}$ and $\mathbf {x}^{{body}}\cdot \boldsymbol {\Gamma }^{(6)}$ . A particle fixed in an imposed uniform flow with velocity, $\mathbf {u}_{{flow}}$ , experiences a force,

(3.5) \begin{equation} \mathbf {f}^{\textit{Stokes}}_{\mathbf {u}_{{flow}}}=\eta _0[\mathbf {p}(\mathbf {p}\cdot \mathbf {u}_{{flow}})(f_3-f_1)+f_1 \mathbf {u}_{{flow}}], \end{equation}

and no torque. In a linear flow with an imposed velocity gradient, $ \boldsymbol {\Gamma } = \mathbf {E}+\boldsymbol {\Omega }$ , with a symmetric (straining) part, $\mathbf {E}$ , and an anti-symmetric (rotational) part, $\boldsymbol {\Omega }$ (such that intrinsic rotation of the imposed linear flow, $\boldsymbol {\omega }_\infty =-0.5\boldsymbol {\epsilon }:\boldsymbol {\Omega }$ ), a fixed particle experiences a torque,

(3.6) \begin{equation} \mathbf {q}_{\boldsymbol {\Gamma }}=\eta _0[\mathbf {p}(\mathbf {p}\cdot \boldsymbol {\omega }_\infty )(q_3-q_1)+q_1\boldsymbol {\omega }_\infty -q_4(\mathbf {E}\cdot \mathbf {p})\times \mathbf {p}], \end{equation}

and no force. Therefore, from (3.4), a torque-free particle in a linear flow rotates with an angular velocity,

(3.7) \begin{equation} \boldsymbol {\omega }^{\textit{Jeffery}}_{\textit{particle}}=\boldsymbol {\omega }_\infty -\frac {q_4}{q_1}(\mathbf {E}\cdot \mathbf {p})\times \mathbf {p}, \end{equation}

that leads to the time rate of change of the particle orientation vector given by,

(3.8) \begin{equation} \dot {\mathbf {p}}^{\textit{Jeffery}}= \boldsymbol {\omega }^{\textit{Jeffery}}_{\textit{particle}} \times \mathbf {p} = \boldsymbol {\omega }_\infty \times \mathbf {p} + \frac {q_4}{q_1}(\mathbf {E}\cdot \mathbf {p})\cdot (\mathbf {I}-\mathbf {p}\mathbf {p}). \end{equation}

Here, $q_4$ is the non-zero component of the torque on a fixed particle in the straining flow given by $\mathbf {x}^{{body}}\cdot \boldsymbol {\Gamma }^{(4)}$ or $-\mathbf {x}^{{body}}\cdot \boldsymbol {\Gamma }^{(5)}$ . The orientation trajectories obtained through the integration of the above equation are known as Jeffery orbits (Jeffery Reference Jeffery1922).

3.2. Fixed particle in stratified fluid

In a coordinate-free form, the $\mathcal {O}(\beta )$ stratification-induced forces and torques on a fixed particle are, using values from table 1,

(3.9) \begin{equation} \mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\Gamma }}= \beta \mathbf {R}\cdot \mathbf {F}_{{strat}}^{{body}}\cdot \mathbf {R}^{{T}}\cdot \mathbf {d},\quad \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{{flow}}}= \beta \mathbf {R}\cdot \mathbf {Q}_{{strat}}^{{body}}\cdot \mathbf {R}^{\text {T}}\cdot \mathbf {d}, \end{equation}

can be expressed as,

(3.10) \begin{equation} \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{{flow}}}= \beta [{Q_{32}^{U_1}} (\mathbf {p} \mathbf {p} \cdot (\mathbf {u}_{{flow}}\times \mathbf {d}))+{Q_{12}^{U_3}} (\mathbf {p} \cdot \mathbf {u}_{{flow}})\mathbf {d}\times \mathbf {p}+ {Q_{23}^{U_1}}(\mathbf {p}\cdot \mathbf {d})\mathbf {p}\times \mathbf {u}_{{flow}}], \end{equation}

in a uniform flow with velocity, $\mathbf {u}_{{flow}}$ , and,

(3.11) \begin{align} \begin{split} &\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\Gamma }}=\beta \mathbf {R}\cdot \mathbf {F}^{{strat}}_{{linear,}}\cdot \mathbf {R}^{\text {T}}\cdot \mathbf {d} \\&=\beta \Big (F_{11}^{\Gamma _1}\mathbf {E}+(F_{33}^{\Gamma _2}+F_{11}^{\Gamma _2}+F_{11}^{\Gamma _1}/2-F^{\Gamma _4}_{13}-F^{\Gamma _4}_{31})(\mathbf {p}\cdot \mathbf {E}\cdot \mathbf {p})\mathbf {pp}\\&\quad +(F_{11}^{\Gamma _1}/2-F_{11}^{\Gamma _2})(\mathbf {p}\cdot \mathbf {E}\cdot \mathbf {p})\mathbf {I}+(F^{\Gamma _4}_{13}-F_{11}^{\Gamma _1})\mathbf {E}\cdot \mathbf {pp}+(F^{\Gamma _4}_{31}-F_{11}^{\Gamma _1})\mathbf {pp}\cdot \mathbf {E}\Big )\cdot \mathbf {d}\\&\quad -\beta \Big (F_{12}^{\Gamma _6}\boldsymbol {\omega }_{{flow}}\!\times\! \mathbf {d} -(F_{12}^{\Gamma _6}-F^{\Gamma _7}_{31})\mathbf {pp}\cdot (\boldsymbol {\omega }_{{flow}}\times \mathbf {d})-(F_{12}^{\Gamma _6}-F^{\Gamma _8}_{23})(\boldsymbol {\omega }_{{flow}}\times \mathbf {p})\mathbf {p}\cdot \mathbf {d}\Big ), \end{split} \end{align}

in a linear imposed flow with $\boldsymbol{\Omega}$ and $\mathbf {E}$ as the vorticity and strain rate tensor. Here, $\boldsymbol {\omega }_{{flow}}=-0.5\boldsymbol {\epsilon }: \boldsymbol{\Omega }$ is the imposed fluid rotation in the perspective of the particle. The physical meaning of the coefficients $Q_{ij}^{U_k}$ ( $i,j,k\in [1,3]$ ) and $F_{ij}^{\Gamma _k}$ ( $i,j\in [1,3],k\in [1,8]$ ), can be inferred from table 1. These are $\beta$ normalized stratification-induced torques, $Q_{ij}^{U_k}$ , and forces, $F_{ij}^{\Gamma _k}$ , along direction $i$ in a fluid with viscosity increasing linearly along $j$ , on a fixed particle with axis of symmetry along direction 3, either in a uniform flow (for $Q_{ij}^{U_k}$ ) with unit velocity along direction $k$ or a linear flow with gradient $\hat {\boldsymbol {\Gamma }}^{(k)}$ (for $F_{ij}^{\Gamma _k}$ ).

3.3. Particle motion in viscosity-stratified fluid

A particle translating with a velocity $\mathbf {u}_{\textit{particle}}$ in a quiescent fluid experiences a stratification-induced torque,

(3.12) \begin{align} \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{\textit{particle}}}&=-\beta [{Q_{32}^{U_1}} (\mathbf {p} \mathbf {p} \cdot (\mathbf {u}_{\textit{particle}}\times \mathbf {d}))+{Q_{12}^{U_3}} (\mathbf {p} \cdot \mathbf {u}_{\textit{particle}})\mathbf {d}\times \mathbf {p}\nonumber \\ &\quad + {Q_{23}^{U_1}}(\mathbf {p}\cdot \mathbf {d})\mathbf {p}\times \mathbf {u}_{\textit{particle}}], \end{align}

in addition to the force, $-\eta _0[\mathbf {p}(\mathbf {p}\cdot \mathbf {u}_{\textit{particle}})(f_3-f_1)+f_1 \mathbf {u}_{\textit{particle}}]$ ( $\mathbf {u}_{{flow}}$ replaced with $-\mathbf {u}_{\textit{particle}}$ in (3.5)) from the uniform viscosity component of the fluid. While a particle rotating with an angular velocity $\boldsymbol {\omega }_{\textit{particle}}$ in a quiescent fluid experiences a torque $-\eta _0[\mathbf {p}(\mathbf {p}\cdot \boldsymbol {\omega }_{\textit{particle}})(q_3-q_1)+q_1\boldsymbol {\omega }_{\textit{particle}}]$ ( $\mathbf {E}=0$ and $\boldsymbol {\omega }_\infty =-\boldsymbol {\omega }_{\textit{particle}}$ in (3.6)) due to the constant part of viscosity. The linear stratification leads to a force,

(3.13) \begin{align} \begin{split} \mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}}&=\beta \Big (F_{12}^{\Gamma _6}\boldsymbol {\omega }_{\textit{particle}}\times \mathbf {d} -(F_{12}^{\Gamma _6}-F^{\Gamma _7}_{31})\mathbf {pp}\cdot (\boldsymbol {\omega }_{\textit{particle}}\times \mathbf {d}) \\&\quad -(F_{12}^{\Gamma _6}-F^{\Gamma _8}_{23})(\boldsymbol {\omega }_{\textit{particle}}\times \mathbf {p})\mathbf {p}\cdot \mathbf {d}\Big ), \end{split} \end{align}

as a relative linear flow of the fluid with (anti-symmetric) velocity gradient $\boldsymbol {\Gamma } = $ $-\boldsymbol {\epsilon }\cdot \boldsymbol {\omega }_{\textit{particle}}$ is experienced by the particle.

The stratification-induced torque on a translating particle and force on a rotating particle lead to a coupling between two types of motion not observed in a uniform viscosity fluid in the absence of inertia (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016; Datt & Elfring Reference Datt and Elfring2019; Anand & Narsimhan Reference Anand and Narsimhan2024; Gong et al. Reference Gong, Shaik and Elfring2024). A particle sediments in a constant viscosity fluid along a linear path that depends on its shape and initial orientation. However, its orientation will change in the presence of viscosity stratification, leading to a different translation path. In a linear flow, while rotating as per Jeffery orbits, a particle translates at a velocity equal to the fluid’s velocity at its centroid. However, a stratification-induced force will lead to a relative translation between the particle and fluid. For example, unlike a uniform viscosity scenario, it can migrate across the streamlines of a simple shear flow or be forced to move relative to the center of a uniaxial extensional flow.

3.3.1. Sedimenting particle

First consider a particle settling under the action of a gravitational force, $\mathbf {g}$ in a quiescent fluid. The particle’s translation and angular velocity are obtained via the following coupled equations,

(3.14) \begin{align} \mathbf {u}_{\textit{particle}}&=\frac {1}{\eta _0}\left [\mathbf {p}(\mathbf {p}\cdot {(\mathbf {g}+\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}})})\frac {f_1-f_3}{f_1f_3}+\frac {1}{f_1} {(\mathbf {g}+\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}})}\right ], \end{align}

(3.15) \begin{align} \boldsymbol {\omega }_{\textit{particle}}&=\frac {1}{\eta _0}\left [\mathbf {p}(\mathbf {p}\cdot { \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{\textit{particle}}} })\frac {q_1-q_3}{q_1q_3}+\frac {1}{q_1} { \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{\textit{particle}}} }\right ], \\[6pt] \nonumber \end{align}

where $\mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{\textit{particle}}}$ (a function of $\mathbf {u}_{\textit{particle}}$ ) and $\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}}$ (a function of $\boldsymbol {\omega }_{\textit{particle}}$ ) are given by (3.12) and (3.13). Upon substitution of relevant variables, the governing equation for angular velocity is,

(3.16) \begin{eqnarray} &\boldsymbol {\omega }_{\textit{particle}} = \frac {\beta }{\eta _0^2}[-t_1 (\mathbf {p} \mathbf {p} \cdot ({\mathbf {g}}\times \mathbf {d}))-t_2 (\mathbf {p} \cdot {\mathbf {g}})\mathbf {d}\times \mathbf {p}-t_3(\mathbf {p}\cdot \mathbf {d})\mathbf {p}\times {\mathbf {g}}]+\mathcal {O}(\beta ^2), \end{eqnarray}

where,

(3.17) \begin{eqnarray} t_1=\frac {Q_{32}^{U_1}}{f_1 q_3}, \quad t_2=\frac {Q_{12}^{U_3}}{f_3 q_1} ,\quad t_3=\frac {Q_{23}^{U_1}}{f_1 q_1}. \end{eqnarray}

The $\mathcal {O}(\beta )$ effect of stratification is to rotate the particle’s centerline, $\mathbf {p}$ . We note that $\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}}$ is $\mathcal {O}(\beta ^2)$ , since, from (3.13), $\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}}\sim \beta \mathcal {O}(\boldsymbol {\omega }_{\textit{particle}})$ and $\boldsymbol {\omega }_{\textit{particle}}$ is itself $\mathcal {O}(\beta )$ , which can be discerned from (3.13) and (3.15). Therefore, up to $\mathcal {O}(\beta )$ , the settling velocity of the particle is simply the Newtonian velocity in (3.3) for the $\mathbf {p}$ altered by the stratification induced torque. At $\mathcal {O}(\beta ^2)$ , a modification in the particle’s translational velocity due to $\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\omega }_{\textit{particle}}}$ is expected. The time rate of change of the particle orientation vector, $\dot {\mathbf {p}}=\boldsymbol {\omega } \times \mathbf {p}$ , is governed by,

(3.18) \begin{align} \begin{split} \dot {\mathbf {p}}=-\frac {\beta }{\eta _0^2}\cdot \left [ t_3(\mathbf {p}\cdot \mathbf {d}){\mathbf {g}}-t_2 (\mathbf {p} \cdot {\mathbf {g}})\mathbf {d}\right ]\cdot (\mathbf {I}-\mathbf {pp})+\mathcal {O}(\beta ^2). \end{split} \end{align}

3.3.2. Freely moving particles in linear flows

The coupling between the rotation and translation of a freely suspended particle in a linear flow with velocity gradient, $\boldsymbol {\Gamma }=\mathbf {E}+\boldsymbol {\Omega }$ leads to the following translation and angular velocity of the particle,

(3.19) \begin{align} &\mathbf {u}_{\textit{particle}}=\mathbf {x}_p\cdot \boldsymbol {\Gamma }+\mathbf {u}^{\textit{Stratified}},\quad \mathbf {u}^{\textit{Stratified}}=\frac {1}{\eta _0}\left [\mathbf {p}(\mathbf {p}\cdot {\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\Gamma }}})\frac {f_1-f_3}{f_1f_3}+\frac {1}{f_1} {\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\Gamma }}}\right ], \end{align}

(3.20) \begin{align} &\boldsymbol {\omega }_{\textit{particle}}=\boldsymbol {\omega }_\infty -\frac {q_4}{q_1}(\mathbf {E}\cdot \mathbf {p})\times \mathbf {p}+ \frac {1}{\eta _0}\left [\mathbf {p}(\mathbf {p}\cdot { ( \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{{flow-particle}}})})\frac {q_1-q_3}{q_1q_3}+\frac {1}{q_1} {( \mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{\text {flow-particle}}}) }\right ]. \end{align}

where, $\mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\Gamma }}$ (given in (3.11) with $\mathbf {E}$ as the symmetric part of imposed velocity gradient, and $\boldsymbol {\omega }_{{flow}}=\boldsymbol {\omega }_\infty -\boldsymbol {\omega }_{\textit{particle}}$ ) is the stratification-induced force on a fixed particle in a linear flow created by the imposed strain and relative rotation of the fluid in the perspective of the particle. The expressions for stratification-induced torques due to relative translation motion between particle and fluid, $\mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{{flow-particle}}}$ is $\mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{\textit{particle}}}$ from (3.12) with $\mathbf {u}_{\textit{particle}}$ replaced with $\mathbf {u}^{\textit{Stratified}}$ . This relative velocity of the particle induced by stratification is governed by,

(3.21) \begin{align} \begin{split} &\mathbf {u}^{\textit{Stratified}}=\frac {\beta }{f_1\eta _0}\left (F_{11}^{\Gamma _1}\mathbf {E}+(F_{11}^{\Gamma _1}/2+F_{11}^{\Gamma _2})(\mathbf {p}\cdot \mathbf {E}\cdot \mathbf {p})\mathbf {I} \right. \\& + \left [\frac {f_1}{f_3}F_{33}^{\Gamma _2}-F_{11}^{\Gamma _2}+F_{11}^{\Gamma _1}/2-F^{\Gamma _4}_{13}-\frac {f_1}{f_3}F^{\Gamma _4}_{31}+\frac {q_4}{q_1}\left (\frac {f_1}{f_3}F_{31}^{\Gamma _7}-F_{23}^{\Gamma _8}\right )\right ](\mathbf {p}\cdot \mathbf {E}\cdot \mathbf {p})\mathbf {pp}\\& \left . +\left (F^{\Gamma _4}_{13}-F_{11}^{\Gamma _1}+\frac {q_4}{q_1}F^{\Gamma _8}_{23}\right )\mathbf {E}\cdot \mathbf {pp}+\left (\frac {f_1}{f_3}F^{\Gamma _4}_{31}-F_{11}^{\Gamma _1}-\frac {q_4}{q_1}\frac {f_1}{f_3}F^{\Gamma _7}_{31}\right )\mathbf {pp}\cdot \mathbf {E}\right )\cdot \mathbf {d}+\mathcal {O}(\beta ^2). \end{split} \end{align}

Since $\mathbf {q}^{\textit{Stratified}}_{\mathbf {u}_{{flow-particle}}}\sim \beta \mathcal {O}(\mathbf {u}^{\textit{Stratified}})\sim \mathcal {O}(\beta ^2)$ , up to $\mathcal {O}(\beta )$ the particle’s centerline follows the Newtonian rotation rate (3.8).

In the next section, the force and torque components listed in table 1 are analytically obtained for spheroids using spheroidal harmonics. These coefficients appear in the particle dynamics equations (3.18) and (3.21). Free settling in gravity using equations (3.3) and (3.18) is explored in § 5 and the motion of particles freely suspended in linear flows calculated using (3.8) and (3.21) is shown in § 6.

4.1. Uniform flow

As mentioned in § 3.2, $Q_{ij}^{U_k}$ (for $i,j,k\in [1,3]$ ) refers to the $\beta$ normalized stratification-induced torque along direction $i$ in a fluid with viscosity increasing linearly along $j$ , on a fixed particle with axis of symmetry along direction 3, in a uniform flow with unit velocity along direction $k$ . Out of the 27 coefficients in $Q_{ij}^{U_k}$ there are only three unique non-zero values. Therefore, the stratification-induced torque on a spheroid fixed in uniform flow of linearly stratified fluid for a general case can be described by the torque induced on the particle in three distinct scenarios which all require the imposed flow and the stratification direction to be mutually perpendicular. In the first case, the particle centerline, denoted by unit vector $\mathbf {p}$ , lies normal to the flow-stratification plane and the stratification-induced torque, denoted by $Q_{32}^{U_1}$ , lies along $\mathbf {p}$ . In the remaining two scenarios, where either the flow direction (torque labeled as $Q_{12}^{U_3}$ ) or the stratification direction (torque denoted by $Q_{23}^{U_1}$ ) is along $\mathbf {p}$ , the stratification-induced torque is normal to $\mathbf {p}$ . In § 4.1.1 we provide a validation of the expressions for $Q_{32}^{U_1}$ , $Q_{23}^{U_1}$ and $Q_{12}^{U_3}$ before discussing the qualitative trends in the variation of these quantities with particle aspect ratio, $\kappa$ , and providing a mechanistic explanation of their origins in § 4.1.2.

Figure 1.

Viscosity-stratification-induced torque in uniform flow past a spheroid with aspect ratio $\kappa$ and major axis equal to 1. The markers for $1\le \kappa \le 30$ are the values obtained from the numerical code.

4.1.1. Validation

The analytical expressions for stratification-induced torques $Q_{32}^{U_1}$ , $Q_{23}^{U_1}$ and $Q_{12}^{U_3}$ for a spheroid fixed in uniform flow ( $\hat {\mathbf {u}}^{(i)}, i\in [1,3]$ from equation (3.1)) are shown in (D4) within appendix D.2 and are plotted in figure 1 as solid curves. The corresponding numerically obtained values for a prolate spheroid with aspect ratio, $1\le \kappa \le 30$ are shown as circular markers in the same figure where a close match between the values obtained from the two different techniques can be observed. Our formulas show that in the limit of a sphere of radius, $l$ , i.e. $\kappa \rightarrow 1$ ,

(4.1) \begin{equation} \lim _{\kappa \rightarrow 1}(Q_{23}^{U_1})=\lim _{\kappa \rightarrow 1}(Q_{12}^{U_3})=-\lim _{\kappa \rightarrow 1}(Q_{32}^{U_1})=2\pi l^3, \end{equation}

which are the same results obtained for a sphere in appendix C (equation (C5)) without using the spheroidal harmonics.

For a slender fiber, i.e., a prolate spheroid with $\kappa \rightarrow \infty$ ,

(4.2) \begin{eqnarray} \begin{split} &\lim _{\kappa \rightarrow \infty }(Q_{32}^{U_1})=-\frac {32\pi }{3}l^3\frac {\frac {\kappa }{\sqrt {\kappa ^2-1}}-1}{2\log (2\kappa )+1}+\mathcal {O}(\kappa ^{-2}),\\ &\lim _{\kappa \rightarrow \infty }(Q_{23}^{U_1})=\frac {8\pi }{3}l^3\frac {1}{\log (2\kappa )-\frac {1}{4\log (2\kappa )}}+\mathcal {O}\Big (\frac {\kappa ^{-2}}{\log (\kappa )}\Big ), \\ &\lim _{\kappa \rightarrow \infty }(Q_{12}^{U_3})=-\frac {8\pi }{3}l^3\frac {1}{(2\log (2\kappa )-1)^2}+\mathcal {O}\Big (\frac {\kappa ^{-2}}{\log (\kappa )}\Big ).\end{split} \end{eqnarray}

Kamal & Lauga (Reference Kamal and Lauga2023) evaluated the torque equivalent to $Q_{23}^{U_1}$ for a slender fiber using resistive force theory. According to their equation (3.21), a slender prolate spheroid of length $2l$ fixed in a fluid with viscosity gradient $\eta _0/(2l)$ parallel to a fiber undergoing a uniform flow with velocity $U_2$ perpendicular to the fiber experiences a torque $4\pi /[3(\log (2\kappa )+1)]l^2 \eta _0 U_2+\mathcal {O}(\log (\kappa )^{-2})$ along the axis normal to the flow-stratification plane. Their torque value re-scaled for $U_2=1$ and viscosity gradient $\beta$ is $(\beta 8\pi /3)l^3/ (\log (2\kappa )+1)+\mathcal {O}(\log (\kappa )^{-2})$ . The equivalent torque (equation (4.2)), from our study in a fluid is $\beta \lim \limits _{\kappa \rightarrow \infty }Q_{23}^{U_1}=(8\pi /3) l^3\beta (1/[\log (2\kappa )-(1/4(\log (2\kappa )))])+\mathcal {O}(\kappa ^{-2}\log (\kappa )^{-1})$ . Hence, in the large $\kappa$ limit our expression matches with that of Kamal & Lauga (Reference Kamal and Lauga2023).

4.1.2. Mechanistic origin of stratification-induced torque

The decomposition of torque provided in (2.15) and (2.16) shows that the stratification-induced torque arises from two distinct fluid stresses (see (2.8)) in a viscosity stratified fluid flow around a particle. The $\beta$ normalized stress component $(\eta '/\eta _0)\boldsymbol {\sigma }^{\textit{Stokes}}$ depicts the stress on a particle in a uniform viscosity fluid, $\boldsymbol {\sigma }^{\textit{Stokes}}$ , but incremented by the local variation in fluid viscosity around the particle. It leads to the stratification-induced torques ${Q_{32}^{U_1}}_{{A}}$ , ${Q_{23}^{U_3}}_{{A}}$ and ${Q_{12}^{U_3}}_{{A}}$ . The remaining part of the stratification-induced torques, ${Q_{32}^{U_1}}_{{B}}$ , ${Q_{23}^{U_3}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ arise due to the fluid stress resulting from the changes in stratification-induced velocity ( $\mathbf {u}^{\textit{Stratified}}$ ) and pressure ( $p^{\textit{Stratified}}$ ) at $\mathcal {O}(\beta )$ .

First consider the case when the imposed uniform flow and viscosity stratification are perpendicular to one another and also to the particle centerline $\mathbf {p}$ . Then the stratification-induced torque, $Q_{32}^{U_1}$ , acts along $\mathbf {p}$ . We find,

(4.3) \begin{equation} {Q_{32}^{U_1}}_{{B}}=0 \rightarrow Q_{32}^{U_1}={Q_{32}^{U_1}}_{{A}}, \end{equation}

for all $\kappa$ . Therefore, $Q_{32}^{U_1}$ is only due the Stokes stress acting in a variable viscosity environment. The fluid stress, $\boldsymbol {\sigma }^{\textit{Stratified}}$ , created by the stratification-induced velocity and pressure does not contribute to this torque. However, both $(\eta '/\eta _0)\boldsymbol {\sigma }^{\textit{Stokes}}$ and $\boldsymbol {\sigma }^{\textit{Stratified}}$ contribute towards $Q_{12}^{U_3}$ and $Q_{23}^{U_1}$ . These are the stratification-induced torques induced normal to $\mathbf {p}$ when either the stratification or the flow direction is along $\mathbf {p}$ . The expressions for decomposed torques ${Q_{23}^{U_1}}_{{A}}$ , ${Q_{12}^{U_3}}_{{A}}$ , ${Q_{23}^{U_1}}_{{B}}$ , and ${Q_{12}^{U_3}}_{{B}}$ are shown in (D5) and plotted in figure 2(a) as a function of particle aspect ratio, $\kappa$ . The torques ${Q_{23}^{U_1}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ are zero for $\kappa =1$ , i.e., the stratification-induced velocity and pressure do not play a role in the stratification-induced torque on a sphere. However, for $\kappa \ne 1$ , ${Q_{23}^{U_1}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ are present in addition to ${Q_{23}^{U_1}}_{{A}}$ and ${Q_{12}^{U_3}}_{{A}}$ .

Figure 2.

Decomposition of stratification-induced torques in uniform flow on spheroids vs. aspect ratio: (a) Decomposition into $\mathbf {q}^{\textit{Stratified}}_{{A}}$ and $\mathbf {q}^{\textit{Stratified}}_{{B}}=\mathbf {q}^{\textit{Stratified}}-\mathbf {q}^{\textit{Stratified}}_{{A}}$ given by (D5) (the markers show values from the numerical code), and (b) Decomposition of $\mathbf {q}^{\textit{Stratified}}_{{B}}$ into that from stratification-induced pressure and strain obtained from the numerical code.

From figure 1, we can observe that the net stratification-induced torques, $Q_{23}^{U_1}$ and $Q_{12}^{U_3}$ change signs at $\kappa \approx 0.55$ and $\kappa \approx 2.0$ respectively. As seen in figure 2(a), the change in signs at $\kappa \approx 0.55$ and $\kappa \approx 2.0$ is due to the competition between the torque arising from $(\eta '/\eta _0)\boldsymbol {\sigma }^{\textit{Stokes}}$ ( ${Q_{23}^{U_1}}_{{A}}$ and ${Q_{12}^{U_3}}_{{A}}$ ) and that from $\boldsymbol {\sigma }^{\textit{Stratified}}$ ( ${Q_{23}^{U_1}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ ). In the upcoming discussion of § 5, this sign change will play an important role in the particle’s rotational dynamics.

The torques ${Q_{23}^{U_1}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ can be further decomposed into those arising from the stratification-induced pressure, ${p}^{\textit{Stratified}}$ , and those from viscous force per unit area $2(\eta _0)\boldsymbol {e}^{\textit{Stratified}}$ , where $\boldsymbol {e}^{\textit{Stratified}}=\nabla \mathbf {u}^{\textit{Stratified}}+(\nabla \mathbf {u}^{\textit{Stratified}})^{\rm T}$ at $\mathcal {O}(\beta )$ . This decomposition can be accessed for the prolate spheroid through the numerical calculation and is shown in figure 2(b) for ${Q_{23}^{U_1}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ , marked with subscripts ‘Pressure’ and ‘Strain.’ Here we observe that the torque due to the stratification-induced pressure and viscous force per unit area, act in the same direction, but the former dominates for prolate spheroids with $\kappa \lessapprox 10$ , beyond which both are of almost equal magnitude. Furthermore, $Q_{12,B,{Pressure}}^{U_3}$ and $Q_{23,B,{Pressure}}^{U_1}$ (as well as $Q_{12,A,{Strain}}^{U_3}$ and $Q_{23,B,{Strain}}^{U_1}$ ) are equal and opposite.

Similarly, the torques ${Q_{23}^{U_1}}_{{A}}$ and ${Q_{12}^{U_3}}_{{A}}$ that arise from the Stokes stress acting in a variable-viscosity environment can be decomposed into that arising from $(\eta '/\eta _0)p^{\textit{Stokes}}$ and $2\eta '\mathbf {e}^{\textit{Stokes}}$ (not shown). Since the pressure acts normal to the surface, and for a sphere the line of action of the pressure force acts through the particle surface, pressure does not contribute to the torque. Therefore, for a sphere the entire stratification-induced torque arises from $2\eta '\mathbf {e}^{\textit{Stokes}}$ , the Stokes strain rate acting in a variable-viscosity environment.

The stratification-induced pressure distribution, $p^{\textit{Stratified}}$ , is the dominant contributor to the torques ${Q_{23}^{U_1}}_{{B}}$ and ${Q_{12}^{U_3}}_{{B}}$ for prolate spheroids with $\kappa \lessapprox 10$ . Using patterns of $\eta 'p^{\textit{Stokes}}$ and $p^{\textit{Stratified}}$ at the surface of prolate spheroids we can understand the mechanisms in various torque distributions observed in figure 2.

Figure 3.

(a) to (c) Stokes pressure $p^{\textit{Stokes}}$ in uniform flow on a particle surface leading to $Q_{23, A}^{U_1}$ via $\eta '/\eta _0p^{\textit{Stokes}}$ on (left) $\kappa =1$ , (middle) $\kappa =2$ and (right) $\kappa =8$ prolate spheroid. (d) to (e) Stratification-induced pressure $p^{\textit{Stratified}}/\beta$ leading to $Q_{23, B}^{U_1}$ for the same $\kappa$ as (a) to (c) respectively. Flow ( $\mathbf {u}_\infty$ ) and stratification ( $\nabla \eta$ ) are respectively perpendicular and parallel to the particle axis of symmetry. Background contours show viscosity variation, $\eta '$ with light (dark) green representing $\eta '\gt 0$ ( $\eta '\lt 0$ ). Streamlines are of the Stokes flow, $\mathbf {u}^{\textit{Stokes}}$ . Red (blue) indicates a positive (negative) surface pressure ( $p^{\textit{Stokes}}$ or $p^{\textit{Stratified}}$ ).

First consider the pressure distribution responsible for torque $Q_{23,B,{Pressure}}^{U_1}$ (and $Q_{23,A,{Pressure}}^{U_1}$ ), i.e., for the scenario when imposed flow and stratification are perpendicular and parallel to $\mathbf {p}$ respectively. Figure 3 shows the pressures $p^{\textit{Stokes}}$ and $p^{\textit{Stratified}}$ on $\kappa =1$ , 2 and 8 prolate spheroids along with the background viscosity variation, $\eta '$ . Here, the variable component of viscosity, $\eta '=0$ at the center of the particle. For every $\kappa$ , the Stokes pressure, $p^{\textit{Stokes}}$ is positive on the particle surface facing the flow. On this flow-facing/ upstream side of the particle, the stratification-induced pressure, $p^{\textit{Stratified}}$ is positive when $\eta '\gt 0$ (second quadrant) and negative when $\eta '\lt 0$ (third quadrant). On the downstream side, $p^{\textit{Stokes}}$ , is negative and the sign of $p^{\textit{Stratified}}$ is also reversed from the upstream side. Overall, the $p^{\textit{Stratified}}$ distribution remains qualitatively similar as the particle aspect ratio is altered (figure 3). The force due to pressure acts along the particle’s center for a sphere leading to a zero moment arm for each surface element or a torque per unit area of zero, i.e. $\mathbf {x}\times \mathbf {n}=0$ on a sphere since $\mathbf {x}=l\mathbf {n}$ for the sphere with radius, $l$ . As $\kappa$ is changed from 1, the surface normal no longer points towards the particle center and there is a non-zero torque per unit area due to pressure. Upon increasing $\kappa$ from 1, the magnitude of $\mathbf {x}\times \mathbf {n}$ near the particle ends increases and the peak locations of $|p^{\textit{Stratified}}|$ move towards the particle ends (figures 3 e, 3 f and 3g) leading to an initial increase in $Q_{23,B,{Pressure}}^{U_1}$ with $\kappa$ in figure 2(b). However, the decreasing surface area and reduction in $|p^{\textit{Stratified}}|$ , upon increasing $\kappa$ for a fixed major axis length, causes the final decrease of $Q_{23,B,{Pressure}}^{U_1}$ with $\kappa$ . Therefore, a maximum in $Q_{23,B,{Pressure}}^{U_1}$ is observed at $\kappa \lessapprox 2$ in figure 2(b). The pressure distribution is such that the torque $Q_{23,B,{Pressure}}^{U_1}$ is along $\nabla \eta \times \mathbf {u}_\infty$ or clockwise in the view shown in figures 3 and is depicted as positive for all $\kappa \gt 1$ in figure 2(b).

Figure 4.

Same legend as figure 3 but for torques $Q_{12, A}^{U_3}$ and $Q_{12, B}^{U_3}$ , i.e., flow ( $\mathbf {u}_\infty$ ) and stratification ( $\nabla \eta$ ) are respectively parallel and perpendicular to the particle axis of symmetry.

While figure 3 depicted the case where the viscosity gradient direction $\mathbf {d}$ and the particle’s centerline $\mathbf {p}$ are aligned, and both perpendicular to $\mathbf {u}_\infty$ , figure 4 shows the case where $\mathbf {u}_\infty$ and $\mathbf {p}$ are aligned, and both are perpendicular to $\mathbf {d}$ . In other words the distinction between figures 4 and 3 is that the viscosity stratification, $\nabla \eta$ , and flow, $\mathbf {u}_\infty$ , directions are swapped. The distribution of $p^{\textit{Stratified}}$ shown in the bottom panels of figure 4 that leads to $Q_{12,B,{Pressure}}^{U_3}$ shown in figure 2(b) corresponds to the scenario when flow and stratification are qualitatively similar to the $p^{\textit{Stratified}}$ distributions shown in the bottom panel of figure 3 that are discussed above. Therefore, while $Q_{23,B,{Pressure}}^{U_1}$ along the direction $\nabla \eta \times \mathbf {u}_\infty$ , $Q_{12,B,{Pressure}}^{U_3}$ is along $\mathbf {u}_\infty \times \nabla \eta$ , leading to an opposite sign of $Q_{23,B,{Pressure}}^{U_1}$ and $Q_{12,B,{Pressure}}^{U_3}$ for each $\kappa$ shown in figure 2(b). The pressure distributions $p^{\textit{Stokes}}$ shown in the top panels of figures 3 and 4 are such that the torque due to $(\eta '/\eta _0)p^{\textit{Stokes}}$ is positive in both cases, which leads to a positive contribution to $Q_{23,A}^{U_1}$ and $Q_{12,A}^{U_3}$ in figure 2(a).

4.2. Linear flows

As stated in § 3.2, $F_{ij}^{\Gamma _k}$ (for $i,j\in [1,3],k\in [1,8]$ ) denotes the stratification-induced force along direction $i$ in a fluid with viscosity increasing linearly along direction $j$ with $\beta =1$ , on a particle with axis of symmetry along direction 3, in linear flows with gradient $\hat {\boldsymbol {\Gamma }}^{k}$ . The flows corresponding to various $\hat {\boldsymbol {\Gamma }}^{k}, k\in [1,8]$ labeled in (3.1) are defined relative to the particle centerline directed along a unit vector $\mathbf {p}$ . These are: planar extensional flow normal to $\mathbf {p}$ ( $k=1$ ), uniaxial extensional flow with extensional axis along $\mathbf {p}$ ( $k=2$ ), planar straining flow in the plane perpendicular ( $k=3$ ) and parallel ( $k=4$ and 5) to $\mathbf {p}$ and purely rotational flows with vorticity directed along ( $k=6$ ) and perpendicular ( $k=7$ and 8) to $\mathbf {p}$ . §§ 4.2.1 and 4.2.2 below discuss the validation and mechanistic origin of $F_{ij}^{U_k}$ (for $i,j\in [1,3],k\in [1,8]$ ).

Figure 5.

Variation of viscosity-stratification-induced forces on a spheroid aspect ratio, $\kappa$ , and major axis of 1 in a linear flow. The markers for $1\le \kappa \le 30$ are the values obtained from the numerical code.

4.2.1. Validation

The analytical expressions for the non-zero $F_{ij}^{\Gamma _k},i,j\in [1,3],k\in [1,8]$ are shown in appendix D.3 (D10). Solid curves in figure 5 show the $\kappa$ variation of these forces along with black symbols obtained from the finite difference based numerical solution of (2.12). Similar to the torques in the previous section (figure 1) a close match between the symbols and solid curves of figure 5 obtained from two different techniques serves as a point of validation. Using resistive force theory for a slender fiber, Kamal & Lauga (Reference Kamal and Lauga2023) obtain the force equivalent to $F_{33}^{\Gamma _2}$ and $F_{23}^{\Gamma _8}$ (their equations 3.28 and 3.32) as $4\pi l^3 /(3\log (2\kappa )-2)$ and $-8\pi l^3 /(3\log (2\kappa )-1)$ with an error of $\mathcal {O}(\log (\kappa )^{-2})$ . From our expressions in (D10), in the limit of a slender prolate spheroid, i.e., $\kappa \rightarrow \infty$ with a major radius, $l$ ,

(4.4) \begin{eqnarray} \begin{split} &\lim \limits _{\kappa \rightarrow \infty }(F_{11}^{\Gamma _2})=\lim \limits _{\kappa \rightarrow \infty }(F_{31}^{\Gamma _4})=\mathcal {O}((\log (\kappa ))^{-2}),\lim \limits _{\kappa \rightarrow \infty }(F_{11}^{\Gamma _1})=\lim \limits _{\kappa \rightarrow \infty }(F_{31}^{\Gamma _7})=\mathcal {O}(\kappa ^{-2}),\\&\lim \limits _{\kappa \rightarrow \infty }(F_{12}^{\Gamma _6})=\mathcal {O}(\kappa ^{-2}(\log (\kappa ))^{-1}), \lim \limits _{\kappa \rightarrow \infty }(F_{23}^{\Gamma _8})=-\frac {8\pi l^3}{3}\frac {1}{\log (2\kappa )-\frac {0.25}{\log (2\kappa )}}+\mathcal {O}(\kappa ^{-2}),\\&\lim \limits _{\kappa \rightarrow \infty }(F_{13}^{\Gamma _4})=\frac {32\pi l^3}{3}\frac {\log (2\kappa )}{\log (2\kappa ^2)^2+\log (4)\log (2\kappa ^2)+\log (2)^2-1}+\mathcal {O}(\kappa ^{-2}),\\ &\lim \limits _{\kappa \rightarrow \infty }(F_{33}^{\Gamma _2})=\frac {4\pi l^3}{3}\frac {1}{\log (2\kappa )\frac {\log (2\kappa )-2}{\log (2\kappa )-1}+\frac {0.75}{\log (2\kappa )-1}}+\mathcal {O}(\kappa ^{-2}). \end{split} \end{eqnarray}

Hence, values from our expressions in the large $\kappa$ limit match those of Kamal & Lauga (Reference Kamal and Lauga2023). In the limit of the sphere of radius $l$ , $\kappa \rightarrow 1$ ,

(4.5) \begin{align} \begin{split} &\lim \limits _{\kappa \rightarrow 1}(F_{11}^{\Gamma _1})=-2\lim \limits _{\kappa \rightarrow 1}(F_{11}^{\Gamma _2})=\lim \limits _{\kappa \rightarrow 1}(F_{33}^{\Gamma _2})=\lim \limits _{\kappa \rightarrow 1}(F_{13}^{\Gamma _4})=\lim \limits _{\kappa \rightarrow 1}(F_{31}^{\Gamma _4})=6\pi l^3,\\ &\lim \limits _{\kappa \rightarrow 1}(F_{12}^{\Gamma _6})=\lim \limits _{\kappa \rightarrow 1}(F_{31}^{\Gamma _7})=\lim \limits _{\kappa \rightarrow 1}(F_{23}^{\Gamma _8})=-2\pi l^3. \end{split} \end{align}

Substituting these formulae into equation (3.11), we find the force on a fixed sphere in linear flow, with gradient $\boldsymbol {\Gamma }=\mathbf {E}+\boldsymbol {\Omega }$ is,

(4.6) \begin{equation} \mathbf {f}^{\textit{Stratified}}_{\boldsymbol {\Gamma }}=\mathbf {f}_{{sphere}}=2\pi \beta l^3(3\mathbf {E}-\boldsymbol {\Omega })\cdot \mathbf {d}. \end{equation}

This is the same expressions as obtained without spheroidal harmonics but directly from the calculation using the flow fields around a sphere in appendix C ((C6) and (C6)). This provides another validation of our use of spheroidal harmonics.

4.2.2. Mechanistic origin of stratification-induced force

For brevity, we only discuss the mechanisms on a sphere in this section. The signs of most of the stratification-induced-torques on a $\kappa \ne 1$ spheroid are the same as those for a sphere (figure 5) and can be explained in a qualitatively similar manner as those for a sphere considered here.

Following the decomposition introduced in (2.15) and (2.16), the stratification-induced force, $\mathbf {f}_{{sphere}}$ from equation (4.6) can be split into two parts. The first arises from the Stokes stress acting in a variable-viscosity environment $(\eta '/\eta _0)\boldsymbol {\sigma }^{\textit{Stokes}}$ , and its force is, $\mathbf {f}_{{A,sphere}}=4\pi \beta l^3(\mathbf {E}+\boldsymbol {\Omega })\cdot \mathbf {d}$ . The second is due to the stratification-induced stress, $\boldsymbol {\sigma }^{\textit{Stratified}}$ , and its force is $\mathbf {f}_{{B,sphere}}=\mathbf {f}_{{sphere}}-\mathbf {f}_{{A,sphere}}$ . Within $\mathbf {f}_{{A,sphere}}$ , the pressure, $(\eta '/\eta _0)p^{\textit{Stokes}}$ contributes an amount $\mathbf {f}_{{A,sphere}}^{{Pressure}}=(8\pi /3)\beta l^3\mathbf {E}\cdot \mathbf {d}$ and the remaining $(4\pi /3)\beta l^3(\mathbf {E}+3\boldsymbol {\Omega })\cdot \mathbf {d}$ arises from the viscous stress, $2\eta '\mathbf {e}^{\textit{Stokes}}$ . A freely suspended sphere’s rotation effectively negates the rotation part ( $\boldsymbol {\Omega }$ ) of the imposed linear flow and the sign of the force arising from the straining part ( $\mathbf {E}$ ) is same from all the decomposed components discussed above ( $\boldsymbol {\sigma }^{\textit{Stratified}}$ , $2\eta '\mathbf {e}^{\textit{Stokes}}$ and $(\eta '/\eta _0)p^{\textit{Stokes}}$ ). Therefore, we use the contours of $p^{\textit{Stokes}}$ along with the background viscosity variation to shed light on the origins of the force arising from $(\eta '/\eta _0)p^{\textit{Stokes}}$ , i.e., $\mathbf {f}_{{A,sphere}}^{{Pressure}}=(8\pi /3)\beta l^3\mathbf {E}\cdot \mathbf {d}$ in uniaxial extension and simple shear flows. The other components of the stratification-induced force in straining flows follow a similar mechanism.

Consider a simple shear flow with strain rate $E_{ij}=0.5(\delta _{i2}\delta _{j1}+\delta _{j2}\delta _{i1})$ such that 1 is the flow, 2 the velocity gradient and 3 the vorticity direction such that, $\mathbf {E}\cdot \mathbf {d}=0.5\begin{bmatrix} d_2&d_1&0 \end{bmatrix}^{\rm T}$ . A sphere in simple shear flow of uniform viscosity fluid is force-free. However, if the fluid’s viscosity increases along the velocity gradient direction ( $d_1=0$ , $d_2=1$ and $d_3=0$ ), the particle will experience a force along the flow direction (perpendicular to the viscosity gradient). Alternatively, if the viscosity increment is along the flow direction ( $d_1=1$ , $d_2=d_3=0$ ), the stratification-induced force is in the velocity gradient direction (again perpendicular to viscosity gradient). The contribution of this force coming from $\eta 'p^{\textit{Stokes}}$ for a horizontal (flow direction) viscosity stratification can be understood through figure 6(a). Relative to the mean, $p^{\textit{Stokes}}$ does not create a force in a constant viscosity fluid as this pressure is both left-right and top-down anti-symmetric. The pressure $(\eta '/\eta _0)p^{\textit{Stokes}}$ is also top-down anti-symmetric, but left-right symmetric. This leads to a net force upwards, i.e., across the streamlines of imposed flow. Similarly if viscosity increases upwards instead ( $d_1=0$ , $d_2=1$ , $d_3=0$ ), $\eta 'p^{\textit{Stokes}}$ , is left right anti-symmetric, but top down symmetric creating a force towards the flow direction. The pressure distribution is the same on the back half of the sphere (not shown). Therefore, if viscosity stratification is entirely along the vorticity direction (perpendicular to the plane of the picture) of the imposed simple shear flow the stratification-induced force is zero.

Stratification also leads to a force on an otherwise force-free sphere in uniaxial extensional flow (strain rate $E_{ij}=\delta _{i1}\delta _{j1}-0.5(\delta _{i2}\delta _{j2}+\delta _{i3}\delta _{j3})$ ). If the viscosity increases along the extensional axis this force is towards the higher viscosity region. This can be explained through the $p^{\textit{Stokes}}$ contours on the particle surface along with background $\eta '$ shown in figure 6(b). The distribution of $p^{\textit{Stokes}}$ on the particle surface is such that it is negative near the extensional axis and positive near the compression plane. Despite local compressive and extensional forces the net force is zero, so a sphere does not move in a constant viscosity fluid. However, if viscosity increases towards the right as in figure 6(b), such that the two horizontal ends that are pulling the surface of the sphere outwards are in different viscosity environments, the force due to pressure $(\eta '/\eta _0)p^{\textit{Stokes}}$ is towards the right. Similarly, if viscosity increases upwards and considering that $p^{\textit{Stokes}}$ pushes the particle inwards on top and bottom, $(\eta '/\eta _0)p^{\textit{Stokes}}$ causes the particle to go downwards. Therefore, in a uniaxial extensional flow if the viscosity gradient lies along the compression (extensional) direction, the sphere will move towards the lower (higher) viscosity region. A sphere freely suspended in uniaxial extensional flow of uniform viscosity fluid has a saddle fixed point at the origin of the imposed flow. However, for the case of viscosity stratified fluid, the saddle point for the particle trajectory is shifted towards the lower viscosity fluid relative to the stagnation point of the imposed flow. The effect of stratification induced pressure (not shown) is similar to that of $(\eta '/\eta _0)p^{\textit{Stokes}}$ described above through figure 6 for both shear and extensional flow.

Figure 6.

Stokes pressure $p^{\textit{Stokes}}$ on the particle surface in (a) simple shear flow and (b) uniaxial extensional flow with viscosity increasing along the horizontal direction. Red (blue) represents a positive (negative) $p^{\textit{Stokes}}$ and light (dark) green represents higher (lower) $\eta '(=\eta -\eta _0)$ . Streamlines are of the Stokes velocity field for the respective flows. In (a) the imposed shear flow is towards the right on the top and left on the bottom half of the figure. The extensional flow in (b) is axisymmetric about the horizontal/ extensional axis. It approaches the particle’s center along the vertical (compression plane) and leaves it along the extensional axis.

In the next two sections, we will consider the motion of a freely suspended spheroid in settling due to gravity and a freely suspended neutrally buoyant particle in a linear flow field due to stratification-induced forces and torques.

5.1. Spheres ( $\kappa =1$ )

The rotational and translational velocities of a sedimenting sphere in a linearly stratified fluid are coupled via (3.14) and (3.16) which in the limit $\kappa \rightarrow 1$ leads to angular and translation velocities $\boldsymbol {\omega }_{\textit{particle}}$ and $\mathbf {u}_{\textit{particle}}$ given by

(5.1) \begin{align} \boldsymbol {\omega }_{\textit{particle}}&=\beta \frac {\mathbf {g}\times \mathbf {d}}{24\pi l \eta _0^2}+\mathcal {O}(\beta ^2), \nonumber \\ \mathbf {u}_{\textit{particle}}&=\frac {\mathbf {g}}{6\pi l \eta _0}-\beta \frac {l^2}{3\eta _0}\boldsymbol {\omega }_{\textit{particle}}\times \mathbf {d}+\mathcal {O}(\beta ^2)=\frac {\mathbf {g}}{6\pi l \eta _0}+\mathcal {O}(\beta ^2). \end{align}

From the $\mathcal {O}(\beta )$ analysis conducted here we observe that the sedimenting velocity of the sphere in a linearly stratified fluid does not change from that in a uniform viscosity fluid. However, at $\mathcal {O}(\beta )$ the stratification induces a rotation to the particle. A fluid rotating relative to a sphere at an angular velocity $\boldsymbol {\omega }_{{flow}}$ experiences a stratification induced force proportional to $\beta \boldsymbol {\omega }_{{flow}}$ ((3.11) and figure 5 b). The relative rotation here arises at $\mathcal {O}(\beta )$ , and part of the stratification induced velocity at $\mathcal {O}(\beta ^2)$ lies along $\boldsymbol {\omega }_{\textit{particle}}\times \mathbf {d}\propto \mathbf {d(d} \cdot \mathbf {g})$ . In other words, beyond the formally valid $\mathcal {O}(\beta )$ effects from the calculation conducting here, we may expect the velocity of a sedimenting sphere to have a finite component proportional to $\beta ^2 d_g \mathbf {d}$ , where the parameter,

(5.2) \begin{equation} d_g=\frac {1}{|| \mathbf {g}||_2}\mathbf {d}\cdot \mathbf {g}, \end{equation}

measures the alignment between gravity and the viscosity variation direction ( $d_g=$ 0, −1 and 1 imply that viscosity increases perpendicular, opposite and towards the gravity direction, respectively). This parameter is qualitatively important in discussion throughout § 5 where it plays a role in the $\mathcal {O}(\beta )$ effect on spheroids. For spheres, when $d_g\ne 0$ and $\mathbf {d}$ is not aligned with gravity $\mathbf {g}$ , we may expect that a sedimenting sphere will fall in a curved path instead of straight line along gravity. Quantitative conclusions about this change in the sphere’s trajectory can not be made from the $\mathcal {O}(\beta )$ calculation conducted in this paper. Two spheres in an inertia-less fluid with uniform velocity fall with no relative motion. However, another impact of the $\mathcal {O}(\beta )$ rotation rate induced by stratification is likely to be a change in the relative motion of the spheres.

5.2. Spheroids with $\kappa \ne 1$

We discussed above that a sedimenting sphere (a spheroid with $\kappa =1$ ) starts rotating due to viscosity stratification at $\mathcal {O}(\beta )$ and may undergo a horizontal drift at $\mathcal {O}(\beta ^2)$ . However, a similar stratification induced rotation on a spheroid with $\kappa \ne 1$ leads to change in settling behavior at $\mathcal {O}(\beta )$ , as its settling velocity (even in a constant viscosity fluid) depends on its centerline orientation, $\mathbf {p}$ . As a consequence, for a non-spherical spheroid, the rotational-translational coupling due to stratification leads to novel sedimenting behavior even at $\mathcal {O}(\beta )$ obtained from equation (3.3) and (3.18). We consider $\eta _0=1$ in the results presented below. The non-linear dynamics defined by the equation (3.18) for the particle rotation suggests a rich set of behaviors where, depending upon the values of $t_2$ and $t_3$ defined in (3.17) and the parameter $d_g$ (5.2), the orientation phase space has neutral orbits, spirals and fixed points (saddle, stable and unstable). Our complete analysis of these dynamics is given in appendix B. We summarize our results below.

The parameters $t_2$ and $t_3$ are only dependent on the particle aspect ratio $\kappa$ , and this variation for a wide range of $\kappa$ is shown in figure 7 (the expressions for $t_2$ and $t_3$ for a spheroid are provided in (D6) and (D7)). We observe that, particle orientation dynamics (and hence the translational dynamics, which is coupled with $\mathbf {p}$ ) depend only on $\kappa$ and $d_g$ . Thus, in $\kappa -d_g$ phase space, we obtain qualitatively different orientation behaviors demarcated by the boundaries shown in figure 8.

Figure 7.

Variation of shape-dependent parameters $t_i,i\in [1,3]$ (defined in (3.17)) with aspect ratio, $\kappa$ , for the rotation of (a) oblate and (b) prolate spheroids during sedimentation in a stratified environment. The parameters $t_2$ and $t_3$ change the particle’s orientation, $\mathbf {p}$ , at $\mathcal {O}(\beta )$ (3.18) analyzed in § 5.2, whereas $t_1$ causes rotation of the particle about it centerline that alters $\mathbf {p}$ if higher orders in $\beta$ are included through full the rotation rate equation (3.18). The focal length is chosen such that the major axis of the particle is $l=1$ for each particle type irrespective of $\kappa$ . The terms are multiplied with $12\pi$ for oblate and $8\pi /\log (2\kappa )$ for prolate particles.

Figure 8.

Phase diagram in $d_g-\kappa$ space of the orientation dynamics described by equation (B1) for spheroidal particles, i.e., $t_2$ and $t_3$ given by equations (D6) and (D7). Inside the boundaries marked $t_2t_3=0$ , $t_2t_3$ is negative, and it is positive otherwise. Here, GS refers to the gravity-stratification plane, and the superscripts $\parallel$ and $\perp$ refer to the final orientation of the particle’s axis of symmetry lying parallel and perpendicular to the GS plane respectively. In the regions labeled with subscript “spiral” the orientation trajectory towards its equilibrium position/orbit occurs in a spiraling motion, instead of a monotonic drift towards these positions.

To describe the particle dynamics in more physical terms than in appendix B, it is useful to define a coordinate system aligned with the gravity vector $\mathbf {g}$ and the stratification direction $\mathbf {d}$ . In this coordinate system, two of the three orthogonal axes are within the gravity-stratification (GS) plane: one of the axes in this plane is along the unit vector along gravity, $\hat {\mathbf {g}}$ , and along the other axis, $\hat {\mathbf {e}}$ , the viscosity gradient is non-decreasing such that $\hat {\mathbf {g}}\cdot \hat {\mathbf {e}}=0$ and $\mathbf {d}\cdot \hat {\mathbf {e}}=||\mathbf {d}\cdot \hat {\mathbf {e}}||_2=\sqrt {1-d_g^2}\gt 0$ . The third orthogonal axis is normal to the gravity-stratification (GS) plane and is defined by the vector $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ . Thus, the gravity and stratification directions can be expressed as

(5.3) \begin{eqnarray} \mathbf {g}=||\mathbf {g}||\hat {\mathbf {g}},\hspace {0.05in}\mathbf {d}=d_g\hat {\mathbf {g}}+\sqrt {1-d_g^2}\hat {\mathbf {e}}, \end{eqnarray}

and the particle orientation vector is expressed in this newly defined coordinate system as

(5.4) \begin{eqnarray} \mathbf {p}=p_g\hat {\mathbf {g}}+p_e\hat {\mathbf {e}}+\sqrt {1-p_g^2-p_e^2}\hat {\mathbf {g}}\times \hat {\mathbf {e}}. \end{eqnarray}

The invariant objects (neutral orbits, spirals, limit cycles and stable/ unstable/ saddle fixed points) of the dynamical system either lie on the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis or within the GS plane. The sign of the parameter $s=d_g(t_3-t_2)$ also plays a key role. The key points from appendix B (the different regions, $L_1$ , $L_2$ and $R_i,i\in [1,6]$ , are shown in figure 8) are:

1. (i) Region $L_1$ : When gravity and stratification are perpendicular, i.e., $d_g=0$ , and $t_2t_3\gt 0$ , the particle’s orientation follows a non-uniform neutral orbit with a time period $2\pi /(\sqrt {t_2t_3}||\mathbf {g}||_2\beta /\eta _0^2)$ . This is similar to the Jeffery orbits observed for particles with Bretherton constant $|q_4/q_1|\lt 1$ in simple shear flows.

2. (ii) Region $L_2$ : If gravity and stratification are collinear, i.e. $d_g=\pm 1$ , the particle aligns with gravity if $s\lt 0$ and perpendicular to gravity if $s\gt 0$ .

3. (iii) Regions $R_1=GS_{{spiral}}^\perp$ and $R_2=GS_{{spiral}}^\parallel$ : For $t_2t_3\gt 0$ and $0\lt |d_g|\lt |d_g|_{{spiral}}=2\sqrt {t_2t_3}/t_2+t_3$ , the GS plane is a limit cycle (stable if $s\gt 0$ and unstable otherwise) and the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis is a spiral (stable if $s\lt 0$ and unstable otherwise). Therefore, the particle’s axis of symmetry spirals towards the gravity-stratification (GS) plane if $s\lt 0$ ( $R_2$ ) or towards the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis if $s\gt 0$ ( $R_1$ ).

4. (iv) Regions $R_3=GS_1^\perp$ and $R_4=GS_3^\parallel$ : When $t_2t_3\gt 0$ and $|d_g|\gt |d_g|_{{spiral}}$ , the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis is a stable node when $s\gt 0$ and unstable otherwise. Hence, for $s\gt 0$ , the particle aligns along the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis irrespective of the initial condition ( $R_3$ ). For $s\lt 0$ ( $R_4$ ), a stable fixed point occurs on the GS plane, attracting all particle orientation trajectories.

5. (v) Regions $R_5=GS_1^\parallel$ and $R_6=GS_2^\parallel$ : The axis $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ is a saddle node for $t_2t_3\lt 0$ , but trajectories approach the node faster along its stable direction than they depart along the unstable direction when $s\gt 0$ ( $R_5$ ). If $s\lt 0$ , the unstable direction is sampled faster ( $R_6$ ). In both scenarios the particle orients at a stable fixed point in the GS plane, but, in region $R_5$ , a larger proportion of the trajectories approach $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ before departing it towards the GS plane.

Since $s/d_g=t_3-t_2\gt 0$ for prolate spheroids and $s/d_g\lt 0$ for oblate spheroids (figure 7), the orientation dynamics behavior observed for a positive $d_g$ for prolate spheroids is qualitatively replicated by a negative $d_g$ in oblate spheroids. As discussed above, the GS plane has fixed points when $|d_g|\gt |d_g|_{{spiral}}$ . The fixed points come in pairs due to the fore-aft symmetry of the particle, so there are an even number (four) of fixed points on the GS plane. When the GS plane is a stable attractor, two of the fixed points are stable nodes and the other two are saddle nodes. In the case when the GS plane is an unstable attractor, two fixed points are saddle nodes and the other two unstable nodes. Figure 9 shows the contours of $p_g^{(0,2)}$ , i.e., the magnitude of the projection of the least unstable fixed points (stable nodes for a stable attractor on GS and saddle nodes for an unstable attractor on GS), $\pm \mathbf {p}^{\text {stable}}_{\text {GS}}$ , projected along the gravity direction, in $\kappa -d_g$ space, where,

(5.5) \begin{equation} \pm p_g^{(0,2)}=\pm \mathbf {p}^{{stable}}_{{GS}}\cdot \hat {\mathbf {g}}. \end{equation}

The 0 in the superscript refers to a fixed point, and 2 refers to the second of the three pairs of fixed points (listed in (B7)) of the corresponding dynamical system given by (3.18).

Figure 9.

Contours of $p_g^{(0,2)}$ in $d_g-\kappa$ phase space for spheroids. As defined in equation (B7), $p_g^{(0,2)}$ is the magnitude of the projection of the location of the stable fixed point (when it exists) within the gravity-stratification (GS) plane along the gravity direction $\hat {\mathbf {g}}$ . Yellow ( $p_g^{(0,2)}\approx 1$ ) indicates the stable fixed point location to be closer to $\hat {\mathbf {g}}$ and blue ( $p_g^{(0,2)}\approx 0$ ) indicates a location closer to $\hat {\mathbf {e}}$ , i.e., perpendicular to $\hat {\mathbf {g}}$ .

Figure 10.

Orientation dynamics in the region $L_1$ , i.e. $d_g=0$ and for (a) $\kappa =0.56$ , (b) $\kappa =0.7$ , (c) $\kappa =1.5$ , and (d) $\kappa =2.0$ . Each curve represent the trajectory of the orientation vector of the particle’s axis. The magnitudes of $\eta _0$ , $\beta$ and $\mathbf {g}$ change the rate of rotation along these orientation trajectories, but not their shape. Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g=0$ considered here the red axis ( $\hat {\mathbf {e}}$ ) correspond to the viscosity stratification direction.

Figure 11.

Translation dynamics of $\kappa =2.0$ particle with $d_g=0$ , $||\mathbf {g}||=1$ and $\beta =0.05$ and 0.5 compared with that for unstratified fluid ( $\beta =0$ ). Each curve represents the trajectory of the particle’s centroid. Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g=0$ considered here the red axis ( $\hat {\mathbf {e}}$ ) correspond to the viscosity stratification direction.

Figure 12.

Orientation dynamics in $R_1$ exemplified with (a) oblate particle, $\kappa =0.65$ with $d_g=-0.5$ and (b) prolate particle, $\kappa =1.5$ with $d_g=0.5$ . The magnitudes of $\eta _0$ , $\beta$ and $\mathbf {g}$ change the rate of rotation along these orientation trajectories, but not their shape. Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g\lt 0$ considered here, the stratification direction lies in the plane of red ( $\hat {\mathbf {e}}$ ) and green ( $\hat {\mathbf {g}}$ ) axes and the viscosity increases in the positive $\hat {\mathbf {e}}$ and negative $\hat {\mathbf {g}}$ directions.

Figure 13.

(a) Orientation dynamics (the magnitudes of $\eta _0$ , $\beta$ , and $\mathbf {g}$ change the rate of rotation along these orientation trajectories, but not their shape) and (b) Translation dynamics in $R_2$ for a prolate spheroid, $\kappa =1.5$ with $d_g=-0.5$ . Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g\lt 0$ considered here, the stratification direction lies in the plane of red ( $\hat {\mathbf {e}}$ ) and green ( $\hat {\mathbf {g}}$ ) axes and the viscosity increases in the positive $\hat {\mathbf {e}}$ and negative $\hat {\mathbf {g}}$ directions.

Figure 14.

(a) Orientation dynamics and (b) Translation dynamics in $R_2$ for an oblate particle, $\kappa =0.65$ with $d_g=0.5$ . Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g\gt 0$ considered here, the stratification direction lies in the plane of red ( $\hat {\mathbf {e}}$ ) and green ( $\hat {\mathbf {g}}$ ) axes and the viscosity increases in the positive $\hat {\mathbf {e}}$ and positive $\hat {\mathbf {g}}$ directions.

Figure 15.

Orientation dynamics in $R_3$ ( $\kappa =0.65$ , $d_g=-0.9$ ) and $R_4$ ( $\kappa =1.5$ , $d_g=-0.9$ ). Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g\lt 0$ considered here, the stratification direction lies in the plane of red ( $\hat {\mathbf {e}}$ ) and green ( $\hat {\mathbf {g}}$ ) axes and the viscosity increases in the positive $\hat {\mathbf {e}}$ and negative $\hat {\mathbf {g}}$ directions.

Figure 16.

Orientation dynamics in $R_5$ ( $\kappa =3.0$ , $d_g=0.9$ ) and $R_6$ ( $\kappa =1/3$ , $d_g=0.9$ ). Green, red and blue axes represent the direction of vectors $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ respectively. For $d_g\gt 0$ considered here, the stratification direction lies in the plane of red ( $\hat {\mathbf {e}}$ ) and green ( $\hat {\mathbf {g}}$ ) axes and the viscosity increases in the positive $\hat {\mathbf {e}}$ and positive $\hat {\mathbf {g}}$ directions.

Figures 10 to 16 show the particle trajectories in black, and axes $\hat {\mathbf {g}}$ , $\hat {\mathbf {e}}$ and $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ are shown in green, red and blue respectively. These figures depict qualitatively different orientation (and hence the resulting translation) dynamics corresponding to regions $L_1$ (figures 10 and 11 for orientation and translation trajectories, respectively) and $R_i,i\in [1,6]$ discussed above. Trajectories in $R_1$ are illustrated in figure 12 and those in $R_2$ in figures 13 and 14. Figure 15 shows orientation trajectories in regions $R_3$ and $R_4$ and figure 16 shows these trajectories in regions $R_5$ and $R_6$ .

The oscillatory behavior of orientation dynamics within $\text {R}_1\cup \text {R}_2\cup L_1$ requires part of the viscosity variation to lie perpendicular to gravity, i.e. $|d_g|\ne 1$ , and this region is centered around $\kappa =1$ , i.e., a sphere. Any vector on a sphere can act as the orientation vector, which must necessarily undergo a neutral periodic orbit as a settling sphere always experiences the same stratification-induced torque, leading to continuous rotation. Therefore, in the $\kappa -d_g$ space plot of figure 8, each point on the vertical dashed line $\kappa =1$ denotes neutral or closed periodic orbits, with the rotation rate decreasing in magnitude as $|d_g|$ increases from 0. At $|d_g|=1$ ( $\mathbf {d}\times \mathbf {g}=0$ ) the sphere experiences no stratification-induced rotation. The sphere line in figure 8, $\kappa =1$ , acts as a bifurcation boundary for different possible behaviors upon changing $\kappa$ .

When the viscosity variation is entirely perpendicular to gravity, i.e., along the dashed horizontal $d_g=0$ line, the condition $t_2t_3\gt 0$ implies closed non-uniform, initial condition dependent periodic orbits for particle orientation with a time period $2\pi \eta _0^2/(\sqrt {t_2t_3}||\mathbf {g}||_2\beta )$ . For spheroids, this occurs when $0.56 \lessapprox \kappa \lessapprox 2$ and is indicated as $L_1$ in figure 8. As discussed in appendix B, these closed orbits are analogous to the Jeffery orbits of a freely rotating spheroid in simple shear flow of a constant viscosity fluid. The GS plane in the stratification-induced rotation of a sedimenting spheroid is analogous to the shearing plane in Jeffery orbits.

Orientation trajectories for four different $\kappa$ in the regime $L_1$ are shown in figure 10 (multiple curves in each plot for a $\kappa$ represents different initial orientations). For a sphere, the orientation trajectories are concentric circles parallel to the GS plane and centered around the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis. For a prolate spheroid in this regime, $1\lt \kappa \lessapprox 2$ , the orientation trajectories deviate from circles and point downwards at the $\hat {\mathbf {g}}$ axis. On the contrary, for oblate particles within $L_1$ ( $0.56\lessapprox \kappa \lessapprox 1$ ) the orientation trajectories point downwards at the $\hat {\mathbf {e}}$ axis. Similar to the Jeffery orbits, for $\kappa \ne 1$ the particle spends different amounts of time in different parts of its orientation trajectory for sedimentation induced rotation within region $L_1$ .

In particular, the particle spends more time (not shown) in the region of the orientation trajectory that points towards the GS plane (similar to Jeffery orbits pointing towards the shearing plane). Therefore, prolate spheroids in the $L_1$ regime spend more of the time with their axes of symmetry aligned with the gravity, $\hat {\mathbf {g}}$ axis, and oblate spheroids spend more time with their axes of symmetry aligned with the viscosity stratification, $\hat {\mathbf {e}}$ axis. In other words, during the majority of their orientation trajectory, the prolate spheroid’s axis of symmetry and the oblate particle’s face are aligned with gravity. This has profound implications for the sedimenting direction of the particle as shown in figure 11, where we compare the motion of a $\kappa =2.0$ particle initially placed at two different non-zero initial angles relative to gravity for unstratified and stratified fluids.

The dashed grey lines in each of the panels of figure 11 show that for a uniform viscosity fluid, the particle sediments at a constant initial orientation-dependent angle relative to gravity. Therefore, along with vertical settling, a spheroid drifts horizontally. However, with a viscosity gradient perpendicular to gravity, the non-uniform periodic nature of the orientation trajectory ensures that the particle falls mostly along the gravity direction without drifting too far horizontally. There is an instantaneous drift, but it is centered about the initial location normal to gravity with the maximum drift bounded, because either the particle’s axis of symmetry (prolate) or its face (oblate) is aligned with the gravity direction for the majority of the time during its orientation trajectory. The time period of the orientation trajectory reduces with increasing $\beta$ , reducing the time spent in orientations other than when its axis of symmetry (prolate) or face (oblate) is aligned with gravity. Therefore, the translation trajectory of the particle’s centroid becomes more confined as $\beta$ , i.e., the magnitude of the viscosity gradient, increases (compare black curves in figure 11 a vs. 11 b). Similar behavior (not shown) is observed for an oblate spheroid where the orientation is mostly along the stratification direction (figure 10 a) which also leads to a translation direction that is more confined and aligned with gravity. This confinement effect is more pronounced for $\kappa$ close to 0.56 and 2 within the regime $L_1$ as closer to these $\kappa$ , the particle’s rotation rate is more non-uniform than for particles with $\kappa$ closer to 1.

When $d_g\ne 0$ , changing $\kappa$ from 1 causes a bifurcation in the orientation trajectories. While the boundaries between regions $R_1$ and $R_2$ ( $\kappa =1$ ) are characterized by neutral closed orbits, spiraling trajectories are observed within regions $R_1$ and $R_2$ , with the rate of spiraling increasing as $\kappa$ deviates from 1 in these regions. A similar bifurcation can also be observed as $d_g$ is altered from 0 when going from $L_1$ to $R_1$ and $R_2$ . A prolate particle ( $\kappa \gt 1)$ spirals towards an orientation perpendicular to the GS plane if viscosity increases along gravity ( $d_g\gt 0$ ), i.e., the region $\text {R}_1={GS}_{{spiral}}^\perp$ . Conversely, if $d_g\lt 0$ , it spirals towards the GS plane in the region $\text {R}_2={GS}_{{spiral}}^\parallel$ . For an oblate particle, region ${GS}_{{spiral}}^\perp$ occurs if $d_g\lt 0$ and ${GS}_{{spiral}}^\parallel$ if $d_g\gt 0$ . The switching of behaviors between prolate and oblate spheroids with the sign of $d_g$ is related to the reversed sign of $t_2-t_3$ (positive for oblate and negative for prolate spheroids) for these particles shown in figure 7. The orientation trajectories for a few starting orientations in region $R_1$ for both an oblate and prolate particle are shown in figure 12. The magnitudes of $\beta$ and $\mathbf {g}$ do not affect the shape of these trajectories, but they change the spiraling rate towards the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis. Since the particle orientation is ultimately aligned normal to gravity, the particle settles parallel to gravity after the initial transient effects (not shown). The duration of the transient is inversely proportional to the rate of spiraling, which itself is proportional to $\beta ||\mathbf {g}||$ . Figure 13(a) shows the orientation trajectories of a $\kappa =1.5$ spheroid within the $R_2$ region, where spiraling towards the GS plane is observed. A portion of each spiral in this figure points downwards towards the GS plane. Similar to the Jeffery orbits or the neutral orbits shown in figure 10 for stratification-induced rotation of a sedimenting spheroid in region $L_1$ , this downward pointing portion is the slowest part of the spiral. The projection of this bottleneck region in the GS plane (of the trajectories shown in figure 13 a) is slightly misaligned with the gravity, $\hat {\mathbf {g}}$ , axis for a $\kappa =1.5$ particle exemplifying prolate spheroids in $R_2$ . For oblate spheroids in $R_2$ , this projects close to, but slightly misaligned from, the $\hat {\mathbf {e}}$ axis (figure 14 a). When the particle orientation ultimately reaches the GS plane, it continues to tumble there, but in a non-uniform fashion. This has a direct consequence on the translation trajectories shown in figure 13(b) and 14(b). The bottleneck or slow region in each spiral leads to an independence of the sedimenting angle from the initial orientation compared to the constant viscosity case (grey lines). Irrespective of initial condition, a spheroid orients within the GS plane and due to the bottleneck, it falls at a similar average angle relative to gravity. The time period for each spiral is inversely proportional to $\beta$ , and at higher $\beta$ the time spent outside of the bottleneck region is smaller, leading to straighter trajectories at larger $\beta$ shown in figures 13(b) and 14(b).

One notable difference appears between the sedimentation of prolate and oblate spheroids in $R_2$ . This requires $d_g\gt 0$ for oblate and $d_g\lt 0$ for prolate spheroids, i.e., it requires stratification to be misaligned with gravity but viscosity to increase along gravity for oblate and decrease for prolate spheroids (figure 8). The direction of $\hat {\mathbf {e}}$ is perpendicular to $\hat {\mathbf {g}}$ and at least a part of viscosity increase is along $\hat {\mathbf {e}}$ (when $|d_g|\ne 1$ ). Prolate particles in $R_2$ migrate towards the positive $\hat {\mathbf {e}}$ (figure 13 b) while the oblate particles towards the negative $\hat {\mathbf {e}}$ (figure 14 b) axis. Therefore, the horizontal drift of prolate spheroids is towards higher viscosity fluid, while that of oblate spheroids is towards lower viscosity fluid in the region $R_2$ .

Within $R_2$ , as a particle is made less spherical, i.e., as $\kappa$ deviates further from one, the spiral points more downwards and, as mentioned above, the spiraling rate increases. Towards the edge of $R_2$ further from $\kappa =1$ , the downwards pointing part of the spirals almost completely touches the GS plane such that at the edge between $R_2$ to $R_4$ in figure 8 another bifurcation is observed. Hence, in $R_4$ , spirals no longer exist, and instead, two fixed points appear on the GS plane as the GS plane bifurcates from a stable limit cycle to a stable sub-space, and the axis $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ bifurcates from an unstable spiral to an unstable node. The orientation trajectories of a $\kappa =1.5$ particle with $d_g=-0.9$ are shown in figure 15(b) and comparing them with figure 13(a) shows the bifurcation. The stable fixed point is closer to the $\hat {\mathbf {g}}$ axis for prolate spheroids, in continuation with the bottleneck region from $R_2$ for such particles. This alignment of the particle closer to gravity is also illustrated by a $p_g^{(0,2)}\approx 1$ (5.5) for prolate spheroids in the region corresponding to $R_4$ in figure 9. A similar topology of orientation trajectories is observed for oblate spheroids (not shown), but the location of the stable fixed point within the GS plane is more sensitive to the values of $\kappa$ and $d_g$ (as shown by the rapidly changing $p_g^{(0,2)}$ in the $R_4$ region for $\kappa \lt 1$ particles than $\kappa \gt 1$ in figure 9). An oblate particle is aligned closer to $\hat {\mathbf {e}}$ axis near the $R_2-R_4$ boundary ( $p_g^{(0,2)} \approx 0$ near this boundary in figure 9), and as $\kappa$ reduces (such that the particle is more disc-like), it moves towards the $\hat {\mathbf {g}}$ axis (indicated in figure 9 by $p_g^{(0,2)}$ being closer to 1 than 0 upon reducing $\kappa$ for oblate particles in $R_4$ ). A similar effect is found upon increasing $d_g$ from 0 to 1 for oblate spheroids. The translation trajectories in region $R_4$ (not shown) follow a similar initial orientation independent trend as discussed above for $R_2$ for prolate and oblate spheroids, except that particles fall in a (perfectly) straight line as their orientation reaches a steady state in $R_4$ .

The bifurcations that occur between $R_1$ and $R_3$ are similar to the $R_2$ to $R_4$ bifurcation just discussed. Region $R_1$ , discussed earlier, is similar to $R_2$ but with spiraling away from the GS plane. Therefore, similar to the bifurcation from $R_2$ to $R_4$ discussed above, a bifurcation occurs at the edge between $R_1$ and $R_3$ in figure 8 as the spiraling rate reaches infinity. Figure 15(a) illustrates the orientation trajectories of a $\kappa =0.65$ particle with $d_g=0.9$ , a case within $R_3$ . In $R_3$ , coming from $R_1$ (comparing figure 12 a with 15 a), the GS plane changes from an unstable limit cycle to an unstable subspace with two fixed points (one saddle and the other unstable node), and the axis $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ bifurcates from a stable spiral node to a stable node. Since the ultimate orientation is perpendicular to gravity, the particle falls along gravity in $R_3$ (not shown).

When gravity and stratification are collinear, i.e., $|d_g|=1$ , a spheroid will orient either parallel or perpendicular to gravity irrespective of its starting orientation. In either case, it sediments along the gravity direction. A prolate spheroid orients perpendicular to gravity if viscosity increases along gravity, i.e., $d_g=1$ , and parallel if $d_g=-1$ . An oblate spheroid shows the opposite trend relative to the sign of $d_g$ . Regions where particles orient perpendicular to gravity are labeled as $L_2$ in figure 8. This behavior is similar to that observed by Anand & Narsimhan (Reference Anand and Narsimhan2024) and can also be deciphered for prolate spheroids from the formulae of Gong et al. (Reference Gong, Shaik and Elfring2024).

In the region $R_5\cup R_6$ , demarcated by vertical lines at $\kappa \lessapprox .56$ and $\kappa \gtrapprox 2$ and horizontal lines for $d_g=1, \kappa \gt 1$ and $d_g=-1, \kappa \lt 1$ , the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis is a saddle node. Transitioning from $R_4$ to $R_6$ (comparing figures 15 b and 16 b) involves a bifurcation where the saddle node on the GS plane becomes an unstable node, and the unstable node at the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis becomes a saddle. This does not alter the topology around the stable fixed point on the GS plane. Similarly, transitioning from $R_3$ to $R_5$ (comparing figures 15 a and 16 a) causes the stable node at $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis to bifurcate into a saddle node, and the saddle node on the GS plane becomes a stable node without altering the topology around the unstable fixed point on the GS plane. Therefore, in both $R_5$ and $R_6$ , the particle ultimately orients within the GS plane, but in a different fashion. In $R_5$ , the stable direction of the saddle at $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ is sampled faster than its unstable direction, so the spheroids first approach a plane that includes the $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ axis before traversing towards the GS plane, as shown in figure 16(a). In $R_6$ , the unstable eigenvalue of the saddle point at $\hat {\mathbf {g}}\times \hat {\mathbf {e}}$ has a greater magnitude than the stable eigenvalue, so the particle reaches its final location in the GS plane faster than in $R_5$ (figure 16 b).

The contours of $p_g^{(0,2)}$ (equation (5.5)) in figure 9 depict the final orientation of the spheroid’s axis of symmetry relative to the gravity direction. As shown in this figure, for $d_g\lt 0$ , the final orientation in region $R_5\cup R_6$ is closer to $\hat {\mathbf {g}}$ for prolate and to $\hat {\mathbf {e}}$ for oblate spheroids. The final orientation is more sensitive to $d_g$ and $\kappa$ within $R_5\cup R_6$ for $d_g\gt 0$ , where both oblate and prolate spheroids orient closer to $\hat {\mathbf {g}}$ (yellow region in figure 9) as $\kappa$ decreases. Increasing $d_g$ , i.e., increasing the alignment between stratification and gravity, within $R_5\cup R_6$ for $d_g\gt 0$ makes prolate and oblate spheroids orient towards $\hat {\mathbf {e}}$ (blue) and $\hat {\mathbf {g}}$ (yellow), respectively. The influence of stratification-induced rotation on sedimentation in $R_5$ and $R_6$ (not shown) is similar to that discussed earlier for $R_3$ and $R_4$ . The particle falls ultimately in a straight, initial condition independent path, with its angle relative to gravity depending on the location of the stable fixed point on the GS plane, $p_g^{(0,2)}$ .

At the intersection of the regions marked $R_5$ and $R_6$ when gravity is perpendicular to stratification, $d_g=0$ , for particles with aspect ratio much different from a sphere Anand & Narsimhan (Reference Anand and Narsimhan2024) also found a single stable orientation that is independent of the initial orientation. The value of $p_g^{(0,2)}$ , from our calculations shown in figure 9 varies from about 0.998 to 0.931 for $2\lessapprox \kappa \lessapprox 10$ compares favorably with an equivalent value of about $|\cos (57\pi /74)|\approx 0.9\approx 0.94$ displayed in figure 14(a) of Anand & Narsimhan (Reference Anand and Narsimhan2024). However, for $0.1\lessapprox \kappa \lessapprox 0.5$ , we find $p_g^{(0,2)}$ to vary between 0.49 and 0.12, whereas for the same $\kappa$ range values displayed in figure 14(b) of Anand & Narsimhan (Reference Anand and Narsimhan2024) vary between 0.5 and 0.25. Furthermore, Anand & Narsimhan (Reference Anand and Narsimhan2024) find stable orientations for oblate spheroids with $\kappa$ as low as $0.01$ , which is well within the $L_1$ (figure 8) region where neutral orbits are predicted by our calculations. Thus, the final orientation predicted by our calculations for large $\kappa$ spheroids agrees well with the previous investigations, but the orientation predicted for oblate spheroids is different. A future numerical investigation of oblate spheroids, akin to that conducted for prolate spheroids in § 4.2.1 (black markers in figure 5), could serve as an independent validation of our findings.

6.1. Spheres ( $\kappa =1$ )

A freely suspended sphere in a linear flow of uniform viscosity fluid simply rotates with the imposed fluid rotation, $\boldsymbol {\omega }_\infty =-0.5\boldsymbol {\epsilon }:\boldsymbol {\Omega }$ , where $\boldsymbol {\Omega }$ is the anti-symmetric part of the imposed velocity gradient, and translates with the local velocity of the imposed flow, $\mathbf {u}_{{fluid}}$ . Coupling between the translational and rotational motion due to viscosity stratification leads to the following particle motion,

(6.1) \begin{align} \begin{split}&\mathbf {u}_{\textit{particle}}=\mathbf {u}_{{fluid}}+\frac { \beta l^2}{\eta _0}\mathbf {E}\cdot \mathbf {d}+\mathcal {O}(\beta ^2). \end{split} \end{align}

In a simple shear flow, with strain rate $\dot {\gamma }$ , the velocity of a sphere relative to the viscosity stratified fluid is ${\dot {\gamma }\beta l^2}/({2\eta _0})\begin{bmatrix}d_2,&d_1,&0\end{bmatrix}$ , where the directions 1, 2 and 3 are respectively in the flow, gradient and vorticity direction of the imposed flow. Thus, for magnitude $\beta$ and direction $\mathbf {d}=\begin{bmatrix} d_1,&d_2,&d_3 \end{bmatrix}$ of viscosity stratification, a sphere of radius $l$ can be moved across the flow streamlines with a speed ${\dot {\gamma }\beta l^2d_1}/({2\eta _0}).$

Observing the effect of stratification in simple shear flow, we may conjecture this effect in particle-filled heated Couette and Poiseuille flows. If the viscosity is uniform, the position of the particles relative to the walls does not change (figures 17 a and 18 a). If the inlet of a channel is at a different temperature than the outlet such that the viscosity increases along the channel length, in Couette flow, as schematically depicted in figure 17, dispersed spheres will migrate towards the wall that moves in the same direction as the increasing viscosity. In the case of Poiseuille flow, if the viscosity increases along the flow direction, the particles migrate towards the center of the channel (figure 18 b). If viscosity increases in the opposite direction, they move towards the walls due to the stratification-induced force and local shear (figure 18 c). These hypotheses ignore inter-particle hydrodynamic interactions and assume the particles to be small enough such that locally they observe a simple shear flow in an unbounded fluid.

Figure 17.

Schematic of spheres dispersed in Couette flow of (a) unstratified and (b) stratified fluid with viscosity increase towards right. Based on the results of a single particle in an unbounded fluid, the particles migrate towards the wall that moves in the direction of increasing viscosity.

Figure 18.

Schematic of spheres dispersed in Poiseuille flow of (a) unstratified and stratified fluid with viscosity (b) increase and (c) decrease towards the flow direction.

In the case of uniaxial extensional flow, spheres move towards more viscous fluid if viscosity varies along the extensional axis, and towards less viscous fluid if viscosity varies along the compression axis. n the case of uniaxial extensional flow, the stratification-induced relative velocity is $\beta l^2/\eta _0 \mathbf {E}\cdot \mathbf {d}+\mathcal {O}(\beta ^2)$ . Therefore, spheres move towards more viscous fluid if viscosity varies along the extensional axis, and towards less viscous fluid if viscosity varies along the compression axis.

6.2. Spheroids with $\kappa \ne 1$

Next, we consider the effect of viscosity stratification on the translation trajectories of a spheroid with $\kappa \ne 1$ , where the particle’s orientation also influences these dynamics (3.21), resulting in a more complex effect. The translation velocity of the particle relative to the fluid for a spheroid is given by

(6.2) \begin{align} \mathbf {u}_{\textit{particle}}-\mathbf {u}_{{fluid}}&=\frac {\beta }{\eta _0} (m_1 \mathbf {E} + m_2 (\mathbf {p}\cdot \mathbf {E}\cdot \mathbf {p})\mathbf {I} +(m_2+m_3) (\mathbf {p}\cdot \mathbf {E}\cdot \mathbf {p})\mathbf {pp}\nonumber \\&\quad + m_4 [\mathbf {E}\cdot \mathbf {pp}+\mathbf {pp}\cdot \mathbf {E}])\cdot \mathbf {d}+\mathcal {O}(\beta ^2), \end{align}

where $m_i, i\in [1,4]$ as a function of $\kappa$ are displayed in figure 19, and their analytical expressions are in (D10) (the functional dependence of $m_i, i\in [1,4]$ on $f_1$ , $f_3$ , $q_1$ , $q_3$ and $F_{ij}^{\Gamma _k},i,j\in [1,4], k\in [1,8]$ can be ascertained by comparing the RHS of the above equation with that of (3.21)). Here, only the first effects of the stratification-induced force on a rotating particle are accounted for, which arise from the rotation rate of the particle’s axis of symmetry up to $\mathcal {O}(\beta )$ . Thus, the particle rotates along Jeffery orbits (3.8).

Figure 19.

Variation of coefficients, $m_i, i\in [1,4]$ in $\mathcal {O}(\beta )$ translation velocity (6.2) in linear flows around spheroids with aspect ratio $\kappa$ and major radius 1.

6.2.1. Spheroids in uniaxial extensional flow

Consider an extensional flow such that the imposed velocity gradient is $\Gamma _{ij}=\delta _{i1}\delta _{j1}-0.5(\delta _{i2}\delta _{j2}+\delta _{i3}\delta _{j3})=E_{ij},$ resulting in $\mathbf {u}_{{fluid}}=\mathbf {r}\cdot \mathbf {E}$ . According to the particle rotation rate in a uniform viscosity fluid given by (3.8), a spheroid obtains a steady state orientation such that $\mathbf {E}\cdot \mathbf {p}=\alpha \mathbf {p}$ , for a constant $\alpha$ . Two possible values of $\alpha$ are + 1 and -1/2, where one corresponds to a stable fixed point and the other to an unstable fixed point, depending upon the sign of $q_4/q_1$ . In other words, a prolate spheroid orients with its major axis along the extensional axis ( $\alpha =1$ is the stable case), and an oblate particle orients its face in a plane consisting of the extensional axis and one of the compressional directions ( $\alpha =-1/2$ is the stable case). At this steady orientation, the spheroids move with the relative translation velocity,

(6.3) \begin{eqnarray} \begin{split} \mathbf {u}_{\textit{particle}}-\mathbf {u}_{{fluid}}=\frac {\beta }{\eta _0 f_1}[m_1\mathbf {E}\cdot \mathbf {d}+ \alpha (m_2\mathbf {d}+(m_2+m_3+m_4)\mathbf {d}\cdot \mathbf {pp})]+\mathcal {O}(\beta ^2),\,\,\,\, \end{split} \end{eqnarray}

where $\alpha =1$ for prolate and –1/2 for oblate spheroids. If the stratification $\mathbf {d}$ is directed along the extensional axis, the translation of the particles follows,

(6.4) \begin{equation} \mathbf {u}_{\textit{particle}}-\mathbf {u}_{{fluid}}=\frac {\beta }{\eta _0 f_1}\begin{cases} (m_1+2m_2+m_3+m_4)\mathbf {d},\hspace {0.1in}\kappa \gt 1\\ (m_1-m_2/2)\mathbf {d},\hspace {0.1in}\kappa \lt 1 \end{cases}. \end{equation}

In the case where the stratification $\mathbf {d}$ is directed along the compressional axis,

(6.5) \begin{equation} \mathbf {u}_{\textit{particle}}-\mathbf {u}_{{fluid}}=\frac {\beta }{\eta _0 f_1}\begin{cases} (-m_1/2+m_2)\mathbf {d},\hspace {0.1in}\kappa \gt 1\\ -1/2(m_1+2m_2+m_3+m_4)\mathbf {d},\hspace {0.1in}\kappa \lt 1 \end{cases}. \end{equation}

The coefficients of $\mathbf {d}$ inside the curly brackets in (6.4) are positive, and those in (6.5) are negative for all $\kappa$ , with the largest magnitude occurring for a sphere (not shown). Hence, a viscosity gradient along the extensional axis makes prolate, oblate and spherical particles move towards more viscous fluid, and a viscosity gradient along a compression axis moves the spheroids towards the lower viscosity region.

6.2.2. Spheroids in simple shear flow

In simple shear flow with velocity gradient $\Gamma _{ij}=\delta _{i2}\delta _{j1}$ , such that at the particle’s center $u_{i,{fluid}}={r}_{2,{spheroid}}$ , and considering just the Jeffery rotation of the spheroid’s orientation, $\mathbf {p}=\begin{bmatrix}p_1&p_2&p_3 \end{bmatrix}^{\rm T}$ , the relative translation velocity is,

(6.6) \begin{align} \begin{split} u_{i,{particle}}-u_{i,{fluid}}=& \frac {\beta }{\eta _0} [0.5 m_1 (d_1\delta _{i2}+d_2\delta _{i1})+p_1p_2 (m_2d_i+(m_2+m_3)p_jd_jp_i) \\&+m_4((p_1\delta _{i2}+p_2\delta _{i1})p_kd_k+p_i(p_2d_1+p_1d_2))]+\mathcal {O}(\beta ^2). \end{split} \end{align}

The Jeffery orbits, or orientation trajectories, of spheroids in simple shear flow of a uniform viscosity fluid can be broadly classified into four types: (a) log-rolling, (b) tumbling, (c) flipping, and (d) wobbling orbits. The effect of stratification on the translation of a particle undergoing these orientation trajectories, along with the orbit descriptions, are discussed below.

Log-rolling orbits

A spheroid initially oriented along the vorticity direction (i.e. $p_1=p_2=0$ and $p_3=1$ ) does not change its orientation with time but simply rolls about its axis at angular velocity equal to half of the shear rate. This motion is thus referred to as log-rolling, where the translation velocity is

(6.7) \begin{equation} u_{i,{particle}}-u_{i,{fluid}}=0.5 m_1\frac {\beta }{\eta _0} ( (d_1\delta _{i2}+d_2\delta _{i1}))+\mathcal {O}(\beta ^2). \end{equation}

Since, $m_1$ does not change sign with $\kappa$ (figure 19), the effect of stratification on a log-rolling spheroid is qualitatively similar to the motion of a sphere discussed in § 6.1. The magnitude of the stratification-induced velocity is largest for the sphere. In the case of a log-rolling spheroid, if stratification is along the vorticity direction of the imposed flow, i.e., $\mathbf {d}=\begin{bmatrix} 0&0&1 \end{bmatrix}$ , the stratification-induced velocity is zero. Otherwise, the trajectory of the particle’s centroid initially located at $\begin{bmatrix}0,&0,&0\end{bmatrix}^{\rm T}$ moves along the parabola $0.5m_1 t \beta \begin{bmatrix} 0.5d_1t + d_2,& d_1, & 0\end{bmatrix}^{\rm T}$ . Hence, a particle placed at the origin moves along the flow direction if stratification is along the velocity-gradient direction of the imposed flow. For stratification along the flow direction, the particle is displaced along the gradient direction and then also swept along the flow direction by the imposed flow. Since $m_1\gt 0$ for all spheroids (figure 19), the particle is ultimately swept towards higher viscosity regions.

Tumbling orbits

A spheroid with initial orientation in the flow-gradient plane ( $p_3=0$ ) remains there and continues to tumble in this tumbling orbit with a time period $T_{\text {Jeffery}}=2\pi (\kappa +1/\kappa )$ (normalized with the shear rate). The stratification-induced translation velocity of a tumbling spheroid is,

(6.8) \begin{align} &\begin{bmatrix} u_{1,{particle}}-u_{1,{fluid}}\\ u_{2,{particle}}-u_{2,{fluid}}\\ u_{3,{particle}}-u_{3,{fluid}} \end{bmatrix}\nonumber \\&\quad =\frac {\beta }{\eta _0} \begin{bmatrix}0.5 m_1 d_2 +m_2d_1p_1p_2+(p_1d_1+p_2d_2)(p_1^2p_2+m_4(p_1+p_2))\\ 0.5 m_1 d_1 +m_2d_2p_1p_2+(p_1d_1+p_2d_2)(p_1p_2^2+m_4(p_1+p_2))\\ m_2d_3p_1p_2 \end{bmatrix}. \end{align}

Considering first the simpler case of stratification along the vorticity direction, i.e., $d_1=d_2=0, d_3=1$ , the particle moves along the viscosity gradient direction with a rate $\beta p_1p_2m_2$ . From figure 19, we observe that $m_2$ changes sign at $\kappa =1$ ; therefore, at a particular orientation within the shearing plane, an oblate spheroid moves in the opposite direction to the prolate spheroid due to stratification effects. However, either particle’s motion along the vorticity direction in this case is reversed when the particle is in the extensional quadrant ( $p_1p_2\gt 0$ ) compared to when it is in the compression quadrant ( $p_1p_2\lt 0$ ). Since a spheroid spends equal amounts of time in these quadrants, its centroid’s translation velocity oscillates about the initial value. Therefore an oblate (prolate) particle, started with its centerline along the gradient direction, makes an excursion towards the higher (lower) viscosity region and periodically returns to its original location. A few examples of the particle’s normalized displacement from its original location in the vorticity direction $\Delta z/\beta$ are shown in the left panel of figure 20.

Not only does $m_2$ decrease for aspect ratios further from the sphere ( $\kappa \ll 1$ and $\kappa \gg 1$ in figure 19), but for the majority of Jeffery orbits, a prolate spheroid spends most of the time oriented along the flow direction ( $p_1=1$ , $p_2=p_3=0$ ) and an oblate with its face in flow-vorticity plane ( $p_2=1$ , $p_1=p_3=0$ ), i.e., the orientations where the induced velocity $\beta p_1p_2m_2$ is zero. This implies that for $\kappa \ll 1$ and $\kappa \gg 1$ , the amplitude of the aforementioned oscillations in the position of the particle’s centroid is smaller for these extreme values of $\kappa$ . Oscillations are zero for a sphere, and hence an optimal $\kappa$ exists for both the $\kappa \lt 1$ and $\kappa \gt 1$ regimes where the particle oscillates with the largest amplitude. The right panel of figure 20 shows the $\kappa$ variation of the location of the particle’s maximum excursion from its original location in the vorticity direction for a viscosity increase along the positive vorticity direction.

Figure 20.

(a) Excursion of spheroids (with major axis, $l=1$ ), $\Delta z/\beta$ , tumbling in the flow-gradient (shearing) plane of simple shear flow along the vorticity direction when viscosity varies along the vorticity direction for different aspect ratio, $\kappa$ . (b) Maximum excursion from the initial location as a function of $\kappa$ . The largest magnitude of $\Delta z$ toward higher viscosity occurs for $\kappa \approx 0.5$ and towards less viscous fluid for $\kappa \approx 20$ .

Similar to the log-rolling case discussed above, for a particle in a tumbling orbit, if the stratification is in the flow-gradient plane (i.e., $d_3=0$ ) the particle will only translate within this plane as from equation (6.8) ${\rm d}{r}_{3,{spheroid}}/{\rm d}t|_{{tumbling}}=0$ . However, unlike the log-rolling scenario, here a stratification purely in the gradient direction ( $d_1=d_3=0$ ) will cause the particle to not only move along the flow but also along the gradient or the shearing direction. The trajectories of $\kappa =1/2, 1, 2$ and 10 particles over one respective Jeffery time period are shown in figure 21(a) when the viscosity increases along the gradient direction and $\beta =0.1$ . Figure 21(b) shows the trajectory for the same parameters, but with viscosity increasing along the flow direction.

Figure 21.

Translation trajectories for various spheroids (with major axis, $l=1$ ) when viscosity varies along (a) gradient direction and (b) flow direction. The viscosity gradient magnitude, $\beta =0.1$ .

Flipping and wobbling orbits

For an initial orientation close to but not on the flow-gradient plane, a spheroid with either large (prolate) or small (oblate) $\kappa$ spends most of its Jeffery orbit near the flow direction. Within a small time frame, it flips from one side of the gradient-vorticity plane to the other, and during this time, it traverses a larger three-dimensional orientation space. These are termed flipping orbits. For initial orientations close to but not on the vorticity axis, Jeffery orbits are also three-dimensional, but in these wobbling trajectories, the rotation of the particle throughout its orbit is more uniform than in the flipping orbits. Due to the three-dimensional rotation of the particle, shown in figure 22(a) for $\kappa =10$ , the stratification-induced force translates the particle in a three-dimensional manner.

Figure 22.

Translation trajectories ((b), (c) and (d)) of a $\kappa =10$ spheroid with $\beta =0.1$ for different flipping (b, c) and wobbling (d) orientation trajectories shown in (a) and different directions of viscosity stratification,  $\mathbf {d}$ .

In § 5.1 we speculated that beyond the $\mathcal {O}(\beta )$ effects considered in this paper, a sedimenting sphere is expected to follow a curved settling trajectory due to the stratification induced force on the rotating particle. There the rotation induced at $\mathcal {O}(\beta )$ might be expected to induce a stratification force and modify the particle’s translation motion at $\mathcal {O}(\beta ^2)$ . Similarly, here, for a spheroid in simple shear flow a stratification induced torque at $\mathcal {O}(\beta ^2)$ may lead to a rotational motion of the particle that deviates from Jeffery orbits.