2.1. Physical model of a spheroid rotating in a linear shear flow

In this study, the flow we considered is formulated by $u=\kappa y+ U_s$ , which is a synthesis of a linear shear flow ( $u_\kappa =\kappa y$ ) and a uniform flow ( $U_s$ ), as depicted in figure 1. Here, $\kappa$ is the shear rate and $u$ is the streamwise velocity. Moreover, the corresponding velocity gradient is $\unicode{x1D63C}^*$ , where the components are zero except for $A_{xy}^*=\kappa$ . Hence, the strain-rate tensor $\unicode{x1D64E}^{\;*}$ and the rotation tensor $\unicode{x1D64A}^*$ are expressed as follows:

(2.1) \begin{equation} \unicode{x1D64E}^{\;*}=\left [\begin{matrix} 0 & \kappa /2 & 0 \\ \kappa /2 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}\right ]\!,\quad \unicode{x1D64A}^{\;*}=\left [\begin{matrix} 0 & \kappa /2 & 0 \\ -\kappa /2 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}\right ]\!. \end{equation}

The rotation of fluid $\boldsymbol \varOmega ^*$ is $[0,0,-\kappa /2]^T$ , which is perpendicular to the flow-gradient plane, i.e. $x\!{-}\!y$ plane and, therefore, coincides with the $z$ -axis. Mostly, the $z$ -axis is also called the vortical axis.

Figure 1.

Schematic of the synthesis shear flow of $u=\kappa y+U_s$ .

In this study, we consider spheroidal particles, which are axisymmetric, as shown in figure 2. The ratio of the symmetry axis ( $2a$ ) to the equatorial axis ( $2b$ ) is denoted as aspect ratio $\lambda =a/b$ . Meanwhile, we also use shape parameter $\varLambda =(\lambda ^2-1)/(\lambda ^2+1)$ to regularize shape within $(-1,1)$ . Therefore, prolate spheroids have $\lambda \gt 1(0\lt \varLambda \lt 1)$ , while oblate spheroids have $0\lt \lambda \lt 1(-1\lt \varLambda \lt 0)$ . For the spherical one, the aspect ratio is $\lambda =1$ and the corresponding shape parameter $\varLambda$ is equal to zero. The orientation of a spheroid is defined by the direction of its symmetry axis, denoted by the unit vector $\boldsymbol{n}$ . Two Euler angles, $\varphi$ and $\vartheta$ , are introduced as shown in figure 2(c). The azimuthal angle $\varphi$ measures the angle between the $x$ -axis and the projection of $\boldsymbol{n}$ onto the $x$ – $y$ plane, while the polar angle $\vartheta$ denotes the angle between $\boldsymbol{n}$ and the vorticity (or $z$ ) axis. In this study, the centre of a spheroidal particle is set at $y=0$ , where a slip velocity between particle and fluid exists, denoted as $u_s\equiv U_s$ . This configuration is consistent with earlier studies (Li et al. Reference Li, Cui, Huang, Zhao and Xu2018; Cui et al. Reference Cui, Liu, Qiu and Zhao2024), which is similar to a spheroid settling in a vertical shear flow with a gravity along to the flow. Nevertheless, the slip velocity is constant in the present study.

Figure 2.

Schematic of (a) prolate and (b) oblate spheroids, where green arrows represent the direction of the symmetry axes. (c) Euler angles ( $\varphi$ , $\vartheta$ ) of symmetry axis.

In the current study, there are two Reynolds numbers. One is the particle Reynolds number based on the slip velocity, namely ${\textit{Re}}_{\!p}=u_sl/\nu$ , where $l$ is the semilength of the major axis and $\nu$ is the fluid kinematic viscosity. The other one is the shear Reynolds number defined as ${\textit{Re}}_\kappa =\kappa l^2/\nu$ , reflecting the effect of fluid inertia from shear. Throughout this study, we primarily consider cases with ${\textit{Re}}_{\!p}, {\textit{Re}}_\kappa \lt 1$ to minimize the influence of high-order fluid inertia. The non-dimensional parameters that correspond to particle inertia are the rotational and translational Stokes numbers. Here, only the rotational Stokes number, denoted as $St_r$ , should be taken into account, because the translational motions are constrained in this study. Note that rotational Stokes number is defined as $St_r=D_{{\kern-1pt}pf}\kappa l^2/\nu \equiv D_{{\kern-1pt}pf} {\textit{Re}}_\kappa$ , where $D_{{\kern-1pt}pf}$ is the density ratio between particle and fluid. For the sake of simplicity, hereinafter, the Stokes number only represents the rotational Stokes number. In addition, we focus on large density ratios ( $D_{{\kern-1pt}pf} \gg 1$ ), such that the fluid inertia associated with the added mass can be neglected.

2.2. Angular dynamics of spheroids

According to Lundell & Carlsson (Reference Lundell and Carlsson2010) and Challabotla et al. (Reference Challabotla, Nilsen and Andersson2015), the angular dynamics of spheroids are computed in the particle coordinate by the following equation:

(2.2) \begin{equation} \unicode{x1D644}^{\,\prime } \frac {\mathrm{d}{{\boldsymbol{\omega}^*}^{\prime }}}{\mathrm{d}t^*}+{\boldsymbol{\omega}^*}^{\prime }\times \unicode{x1D644}^{\;\prime }\boldsymbol{\cdot }{\boldsymbol{\omega}^*}^{\prime } = {\boldsymbol{T}}^{\prime }. \end{equation}

Here, $\unicode{x1D644}^{\;\prime}$ is the moment of inertia, ${\boldsymbol{\omega}^*}^{\prime }$ is the rotation of particle and $\boldsymbol{T'}$ is the resultant torque acting on spheroid. Note that $(\boldsymbol{\cdot })'$ denotes variables expressed in the particle coordinate, where the vectors (or tensors) are defined with respect to the laboratory reference frame but their components are written in a Cartesian basis attached to the particle. The left-hand side of (2.2) represents the dynamics of the rigid body, referring to the particle inertia. For a spheroidal particle, the moment of inertia is expressed as

(2.3) \begin{equation} \unicode{x1D644}^{\;\prime}=\frac {1}{5}m_{\!p}\left [\begin{matrix} 2b^2&0&0\\ 0&a^2+b^2&0\\ 0&0&a^2+b^2 \end{matrix}\right ]\!, \end{equation}

where $m_{\!p}$ is the mass of spheroidal particle, i.e. $m_{\!p}=({4}/{3})\rho _{\!p}\pi a b^2$ . Here, $\rho _{\!p}$ is the density of the particle. To describe quantities in different frames, a vector $\boldsymbol{p}$ in the laboratory reference frame can be transformed from the particle coordinate by $\boldsymbol{p} = \unicode{x1D64D} \boldsymbol{\cdot }\boldsymbol{p'}$ , while a matrix $\unicode{x1D63D}$ in the laboratory reference frame can be transformed from the particle coordinate via $\unicode{x1D63D}=\unicode{x1D64D} \boldsymbol{\cdot }\unicode{x1D63D}' \boldsymbol{\cdot }\unicode{x1D64D}^T$ . Here, $\unicode{x1D64D}$ is a rotation matrix, which is defined as follows:

(2.4) \begin{equation} \unicode{x1D64D}=\left [\begin{matrix} e_0^2+e_1^2-e_2^2-e_3^2 & 2(e_1e_2-e_0e_3) & 2(e_1e_3+e_0e_2) \\ 2(e_1e_2+e_0e_3) & e_0^2-e_1^2+e_2^2-e_3^2 & 2(e_2e_3-e_0e_1) \\ 2(e_1e_3-e_0e_2) & 2(e_2e_3+e_0e_1) & e_0^2-e_1^2-e_2^2+e_3^2 \end{matrix}\right ]\!, \end{equation}

where $\boldsymbol e=[e_0,e_1, e_2, e_3]$ are unit quaternions. The corresponding evolution of $\boldsymbol{e}$ can be expressed by

(2.5) \begin{equation} \frac {\mathrm{d}{\boldsymbol{e}}}{\mathrm{d}t^*}=\unicode{x1D642}^T\boldsymbol{\cdot }\boldsymbol{\omega}^{* \prime }. \end{equation}

Here,

(2.6) \begin{equation} \unicode{x1D642}=\left [\begin{matrix} -e_1 & e_0 & -e_3 & e_2 \\ -e_2 & -e_3 & e_0 & e_1 \\ -e_3 & e_2 & -e_1 & e_0 \end{matrix}\right ]\!. \end{equation}

In this study, the resultant torque $\boldsymbol{T'}$ consists of Jeffery torque $T_{\!J}$ and fluid-inertial torque $T_I$ . The Jeffery torque in any frame satisfies the following form (Brenner Reference Brenner1963):

(2.7) \begin{equation} \boldsymbol T_{\!J}=\unicode{x1D646}_R\boldsymbol{\cdot }(\boldsymbol \varOmega ^*-\boldsymbol \omega ^*)+\unicode{x1D643} : \unicode{x1D64E}^{\;*}, \end{equation}

where $\unicode{x1D646}_R$ and $\unicode{x1D643}$ are the resistant coefficient tensors of rotation and shear, respectively, only related to the particle shape and the viscosity of the fluid. Note that $\unicode{x1D643}$ is a third-order tensor. In the particle coordinate,

(2.8) \begin{align} T^{\prime}_{J x} &=\frac {16\pi \mu b^3}{3\chi _\perp }\big(\varOmega ^{*\prime }_x-\omega ^{*\prime }_x\big),\nonumber\\ T^{\prime}_{J y} &=\frac {16\pi \mu b^3}{3(\lambda ^2\chi _\parallel +\chi _\perp )}\big[(1-\lambda ^2)S^{*\prime }_{xz}+(1+\lambda ^2)\big(\varOmega ^{*\prime }_y-\omega ^{*\prime }_y\big)\big],\nonumber\\T^{\prime}_{J z} &=\frac {16\pi \mu b^3}{3(\lambda ^2\chi _\parallel +\chi _\perp )}\big[(\lambda ^2-1)S^{*\prime }_{xy}+(\lambda ^2+1)\big(\varOmega ^{*\prime }_z-\omega ^{*\prime }_z\big)\big]. \end{align}

Here, $\chi _\perp$ and $\chi _\parallel$ are the integral shape parameters as listed in table 1. Note that the subscript $\parallel$ means the coefficients along the symmetry axis while the subscript $\perp$ represents the coefficients normal to the symmetry axis. The fluid-inertial torque in the Eulerian frame is expressed as

(2.9) \begin{equation} \boldsymbol T_I =-\rho _{\!f}r_{\!p}^3 F_\lambda (\boldsymbol u_s \boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol u_s\times \boldsymbol n), \end{equation}

where $r_{\!p}$ represents the semilength of the major axis of the spheroid, i.e. $r_{\!p}=a$ for probate spheroids and $r_{\!p}=b$ for oblate spheroids. In the particle coordinate,

(2.10) \begin{equation} \boldsymbol T^{\prime}_I =-\rho _{\!f}r_{\!p}^3 F_\lambda u^{\prime}_{s 1}\left (\begin{matrix} 0\\ u^{\prime}_{s 3} \\[3pt] -u^{\prime}_{s 2} \end{matrix}\right )\!, \end{equation}

where $u^{\prime}_{s 1}$ , $u^{\prime}_{s 2}$ and $u^{\prime}_{s 3}$ are the three components of slip velocity in the particle coordinate, respectively. Note that $F_\lambda$ is a function of fluid-inertial torque depending on particle shape (Dabade et al. Reference Dabade, Marath and Subramanian2015) (see Appendix A).

Table 1.

Shape parameters of a spheroid.

In the shear flow, we use the shear rate $\kappa$ to non-dimensionlize the angular dynamics of spheroids, i.e. particle rotation $\boldsymbol \omega =\boldsymbol \omega ^*\!/\!\kappa$ , fluid rotation $\boldsymbol \varOmega =\boldsymbol \varOmega ^*\!/\!\kappa$ , strain-rate $\unicode{x1D64E}=\unicode{x1D64E}^{\;*}\!/\!\kappa$ , fluid rotation tensor $\unicode{x1D64A}=\unicode{x1D64A}^*\!/\!\kappa$ and time $t=t^*\kappa$ . The non-dimensional governed equations of angular dynamics of inertial spheroidal particles in the Eulerian frame are expressed as follows:

(2.11) \begin{align} St_r \unicode{x1D63E}\boldsymbol{\cdot }\frac {\mathrm{d}{\boldsymbol \omega }}{\mathrm{d}t}+\boldsymbol{n}\times (St_r C_\perp \varLambda \boldsymbol \omega \boldsymbol \omega )\boldsymbol{\cdot }\boldsymbol n&=\boldsymbol \varOmega -\boldsymbol \omega +\boldsymbol n\times \varLambda \unicode{x1D64E} \boldsymbol{\cdot }\boldsymbol n+\boldsymbol n\times M\hat {\boldsymbol u}_s\hat {\boldsymbol u}_s\boldsymbol{\cdot }\boldsymbol n,\nonumber \\ \frac {\mathrm{d}{\boldsymbol n}}{\mathrm{d}t}&=\boldsymbol \omega \times \boldsymbol n. \end{align}

Here, $\boldsymbol n$ is the direction of the symmetry axis of the spheroid, and the left-hand side equation refers to the particle inertia term. In the right-hand side, $\boldsymbol \varOmega -\boldsymbol \omega +\boldsymbol n\times \varLambda \unicode{x1D64E} \boldsymbol{\cdot }\boldsymbol n$ is the term of Jeffery torque, while $\boldsymbol n\times M\hat {\boldsymbol u}_s\hat {\boldsymbol u}_s\boldsymbol{\cdot }\boldsymbol n$ is the term of fluid-inertial torque. Moreover, $\hat {\boldsymbol u}_s$ is a unit vector of slip velocity, which directs to the streamwise direction. Meanwhile, $\unicode{x1D63E}$ , $C_\perp$ and $\varLambda$ are coefficients related to particle shape, i.e.

(2.12) \begin{align} \unicode{x1D63E} & =C_\parallel \boldsymbol n\boldsymbol n+C_\perp (\unicode{x1D640}-\boldsymbol n \boldsymbol n), \nonumber \\ C_\perp & =\frac {k_\perp ^2k_\parallel (k_\perp ^2\chi _\perp +k_\parallel ^2\chi _\parallel )}{20}, \nonumber \\ C_\parallel & =\frac {2k_\perp ^4k_\parallel \chi _\perp }{20}, \nonumber \\ \varLambda & =\frac {k_\parallel ^2-k_\perp ^2}{k_\parallel ^2+k_\perp ^2}\equiv \frac {\lambda ^2-1}{\lambda ^2+1}. \end{align}

Here, $k_\perp$ , $k_\parallel$ are shape parameters as listed in table 1, which are defined as $k_\perp =b/r_{\!p}$ and $k_\parallel =a/r_{\!p}$ , respectively. Note that $r_{\!p}$ is the semilength of the major axis of spheroids, depending on whether the shape is prolate or oblate. Meanwhile, $M$ reflects the direction and the strength of fluid-inertial torque, which is linked to particle shape, slip velocity and shear rate, i.e.

(2.13) \begin{equation} M =\frac {3 F_\lambda \big(k_{\perp }^2\chi _\perp +k_{\parallel }^2\chi _\parallel \big) }{16\pi \big(k_{\perp }^2+k_{\parallel }^2\big)}\frac {u_s^2}{\nu \kappa }. \end{equation}

Here, the ratio $u_s^2/\nu \kappa$ , denoted as $\zeta ^2$ , reflects the relative importance of the fluid-inertial torque compared with the Jeffery torque (Sheikh et al. Reference Sheikh, Gustavsson, Lopez, Lévêque, Mehlig, Pumir and Naso2020), i.e. $\zeta ^2\sim |\boldsymbol T_I|/|\boldsymbol T_{\!J}|$ . A larger $\zeta$ indicates a stronger effect of fluid-inertial torque because $|M|\propto \zeta ^2$ . Moreover, the definition $\zeta ^2=u_s^2/\nu \kappa$ can also be reformulated as $\zeta ^2 ={\textit{Re}}_{\!p}^2 /{\textit{Re}}_\kappa$ , leading to $\zeta ={\textit{Re}}_{\!p}/\sqrt {{\textit{Re}}_\kappa }$ .

It is worth clarifying the regime of validity of the PP model described in (2.11). Although the last term on the right-hand side of (2.11) includes the slip-induced torque contribution, torque contributions due to the fluid-inertial effects associated with shear alone have been neglected. While in the general case, the first such shear-induced inertial contribution arises at $O(\sqrt {{\textit{Re}}_\kappa })$ (Candelier, Mehlig & Magnaudet Reference Candelier, Mehlig and Magnaudet2019), for the axisymmetric spheroidal particles of interest, this contribution vanishes, and the first inertial correction to the Jeffery torque is $O({\textit{Re}}_\kappa )$ . The relative magnitudes of the slip- and shear-induced fluid-inertial contributions are determined by the ratio $\zeta ^2/{\textit{Re}}_\kappa = ({\textit{Re}}_{\!p}/{\textit{Re}}_\kappa )^2$ . It is worth noting that, in the limit that both these contributions remain small compared with the Jeffery torque, they may be linearly superimposed regardless of their relative magnitudes, that is to say, regardless of the ratio ${\textit{Re}}_{\!p}/{\textit{Re}}_\kappa$ . This superposition may be justified based on symmetry arguments that rely on both these contributions being inner-region contributions, as a result of which the cross-coupling terms are identically zero at leading order. Since we neglect the $O({\textit{Re}}_\kappa )$ contribution in (2.11), the current PP model is restricted to $\zeta ^2/{\textit{Re}}_\kappa \gg 1$ , or ${\textit{Re}}_\kappa \ll {\textit{Re}}_{\!p}$ . In the opposite limit where the slip-induced fluid-inertial torque becomes comparable to the leading-order Jeffery torque ( $\zeta \gt O(1)$ ), the nature of the shear-inertial torque contribution is not clear. This is because while expressions for this contribution have been obtained by Einarsson et al. (Reference Einarsson, Candelier, Lundell, Angilella and Mehlig2015b ) and Dabade et al. (Reference Dabade, Marath and Subramanian2016), these are for the specific case where the spheroid, at leading order, rotates along Jeffery orbits, on account of it being neutrally buoyant. This rotation then determines the torque contribution of the unsteady acceleration term. When the spheroid rotation is modified due to the effects of the slip-induced fluid-inertial torque, the unsteady acceleration contribution will also be changed, and the resulting expression for the shear-induced inertial contribution has not yet been calculated. Nevertheless, we perform PR simulations, as detailed in § 4, to assess the validity of the present PP model in (2.11).

For stability analysis, (2.11) is more convenient and concise than (2.2) where the transformation between particle coordinate and laboratory reference frame is frequent. To investigate weak inertial spheroidal particles, (2.11) can be degenerated to a weakly inertial form by a power series of $\boldsymbol \omega$ with a small parameter $St_r$ , namely $\boldsymbol \omega =\boldsymbol \omega _0+St_r \boldsymbol \omega _1+O(St_r^2)$ . Therefore, the weakly inertial form of the orientational dynamics of a spheroid is

\begin{align*} \frac {\mathrm{d}{\boldsymbol n}}{\mathrm{d}t}=\,&\unicode{x1D64A}\boldsymbol{\cdot }\boldsymbol n+\unicode{x1D63D}_s\boldsymbol{\cdot }\boldsymbol n-(\unicode{x1D63D}_s:\boldsymbol n \boldsymbol n)\boldsymbol n+\widetilde {St_r} \lbrace (\unicode{x1D63D}_s:\boldsymbol n\boldsymbol n) \unicode{x1D64A} \boldsymbol{\cdot }\boldsymbol n\nonumber \\ &+[-\unicode{x1D63D}_s\boldsymbol{\cdot }\unicode{x1D64A}-\unicode{x1D63D}_s \boldsymbol{\cdot }\unicode{x1D63D}_s+2(\unicode{x1D63D}_s:\boldsymbol n \boldsymbol n)\unicode{x1D63D}_s]\boldsymbol{\cdot }\boldsymbol n\nonumber \\ &-\boldsymbol n\boldsymbol{\cdot }[-\unicode{x1D63D}_s\boldsymbol{\cdot }\unicode{x1D64A}-\unicode{x1D63D}_s \boldsymbol{\cdot }\unicode{x1D63D}_s+2(\unicode{x1D63D}_s:\boldsymbol n \boldsymbol n)\unicode{x1D63D}_s] \boldsymbol{\cdot }\boldsymbol n\boldsymbol n \end{align*}

(2.14) \begin{align} & -\varLambda [\unicode{x1D64A}\boldsymbol{\cdot }\unicode{x1D64A}+\unicode{x1D64A} \boldsymbol{\cdot }\unicode{x1D63D}_s-(\unicode{x1D63D}_s:\boldsymbol n\boldsymbol n)\unicode{x1D64A}]\boldsymbol{\cdot }\boldsymbol n+\varLambda \boldsymbol n\boldsymbol{\cdot }(\unicode{x1D64A}\boldsymbol{\cdot }\unicode{x1D64A}+\unicode{x1D64A}\boldsymbol{\cdot }\unicode{x1D63D}_s)\boldsymbol{\cdot }\boldsymbol n\boldsymbol n\rbrace \nonumber \\ &+O\big(\widetilde {St_r}^2\big), \end{align}

where $\widetilde {St_r}\equiv St_r C_\perp$ , and $\unicode{x1D63D}_s$ is a symmetrical tensor and expressed as $\unicode{x1D63D}_s=\varLambda \unicode{x1D64E} +M\hat {\boldsymbol u}_s\hat {\boldsymbol u}_s$ . The details of the deviation for (2.14) can be found in Appendix B. In addition to the leading-order term, we also retain $\widetilde {St_r}$ -order terms, which primarily reflect the particle inertial effect on the spheroid orientation. The comparison between the results of (2.11) and (2.14) at very small $St_r$ is shown in Appendix B, reflecting a validity of this weakly inertial form of spheroid orientation.

In this study, we define two Stokes numbers: $St_r$ and $\widetilde {St_r}$ . The former, $St_r$ , follows the conventional definition widely used in previous studies, while the latter is defined as $\widetilde {St_r}\equiv St_r C_\perp$ . To illustrate the relationship between these two Stokes numbers, we plot the factor $1/C_\perp$ , which is the scaling factor between $St_r$ and $\widetilde {St_r}$ , as a function of the aspect ratio $\lambda$ in figure 3. According to figure 3, $1/C_\perp$ increases exponentially as the aspect ratio $\lambda$ deviates from unity. This indicates that the amplitude of $St_r$ can be large even for a small value of $\widetilde {St_r}$ for particles with a shape deviating far away from a sphere. Nevertheless, $\widetilde {St_r}$ serves only as an approximate measure of particle inertia; thus, a large $\widetilde {St_r}$ indicates dominant particle-inertia effects, and vice versa. Hereinafter, we use $\widetilde {St_r}$ instead of $St_r C_\perp$ for the sake of the simplicity of the expression.

Figure 3.

Relationship between $St_r/\widetilde {St_r}=1/C_\perp$ and aspect ratio $\lambda$ .

2.3. Equilibrium rotation modes of spheroids in shear flow

There are three equilibrium rotation modes of spheroids, i.e. logrolling, tumbling and aligning, as shown in figure 4. The transitional rotation modes are not considered in the current study. The logrolling mode signifies spheroids rotating about the symmetry axis, which orients to the vortical axis. The tumbling mode represents a spheroid rotating about its equatorial axis with its symmetry axis in the flow-gradient plane. The aligning mode means that the spheroid cannot rotate, but its symmetry axis keeps an angle relative to the streamwise direction.

Figure 4.

Rotation modes of spheroids: (a) logrolling; (b) tumbling; (c) aligning.

For the logrolling, which is the simplest mode, the equilibrium state is $\boldsymbol n^e=[0,0,\pm 1]^T$ .

For the steady aligning mode, the corresponding orientation $\boldsymbol n^e$ is derived from (2.11) by setting $\boldsymbol \omega ^e=0$ and $({\mathrm{d}{\boldsymbol \omega ^e}}/{\mathrm{d}t})=0$ . Then (2.11) reduces to

(2.15) \begin{equation} \boldsymbol \varOmega =\boldsymbol n^e \times (\varLambda \unicode{x1D64E}+M\hat {\boldsymbol u}\hat {\boldsymbol u})\boldsymbol{\cdot }\boldsymbol n^e. \end{equation}

In the current shear flow, $\boldsymbol u=\kappa \!\boldsymbol y+\boldsymbol u_s$ , the aligning mode corresponds to the particle orientation vector $\boldsymbol n^e$ lying in the $x$ – $y$ plane. Under this condition, only the third component of the flow rotation $\boldsymbol \varOmega$ remains,

(2.16) \begin{equation} -\!\frac {1}{2}=-\frac {\varLambda }{2}{n^e_x}^2+Mn^e_x n^e_y+\frac {\varLambda }{2}{n^e_y}^2. \end{equation}

Meanwhile, the unit orientation vector $\boldsymbol n^e$ lying in the $x$ – $y$ plane satisfies the constraint ${n^e_x}^2+{n^e_y}^2=1$ . By solving above equations, two pairs of solutions for $\boldsymbol n^e$ are obtained as follows:

(2.17) \begin{align} \left \{ \begin{aligned} {n^e_x}^2_{(1)}&=\dfrac {\varLambda ^2+M^2+\varLambda +\sqrt {\varLambda ^2M^2+M^4-M^2}}{2\varLambda ^2+2M^2} \\[4pt] {n^e_y}^2_{(1)}&=1-{n^e_x}^2_{(1)}=\dfrac {\varLambda ^2+M^2-\varLambda -\sqrt {\varLambda ^2M^2+M^4-M^2}}{2\varLambda ^2+2M^2}\nonumber \\ \end{aligned}\right . , \nonumber \\ \left \{ \begin{aligned} {n^e_x}^2_{(2)}&=\dfrac {\varLambda ^2+M^2+\varLambda -\sqrt {\varLambda ^2M^2+M^4-M^2}}{2\varLambda ^2+2M^2} \\[4pt]{n^e_y}^2_{(2)}&=1-{n^e_x}^2_{(2)}=\dfrac {\varLambda ^2+M^2-\varLambda +\sqrt {\varLambda ^2M^2+M^4-M^2}}{2\varLambda ^2+2M^2} \end{aligned}\right . . \end{align}

To be concise, the component $n_x^e$ is denoted as $p$ , leading to the following two solutions of $p$ satisfying

(2.18) \begin{align}& p_{(1)}^2=\frac {\varLambda ^2+M^2+\varLambda +\sqrt {\varLambda ^2M^2+M^4-M^2}}{2\varLambda ^2+2M^2},\nonumber \\ &p_{(2)}^2=\frac {\varLambda ^2+M^2+\varLambda -\sqrt {\varLambda ^2M^2+M^4-M^2}}{2\varLambda ^2+2M^2}. \end{align}

Finally, the equilibrium orientation $n^e$ of aligning mode satisfies

(2.19) \begin{equation} \boldsymbol n^e=\left[p, \frac {2\varLambda p^2-\varLambda -1}{2Mp},0\right]^T. \end{equation}

For each $p_{(1)}^2$ or $p_{(2)}^2$ , there is a pair of positive and negative values, representing two opposite vectors, i.e. $\boldsymbol n_+^e=-\boldsymbol n_{-}^e$ . Note that the subscripts ‘ $+$ ’ and ‘ $-$ ’ represent the directions computed by positive and negative $p$ , respectively. Because of the symmetry of the spheroid, the behaviours of the spheroid corresponding to these two vectors are identical.

Regarding tumbling, which is the most complicated mode, the equilibrium state is a periodic function with time in the flow-gradient plane, denoted as $\boldsymbol n^e(t)$ . According to Cui et al. (Reference Cui, Liu, Qiu and Zhao2024), when $\widetilde {St_r}$ is small, the non-dimensional period of tumbling rotation is $T_{p}^H=4\pi /\sqrt {1-\varLambda ^2-M^2}$ , while close to $T_{p}^L=4 \pi$ for large $\widetilde {St_r}$ . Moreover, around $\widetilde {St_r}\approx 1$ , the period of tumbling rotation transients from $T_{p}^H$ to $T_{p}^L$ . Meanwhile, when Stokes number $\widetilde {St_r}$ is non-zero, the periodic motions of the tumbling are only acquired numerically by solving a quasi-two-dimensional equation (Cui et al. Reference Cui, Liu, Qiu and Zhao2024), as expressed by

(2.20) \begin{equation} \widetilde {St_r}\frac {\mathrm{d}{\omega _z}}{\mathrm{d}t}=\frac {1}{2}\sqrt {\varLambda ^2+M^2}\cos (2\varphi +2\varphi _0)+\left(-\frac {1}{2}-\omega _z\right)\!, \end{equation}

where $2\varphi _0=\arctan (M/\varLambda )$ , $\omega _z$ and $\varphi$ are the angular velocity and displacement of a spheroid rotating in the flow-gradient plane, respectively. The details can be found in Cui et al. (Reference Cui, Liu, Qiu and Zhao2024). For finite $\widetilde {St_r}$ , there are no analytical solutions for the tumbling mode. Therefore, we primarily use numerical solutions of the tumbling mode as equilibrium states. Note that the tumbling mode always exists when $\varLambda ^2 + M^2 \lt 1$ . When $\varLambda ^2 + M^2 \geqslant 1$ , the existence of the tumbling mode is no longer guaranteed. As shown in figure 5(a), there exists an upper critical value $A_c$ of $\sqrt {\varLambda ^2 + M^2}$ that depends on $\widetilde {St_r}$ . According to figure 5(b), when the parameter $\zeta$ exceeds the upper critical threshold, indicated by the green lines, only the aligning mode exists, regardless of the initial conditions. Conversely, when $\zeta$ falls below the magenta lines corresponding to $\varLambda ^2 + M^2 = 1$ , only the tumbling mode remains. Between the green and magenta lines, there is a transition region where both aligning and tumbling modes coexist. In this region, the final equilibrium state is sensitive to the initial conditions, including both the orientation and angular velocity of the spheroid.

Figure 5.

(a) Upper critical values $A_c=\sqrt {\varLambda ^2+M^2}$ versus Stokes number $\widetilde {St_r}$ . (b) Regimes of equilibrium tumbling and aligning modes at $\widetilde {St_r}=5$ , where green lines represent $\varLambda ^2+M^2=3.18^2$ and magenta lines stand for $\varLambda ^2+M^2=1$ .

3.1. The case when $\widetilde {St_r}\ll 1$

When $\widetilde {St_r} \ll 1$ , (2.14) reflects the rotation motions of spheroids. In this section, (2.14) has the form of $({\mathrm{d}{\boldsymbol n}}/{\mathrm{d}t})= \mathcal{N}(\boldsymbol n)$ and the corresponding perturbation equation is $({\mathrm{d}{\delta \boldsymbol n}}/{\mathrm{d}t})=({\partial }/{\partial \boldsymbol n}) \mathcal{N}(\boldsymbol n^e)\boldsymbol{\cdot }\delta \boldsymbol n$ , where $\boldsymbol n^e$ represents an orientation of equilibrium state and $({\partial }/{\partial \boldsymbol n}) \mathcal{N}(\boldsymbol n^e)$ refers to a stability matrix, denoted as $\unicode{x1D64C}(\boldsymbol n^e)$ . Due to the constraint $n_x^2+n_y^2+n_z^2=1$ , the three-dimensional orientational dynamics of the particle are effectively reduced into a two-dimensional system in $(n_x, n_y)$ subspace by substituting $n_z=\pm \sqrt {1-n_x^2-n_y^2}$ into the orientational dynamics equations, i.e.

(3.1) \begin{equation} \frac {\mathrm{d}{}}{\mathrm{d}t} \left [\begin{matrix} n_x\\n_y\\ \end{matrix}\right ]=\mathcal{N}(n_x,n_y), \end{equation}

For logrolling and aligning states, $\boldsymbol n^e$ is independent of time, resulting in a constant stability matrix $\unicode{x1D64C}$ . Therefore, we can use a stability theory of fixed points to analyse logrolling and aligning modes. However, the tumbling modes are periodic with time, which should be analysed by Floquet theory (Joseph Reference Joseph1976; Chicone Reference Chicone2006).

3.1.1. Logrolling mode

For the logrolling mode, the orientation of the spheroid is $\boldsymbol n^e=[0,0,\pm 1]^T$ . The perturbation equation of logrolling mode is

(3.2) \begin{equation} \frac {\mathrm{d}{}}{\mathrm{d}t}\left [\begin{matrix} \delta n_x\\\delta n_y\\ \end{matrix}\right ]=\left [\begin{matrix} \dfrac {\left(-\varLambda ^2-2M^2+\varLambda \right)\widetilde {St_r}}{2}+M & -\dfrac {(1+\varLambda )\left(M\widetilde {St_r}-1\right)}{2}\\[8pt] -\dfrac {1}{2}+\dfrac {\varLambda }{2}&0\\ \end{matrix}\right ]\left [\begin{matrix} \delta n_x\\ \delta n_y\\ \end{matrix}\right ]\!. \end{equation}

Consequently, the eigenvalues of the corresponding stability matrix $\unicode{x1D64C}$ are shown as follows:

(3.3) \begin{align} b_{1,2}&=\frac {M}{2}+\widetilde {St_r} \frac {\left(-\varLambda ^2-2M^2+\varLambda \right)}{4}\nonumber\\ &\quad \pm \frac {\sqrt {(2M^2+\varLambda ^2-\varLambda )^2\widetilde {St_r}^2-4(2M^2+2\varLambda ^2-\varLambda -1)M \widetilde {St_r}+4M^2+4\varLambda ^2-4}}{4}\nonumber\\ &\quad +O\left(\widetilde {St_r}^2\right)\!, \end{align}

where the signs of $b_{1,2}$ determine the stability of the logrolling mode. When $\widetilde {St_r}=0$ , (3.3) degenerates into

(3.4) \begin{align} &b_{1,2}=\frac {M}{2}\pm \frac {\sqrt {M^2+\varLambda ^2-1}}{2}, \end{align}

which is consistent with the results in Cui et al. (Reference Cui, Liu, Qiu and Zhao2024). For prolate spheroids, when the fluid-inertial torque is considered, $M$ is negative, and the maximum real part of $b_{1,2}$ is also negative, reflecting that the logrolling of the prolate spheroid is stable. For oblate spheroids, $M$ is positive, and the maximum real part of $b_{1,2}$ is positive, indicating that the logrolling of an oblate spheroid is unstable. Note that $\varLambda ^2$ is linked to shape and is always below unity. When particle inertia is taken into consideration, the stability map changes with increasing $\widetilde {St_r}$ and varying $\zeta$ , as depicted in figure 6. For $\widetilde {St_r}=0.1$ (see figure 6 a), the logrolling mode for prolate spheroids is unstable but it is stable for oblate spheroids when $\zeta$ is small. It is obvious that the critical $\zeta _c$ for neutral stability of thin oblate ones is nearly independent of $\lambda$ and larger than that of prolate ones, which is attributed to the nearly constant and large $F_\lambda$ for oblate spheroids. Meanwhile, the critical $\zeta _c$ decreases as the aspect ratio $\lambda$ increases for prolate spheroids. This indicates that the stability of the logrolling of the prolate spheroids is easily affected by the fluid-inertial torque from stable state to unstable state, but for the oblate spheroids, $\zeta _c$ is relatively higher. At the same time, for the fixed aspect ratios, e.g. $\lambda =0.1$ and $10$ as shown in figure 6(b), the critical $\zeta _c$ and $\widetilde {St_r}$ exhibit a power-law relationship $\widetilde {St_r}\sim \zeta _c^2$ , indicating that as particle inertia increases, a larger fluid-inertial torque is required to overcome the effect of the particle inertia.

Figure 6.

Stability of the logrolling of spheroidal particle computed by (3.3). (a) The maximum real parts of eigenvalues $b_{1,2}$ with different $\zeta$ and $\lambda$ at $\widetilde {St_r}=0.1$ , and (b) boundaries between stable and unstable regions of spheroids versus particle inertia $\widetilde {St_r}$ within $0\sim 0.5$ for $\lambda =0.1$ (solid line) and $\lambda =10$ (dashed line). In panel (a), solid lines represent zero value, and magenta dashed–dotted lines represent $\varLambda ^2+M^2=1$ . Note that when the maximum real part of the eigenvalue is negative, the logrolling mode is stable, and vice versa. Here ‘US’ and ‘S’ represent unstable and stable states, respectively.

3.1.2. Aligning mode

The aligning mode only exists when $A^2=M^2+\varLambda ^2\geqslant 1$ , and there are two pairs of orientations referring to the equilibrium states. The first pair satisfies $p_{(1)}^2= ({\varLambda ^2+M^2+\varLambda +\sqrt {\varLambda ^2M^2+M^4-M^2}})/({2\varLambda ^2+2M^2})$ , and the corresponding eigenvalues of the stability matrix $\unicode{x1D64C}$ are

(3.5) \begin{align} b_1= & -\frac {M^2-|M|M +2\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{2M} \nonumber \\ & -\frac {\widetilde {St_r}}{4} \left [\frac {3\varLambda ^2M + \varLambda ^2|M|+ 4M^3 - 3M - |M|}{M}\right . \nonumber \\ &\left .+\frac { (2M-|M|)\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{M}\right ]+O\big(\widetilde {St_r}^2\big), \end{align}

(3.6) \begin{align} b_2= & -\frac {M^2+|M|M +2\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{2M} \nonumber \\ & -\frac {\widetilde {St_r}}{4}\left [\frac {3\varLambda ^2M-\varLambda ^2|M|+ 4M^3 - 3M + |M|}{M}\right . \nonumber \\ &\left . +\frac { (2M+|M|)\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{M}\right ]+O\big(\widetilde {St_r}^2\big). \end{align}

For the second one, $p_{(2)}^2=({\varLambda ^2+M^2+\varLambda -\sqrt {\varLambda ^2M^2+M^4-M^2}})/({2\varLambda ^2+2M^2})$ , and the corresponding eigenvalues of the stability matrix $\unicode{x1D64C}$ are

(3.7) \begin{align} b_1= & -\frac {M^2-|M|M -2\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{2M} \nonumber \\ & -\frac {\widetilde {St_r}}{4} \left [\frac {3\varLambda ^2M + \varLambda ^2|M|+ 4M^3 - 3M - |M|}{M}\right . \nonumber \\ &\left .-\frac { (2M-|M|)\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{M}\right ]+O\big(\widetilde {St_r}^2\big), \end{align}

(3.8) \begin{align} b_2= & -\frac {M^2+|M|M -2\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{2M} \nonumber \\ & -\frac {\widetilde {St_r}}{4}\left [\frac {3\varLambda ^2M-\varLambda ^2|M|+ 4M^3 - 3M + |M|}{M}\right . \nonumber \\ &\left . -\frac { (2M+|M|)\sqrt {\varLambda ^2M^2 + M^4 - M^2}}{M}\right ]+O\big(\widetilde {St_r}^2\big). \end{align}

The stability of the aligning mode is shown in the regions of $\varLambda ^2+M^2\geqslant 1$ in figure 7(a) at $\widetilde {St_r}=0.1$ . A stable aligning mode for oblate spheroids and an unstable align mode for prolate ones at small Stokes number with a pronounced slip velocity are observed, which explains why we can only observe the aligning mode for oblate spheroids instead of prolate ones at small $\widetilde {St_r}$ in our previous PP simulation results (Cui et al. Reference Cui, Liu, Qiu and Zhao2024).

Figure 7.

Stability of the aligning (in $\varLambda ^2+M^2\geqslant 1$ ) and the tumbling (in $\varLambda ^2+M^2\lt 1$ ) of spheroidal particle computed by (3.5)–(3.7) and (3.14). (a) The maximum real parts of eigenvalues in (3.5–3.7) and the values of Floquet exponent in (3.14) with different $\zeta$ and $\lambda$ at $\widetilde {St_r}=0.1$ , and (b) boundaries between stable and unstable regions of spheroids versus particle inertia $\widetilde {St_r}$ within $0\sim 0.5$ for $\lambda =0.1$ (solid line) and $\lambda =10$ (dashed line). In panel (a), solid lines represent zero values, and magenta dashed–dotted lines represent $\varLambda ^2+M^2=1$ . Note that when the contour values shown in (a) are negative, the corresponding rotation mode is stable, and vice versa. Here ‘US’ and ‘S’ represent unstable and stable states, respectively.

3.1.3. Tumbling mode

For the tumbling mode ( $M^2+\varLambda ^2\lt 1$ ), the equilibrium motion of $\boldsymbol n^e$ is a periodic function with the time. Inspired by the study of Einarsson et al. (Reference Einarsson, Angilella and Mehlig2014), analytical solutions of the eigenvalues of the stability matrix are conducted via Floquet theory with $\widetilde {St_r} \ll 1$ . By substituting unit vector $\boldsymbol n=[\sin \vartheta \cos \varphi ,\sin \vartheta \sin \varphi ,\cos \vartheta ]$ into (2.14), the orientational dynamics of the spheroid is expressed by

(3.9) \begin{align} \frac {\mathrm{d}{\varphi }}{\mathrm{d}t}= & f(\varphi ,\vartheta ) \nonumber \\ = & \frac {1}{2}(-1+\varLambda \cos 2\varphi -M\sin 2\varphi )+\frac {\widetilde {St_r}}{8}\lbrace [2M^2\sin 2\varphi \nonumber \\ & -M(\varLambda -1)\cos 2\varphi +M^2\sin 4\varphi -M(2\varLambda \cos 4\varphi -\varLambda -1)]\cos 2\vartheta \nonumber \\ & +2M^2\sin 2\varphi -M(\varLambda +3)\cos 2\varphi -M^2\sin 4\varphi \nonumber \\ & +M(2\varLambda \cos 4\varphi +\varLambda +1)+2\varLambda ^2\sin 4\varphi \sin ^2\vartheta \nonumber \\ & +\sin 2\varphi [\varLambda (\varLambda +1)\cos 2\vartheta +\varLambda ^2-3\varLambda ]\rbrace +O\big(\widetilde {St_r}^2\big), \end{align}

(3.10) \begin{align} \frac {\mathrm{d}{\vartheta }}{\mathrm{d}t}= &\, g(\varphi ,\vartheta ) \nonumber \\ = & \sin \vartheta \cos \vartheta \cos \varphi (M\cos \varphi +\varLambda \sin \varphi ) \nonumber \\ & + \widetilde {St_r}\sin \vartheta \cos \vartheta \cos ^2\varphi \left \lbrace 2 M^2\sin ^2\vartheta \cos ^2\varphi \vphantom{\frac {\varLambda }{2}}\right . \nonumber \\ & +4M\varLambda \sin \varphi \cos \varphi \sin ^2\vartheta -M^2-\frac {M(\varLambda +1)}{2}\tan \varphi \nonumber \\ &\left . +\,2\varLambda ^2\sin ^2\vartheta \sin ^2\varphi -\frac {\varLambda ^2}{2}+\frac {\varLambda }{2}\right \rbrace +O\big(\widetilde {St_r}^2\big). \end{align}

Here, $\varphi$ denotes the relative angle between the x-axis and the projection of $\boldsymbol n$ in the x–y plane, while $\vartheta$ represents the angle of $\boldsymbol n$ with respect to the vortical axis, i.e. z-axis. According to (3.9), the evolution of $\varphi$ and $\vartheta$ are decoupled only if $\widetilde {St_r}=0$ . However, the presence of particle inertia makes $\varphi$ and $\vartheta$ coupled with each other, leading to a strong nonlinearity. When a spheroid is tumbling in the x–y plane, its equilibrium $\vartheta _0=\pi /2$ , while $\varphi _0(t)$ is a periodic function with a period $T_{\!p}=4\pi /\sqrt {1-\varLambda ^2-M^2}$ when $\widetilde {St_r}\ll 1$ (Cui et al. Reference Cui, Liu, Qiu and Zhao2024). By considering small perturbations on the equilibrium, i.e. $\vartheta =\vartheta _0+\delta \vartheta$ and $\varphi (t)=\varphi _0(t)+\delta \varphi (t)$ , the corresponding linear perturbation equation of $\delta \vartheta$ is

(3.11) \begin{equation} \frac {\mathrm{d}{}}{\mathrm{d}t}\delta \vartheta =\delta \vartheta \frac {\partial g}{\partial \vartheta }\Bigg|_{(\varphi _0(t),\vartheta _0)}+\delta \varphi \frac {\partial g}{\partial \varphi }\Bigg|_{(\varphi _0(t),\vartheta _0)}. \end{equation}

When $\vartheta _0=\pi /2$ , $g$ is independent of $\varphi$ , namely $({\partial g}/{\partial \varphi })=0$ . Hence, the above equation is simplified to

(3.12) \begin{equation} \frac {\mathrm{d}{}}{\mathrm{d}t}\ln \delta \vartheta =\left . \frac {\partial g}{\partial \vartheta }\right|_{(\varphi _0(t),\vartheta _0)}. \end{equation}

According to the Floquet theory, we integrate (3.12) along the trajectory ( $\varphi _0, \vartheta _0$ ) with a period $T_{\!p}$ and average the Floquet exponent by $T_{\!p}$ , i.e.

(3.13) \begin{equation} \gamma _{\pi /2}=\frac {1}{T_{\!p}}\int _0^{T_{\!p}}{\rm d}t\frac {\partial g}{\partial \vartheta }\Bigg|_{(\varphi _0(t),\vartheta _0)}=-\frac {1}{T_{\!p}}\int _0^{2\pi }\frac {{\rm d} \varphi }{\dot \varphi }\frac {\partial g}{\partial \vartheta }\Bigg|_{(\varphi _0(t),\vartheta _0)}. \end{equation}

Here, $\dot {\varphi }={\mathrm{d}{\varphi }}/{\mathrm{d}t}$ . After integrating above equation and retaining the first-order term of $\widetilde {St_r}$ , the Floquet exponent $\gamma _{\pi /2}$ of the tumbling mode is expressed as

(3.14) \begin{equation} \gamma _{\pi /2}=-\frac {M}{2}+\frac {\widetilde {St_r}}{4}\left [M^2+\varLambda ^2+(1-\varLambda )\sqrt {1-M^2-\varLambda ^2}-1\right ]+O\big(\widetilde {St_r}^2\big). \end{equation}

When $\gamma _{\pi /2}\lt 0$ , the tumbling mode is stable and vice versa. When $M=0$ , the simplified form of (3.14) is consistent with the results reported in Einarsson et al. (Reference Einarsson, Angilella and Mehlig2014). The corresponding results of (3.14) are shown in figure 7. When the fluid inertia effect is weaker than the particle inertia effect, the tumbling of prolate spheroids is stable but not for oblate ones. Contrarily, when the fluid inertia effect is greater than the particle inertia effect, the tumbling of prolate spheroids becomes unstable, but that of the oblate spheroids becomes stable. These results are consistent with our previous simulations (Cui et al. Reference Cui, Liu, Qiu and Zhao2024). Besides, for fixed aspect ratios, such as $\lambda =0.1$ and $10$ , the critical $\zeta _c$ for neutral stability still shows $\widetilde {St_r}\sim \zeta _c^2$ with $\widetilde {St_r}\ll 1$ , as depicted in figure 7(b).

Meanwhile, the Floquet exponent $\gamma _{\pi /2}$ indicates the rate at which perturbed trajectories converge to or diverge from the tumbling equilibrium mode. A Floquet exponent with a value close to zero indicates a slow convergence (negative exponent) or divergence (positive exponent) and vice versa. According to (3.14), the Floquet exponents associated with the tumbling modes at $\zeta =0$ fall within the range $[-0.08,0.18] \widetilde {St_r}$ , which are one to two orders of magnitude smaller than $\widetilde {St_r}$ (Einarsson et al. Reference Einarsson, Angilella and Mehlig2014). This indicates that, for small $\widetilde {St_r}$ , the evolution of perturbations near the tumbling mode is slow. Note that an $O(1)$ or smaller value of $\widetilde {St_r}$ may correspond to a large $St_r$ for extreme aspect ratios. These findings are similar to those of Dabade et al. (Reference Dabade, Marath and Subramanian2016), although their study focused exclusively on neutrally buoyant spheroids, whereas the present work investigates spheroids with a sufficiently large particle-to-fluid density ratio.

3.2. Finite $\widetilde {St_r}$

When $\widetilde {St_r}$ is finite, such as $\widetilde {St_r}\sim 1$ or $\widetilde {St_r}\gt 1$ , the weakly inertial form of the orientational dynamics of the spheroid (see (2.14)) becomes invalid. For finite $\widetilde {St_r}$ , the stability analysis should be performed based on (2.11), whose perturbation equation is of the following form:

(3.15) \begin{equation} \frac {\mathrm{d}{}}{\mathrm{d}t}\left [\begin{matrix} \delta \boldsymbol n \\ \delta \boldsymbol \omega \\ \end{matrix}\right ] = \unicode{x1D64C}\boldsymbol{\cdot }\left [\begin{matrix} \delta \boldsymbol n \\ \delta \boldsymbol \omega \\ \end{matrix}\right ]\!. \end{equation}

Here, $\delta \boldsymbol n$ and $\delta \boldsymbol{\omega}$ are the perturbations of orientation and angular velocity, respectively. Here $\unicode{x1D64C}$ is the stability matrix related to $\boldsymbol n^e$ and $\boldsymbol \omega ^e$ . According to Cui et al. (Reference Cui, Liu, Qiu and Zhao2024), the boundary of the equilibrium aligning and the tumbling modes is vague due to the particle inertia, i.e. near $A^2=M^2+\varLambda ^2=1$ spheroids may align or continuously tumble depending on initial conditions. Consequently, it is difficult to determine the boundaries for only the aligning mode or tumbling mode directly. Nevertheless, the aligning mode can be thought of a special tumbling mode with an infinitely long period. In this way, we can simplify aligning mode and tumbling mode into one case. Therefore, we just discuss the stability of logrolling, and tumbling–aligning modes.

3.2.1. Logrolling mode

For the logrolling mode, $\boldsymbol n^e=[0,0,\pm 1]^T$ . Due to the constraint for the unit vector $\boldsymbol n$ , the perturbation along $z$ direction $\delta n_z$ near the logrolling mode $\boldsymbol n^e$ is zero, i.e. $\delta n_z =-n_x^e/n_z^e\delta n_x-n_y^e/n_z^e\delta n_y=0$ . Therefore, $\delta \boldsymbol n$ only has the components of $\delta n_x$ and $\delta n_y$ , and the corresponding perturbation equations with constant coefficients are divided into two independent equation sets, i.e.

(3.16) \begin{align}& \frac {\mathrm{d}{ }}{\mathrm{d}t}\left [\begin{matrix} \delta \omega _x \\ \delta \omega _y \\ \delta n_x \\ \delta n_y \end{matrix}\right ] =\left [\begin{matrix} -\dfrac {1}{St_r C_\perp } & \dfrac {\varLambda }{2} & -\dfrac {\varLambda }{2St_r C_\perp } & -\dfrac {\varLambda }{4} \\[10pt] \dfrac {\varLambda }{2} & -\dfrac {1}{St_r C_\perp } & \dfrac {\varLambda }{4}+\dfrac {M}{St_r C_\perp } & \dfrac {\varLambda }{2St_r C_\perp } \\[10pt] 0 & 1 & 0 & \dfrac {1}{2} \\[10pt] -1 & 0 & -\dfrac {1}{2} & 0 \end{matrix}\right ] \left [\begin{matrix} \delta \omega _x \\ \delta \omega _y \\ \delta n_x \\ \delta n_y \end{matrix}\right ]\!, \end{align}

(3.17) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \frac {\mathrm{d}{\delta \omega _z}}{\mathrm{d}t} =-\frac {1}{St_rC_\parallel }\delta \omega _z . \end{align}

Here, (3.17) has a negative eigenvalue $-1/(St_rC_\parallel )$ , and, therefore, only (3.16) are viewed as primary perturbation equations to reflect the stability of the logrolling mode. The real part of the maximum eigenvalue of the stability matrix in (3.16), denoted as $\xi _m$ , serves as the criterion for the stability. The distributions of $\xi _m$ with different Stokes numbers, aspect ratios and $\zeta$ are shown in figure 8. Similar to the weakly inertial results, the black lines for the neutral stable logrolling mode rise with increasing particle inertia, reflecting that greater fluid inertia is needed to overcome the effect of particle inertia.

Figure 8.

Stability of the logrolling mode of the spheroidal particle computed by (3.16) with (a) $\widetilde {St_r}=1$ and (b) $\widetilde {St_r}=5$ . Here, black solid lines represent zero value of eigenvalues representing boundaries between stable and unstable regions, and magenta dashed–dotted lines represent $\varLambda ^2+M^2=1$ . Here ‘US’ and ‘S’ represent unstable and stable states, respectively.

3.2.2. Aligning and tumbling modes

Regarding the aligning and the tumbling modes in the flow-gradient plane, the primary perturbation equation is

(3.18) \begin{equation} \frac {\mathrm{d}{ }}{\mathrm{d}t}\left [\begin{matrix} \delta \omega _x \\ \delta \omega _y \\ \delta n_z \end{matrix}\right ] =\unicode{x1D64C}(\boldsymbol \omega ^e(t),\boldsymbol n^e(t)) \left [\begin{matrix} \delta \omega _x \\ \delta \omega _y \\ \delta n_z \end{matrix}\right ]\!, \end{equation}

where the components of the stability matrix $\unicode{x1D64C}$ are listed as

(3.19) \begin{align} Q_{11} & =\left(\frac {1}{St_r C_\perp }-\frac {1}{St_r C_\parallel }\right){n_x^e}^2-\varLambda n_x^e n_y^e \omega _z^e-\frac {1}{St_r C_\perp }, \nonumber \\ Q_{12} & =-\varLambda {n_y^e}^2\omega _z^e+\left(\frac {1}{St_r C_\perp }-\frac {1}{St_r C_\parallel }\right)n_x^e n_y^e, \nonumber \\ Q_{13} & =-\frac {\varLambda }{2St_r C_\perp }n_x^e+\left(\frac {1}{St_r C_\perp }-\frac {1}{St_r C_\parallel }\right)\left(\omega _z^e+\frac {1}{2}\right)n_x^e-\varLambda n_y^e {\omega _z^e}^2, \nonumber \\ Q_{21} & =\varLambda {n_x^e}^2\omega _z^e+\left(\frac {1}{St_r C_\perp }-\frac {1}{St_r C_\parallel }\right)n_x^e n_y^e, \nonumber \\ Q_{22} & =\left(\frac {1}{St_r C_\perp }-\frac {1}{St_r C_\parallel }\right){n_y^e}^2+\varLambda n_x^e n_y^e \omega _z^e-\frac {1}{St_r C_\perp }, \nonumber \\ Q_{23} & =\frac {\varLambda }{2St_r C_\perp }n_y^e+\left(\frac {1}{St_r C_\perp }-\frac {1}{St_r C_\parallel }\right)\left(\omega _z^e+\frac {1}{2}\right)n_y^e+\varLambda n_y^e {\omega _z^e}^2+\frac {M}{St_r C_\perp }n_x^e, \nonumber \\ Q_{31} & =n_y^e,\qquad Q_{32}=-n_x^e,\qquad Q_{33}=0. \end{align}

Other independent equation set of perturbation equations is

(3.20) \begin{equation} \frac {\mathrm{d}{}}{\mathrm{d}t}\left [\begin{matrix} \delta \omega _z \\ \delta n_x \\ \delta n_y \end{matrix}\right ] = \left [\begin{matrix} -\dfrac {1}{St_r C_\perp } & \dfrac {\varLambda }{St_r C_\perp }n_x^e-\dfrac {M}{St_r C_\perp }n_y^e & \dfrac {\varLambda }{St_r C_\perp }n_y^e-\dfrac {M}{St_r C_\perp }n_x^e \\[10pt] -n_y^e & 0 & -\omega _z^e \\[6pt] n_x^e & \omega _z^e & 0 \\ \end{matrix}\right ] \left [\begin{matrix} \delta \omega _z \\ \delta n_x \\ \delta n_y \end{matrix}\right ]\!, \end{equation}

which only characterizes the stability of the aligning and tumbling modes confined to the flow-gradient plane, but does not capture departures of the spheroid orientation from the flow-gradient plane. Therefore, we consider the stability of (3.18). According to our previous study (Cui et al. Reference Cui, Liu, Qiu and Zhao2024), the rotation modes in the flow-gradient plane for the spheroids with finite Stokes number depend on the initial orientations and angular velocities within certain ranges of $\sqrt {\varLambda ^2+M^2}, \widetilde {St_r}$ and $\zeta$ . Nevertheless, those results in Cui et al. (Reference Cui, Liu, Qiu and Zhao2024) were obtained without perturbations along the vortical axis. In this study, differing from the quasi-two-dimensional results, the perturbations along the vortical axis are considered. Nevertheless, we still reckon the aligning and tumbling modes in the flow-gradient plane are two equilibria of spheroids in three-dimensional cases. However, it is unfeasible to compute (2.20) by traversing all initial values. To acquire final rotation modes as fast as possible, similar to Cui et al. (Reference Cui, Liu, Qiu and Zhao2024), we set $(\omega _z=0, \varphi +\varphi _0=0)$ and $(\omega _z=-0.5, \varphi +\varphi _0=\pi /2)$ as initial values to obtain time series of spheroids orientations and angular velocities by (2.20). After that, by using Floquet theory numerically, as shown in Appendix C, we can attain a maximum eigenvalue ${\varrho }_m$ as a criterion for stability. Note that although the aligning mode is not periodic, we can still use Floquet theory by considering the period of the aligning mode set as any non-zero values. The phase diagrams of the stability of the tumbling and the aligning modes are shown in figure 9. When using the initial values $(\omega _z=0, \varphi +\varphi _0=0)$ , the corresponding equilibrium is the aligning mode with $M^2+\varLambda ^2\geqslant 1$ , as shown in figure 9(a,b), suggesting stable aligning modes for oblate spheroids and unstable aligning modes for prolate ones. In figure 9(c,d), the results are obtained using initial conditions $(\omega _z=-0.5, \varphi +\varphi _0=\pi /2)$ within the region defined by $M^2+\varLambda ^2\lt A_c^2$ (tumbling mode), where the critical value are $A_c=1.05$ at $\widetilde {St_r}=1$ and $A_c=3.18$ at $\widetilde {St_r}=5$ . The stability is determined by the relative importance of fluid inertia to particle inertia. For small $\zeta$ , namely a weak fluid inertia effect compared with particle inertia, the tumbling mode for oblate spheroids is unstable but remains stable for prolate spheroids. As $\zeta$ increases, for moderate aspect ratios close to unit, the tumbling mode becomes stable for oblate spheroids while unstable for prolate spheroids. Note that as particle inertia increases, the critical $\zeta _c$ , which separates the stable and unstable regimes for the tumbling mode, also increases.

Figure 9.

Stability of (a,b) the aligning and (c,d) the tumbling modes of spheroidal particle computed by (3.18) with finite Stokes numbers (a,c) $\widetilde {St_r}=1$ and (b,d) $\widetilde {St_r}=5$ . Note that the contour represents Floquet exponent $\gamma _m =\ln (\varrho _m)/T_{\!p}$ . The equilibrium rotations in the x–y plane are obtained by numerically integrating (2.20). (a,b) The aligning mode is obtained by using the initial condition ( $\omega _z=0, \varphi +\varphi _0=0$ ), (c,d) while the tumbling mode is evaluated using ( $\omega _z=-0.5, \varphi +\varphi _0=\pi /2$ ) as initial values. Black solid lines represent $\gamma _m=0$ , and magenta dashed–dotted lines represent $\varLambda ^2+M^2=1$ . The green solid lines in (a) and (c) represent $\varLambda ^2+M^2=1.05^2$ and in (b) and (d) represent $\varLambda ^2+M^2=3.18^2$ . Between the green and magenta lines lies a coexistence regime, where both aligning and tumbling modes are possible, and the final state depends sensitively on the initial orientation and angular velocity. Here ‘US’ and ‘S’ represent unstable and stable states, respectively.

3.3. Phase diagrams of stable rotations of spheroids

To summarize, we integrate all stable modes into phase diagrams with $\widetilde {St_r}=0$ , 0.1, 1 and 5 and the aspect ratios $\lambda$ ranging from $0.02$ to $30$ , as shown in figure 10. Note that the results of $\widetilde {St_r}=0$ and 0.1 are computed by (3.3)–(3.14), while those of $\widetilde {St_r}=1$ and $5$ are obtained from the numerical stability analysis based on (3.16) and (3.18).

Figure 10.

Phase diagrams of spheroid stable rotations with Stokes numbers $\widetilde {St_r}$ equal to (a) 0, (b) 0.1, (c) 1 and (d) 5, respectively. Here ‘A’, ‘T’ and ‘L’ represent the aligning, tumbling and logrolling modes, respectively. Black lines denote contours where the corresponding maximum eigenvalue is zero, indicating the neutral stability. The green and magenta lines separate distinct equilibrium-mode regimes. The colour-filled areas indicate the corresponding stable equilibrium modes. The abbreviations are as follows: A, aligning; L, logrolling; T, tumbling.

Figure 10 indicates the complexity of stable rotation modes of spheroids at different $\widetilde {St_r}$ . This is because of the competition between the fluid-inertial effect and the particle inertia. Here, the combinations of characters in different parameter regions represent the possible stable rotation modes of the spheroids, e.g. TL represents the rotation modes in which the spheroids have tumbling or logrolling stable modes. Note that the multistable rotation modes represent there is more than one stable equilibrium rotation modes for certain parameters. This indicates that, with different initial conditions, the rotation of a spheroid may converge to different stable rotation modes. At $\widetilde {St_r}=0$ in figure 10(a), inertia-free spheroids only have monostable modes. As is well known, when $\zeta =0$ , i.e. in the absence of fluid-inertial torque, spheroids follow Jeffery orbits. For $\zeta \gt 0$ , logrolling mode is stable for prolate spheroids, while aligning and tumbling modes separated by $\varLambda ^2+M^2=1$ are stable for oblate spheroids. These results are consistent with our previous study (Cui et al. Reference Cui, Liu, Qiu and Zhao2024). In figure 10(b), where $\widetilde {St_r}=0.1$ , satisfying $\widetilde {St_r}\ll 1$ , we observe that prolate spheroids have tumbling, TL and logrolling stable modes as $\zeta$ increases, while oblate spheroids with moderate aspect ratios have logrolling, TL, tumbling and aligning stable modes as $\zeta$ increases. For $\lambda \lesssim 0.05$ , there is an aligning–logrolling (AL) bistable mode between the regions of aligning and TL modes. Meanwhile, the region of the TL mode of prolate spheroids is much narrower than that of oblate ones, indicating that the condition of the TL mode for the prolate spheroids is more rigorous.

As particle inertia increases (see figure 10 c,d), the phase diagram for different stable rotation modes becomes more complicated. Regarding the oblate spheroids, the tristable ATL and bistable AT regions emerge and are enlarged as Stokes number increases, especially for the oblate spheroids with small aspect ratios. Meanwhile, the boundaries between logrolling and TL regions are prominently moved towards larger $\zeta$ and intersect with magenta lines $\varLambda ^2+M^2=1$ . For prolate spheroids, the TL region is significantly enlarged compared with that at $\widetilde {St_r}=0.1$ . According to figure 10, it is pronounced that inertial prolate spheroids experience tumbling, TL, and logrolling stable modes as $\zeta$ increases. Nevertheless, the parameter space of the TL region is significantly enlarged as Stokes number $\widetilde {St_r}$ increases.

In summary, when the fluid-inertial torque is considered, it is common to find multistable rotation modes of spheroids. These phenomena are mainly observed for oblate spheroids with large fluid inertia, while prolate spheroids only have a bistable TL region. Moreover, as particle inertia increases, the TL region of prolate spheroids is significantly enlarged and overwhelms that of oblate ones. This phase diagram indicates that when the fluid-inertial torque is taken into consideration, the rotations of oblate spheroids are more sensitive to the initial values, and exhibit variety compared with prolate spheroids.