1. Introduction
Particle-laden flow exists widely in nature and engineering. From a fundamental perspective, many particles possess a variety of shapes, e.g. bacteria, viruses, biconcave disc-like human red blood cells and cylindrical Escherichia coli (Daniel, Mehmet & Patrick Reference Daniel, Mehmet and Patrick2007; Park et al. Reference Park, Kwon, Chung, Lee and Shin2008; Hu et al. Reference Hu, Kang, Lin, Lin, Bao and Zhu2023a). Spheroids are the simplest kind of non-spherical particle. Jeffery (Reference Jeffery1922) first studied the motion of an ellipsoidal particle in a Couette shear flow at a century ago. This classical study has received much attention since 2010 due to the recent progress in fluid dynamics in microfluidics, as well as microrobots, and their applications to biological and medical problems (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Kroo et al. Reference Kroo, Jeremy, Noah, Manu and Eric2022; Ouyang et al. Reference Ouyang, Lin, Lin, Yu and Nhan2023). In practice, most fluids in nature show rheological properties, which influence the motion of particles. For example, the transport of smog and virus particles affects health; the sand-carrying technology of drilling fluid during fracturing affects oil production in the petroleum industry (Bird Reference Bird1976; Manoorkar & Morris Reference Manoorkar and Morris2021; Sojwal & Morris Reference Sojwal and Morris2021; Ferroudji et al. Reference Ferroudji, Rahman, Hadjadj, Ofei, Khaled, Rushd and Gajbhiye2022; Hu et al. Reference Hu, Lin, Lin, Zhu and Yu2022). Some fundamental issues, such as the effects of fluid rheology and particle shape on motion behaviour, have always been a focus of attention (Reddy, Tiwari & Singh Reference Reddy, Tiwari and Singh2019; Li et al. Reference Li, Abbas, Morris, Climent and Magnaudet2020; Saccone, Marchioli & Marchis Reference Saccone, Marchioli and Marchis2022).
Lateral migration of particles in channel flow has gained much attention in high-throughput cytometry (Liu et al. Reference Liu, Hu, Jiang and Sun2015; Haddadi & Di Carlo Reference Haddadi and Di Carlo2017; Majji, Banerjee & Morris Reference Majji, Banerjee and Morris2018; Nakayama et al. Reference Nakayama, Yamashita, Yabu, Itano and Sugihara-Seki2019). The classical Segre–Silberberg effect shows that the spherical particles in a circular channel flow of Newtonian fluid migrates radially towards the wall to seek an equilibrium position where the total radial force is zero (Segré & Silberberg Reference Segré and Silberberg1961; Yu & Shao Reference Yu and Shao2010; Hamed, Lee & Di Carlo Reference Hamed, Lee and Di Carlo2014; Yang, Huang & Lu Reference Yang, Huang and Lu2017). The migration of spherical particles has been extensively studied in Newtonian fluid (Shao, Yu & Sun Reference Shao, Yu and Sun2008; Matas, Morris & Élisabeth Reference Matas, Morris and Élisabeth2009; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014). Hood, Lee & Roper (Reference Hood, Lee and Roper2015) and Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) proposed the lift force equation under low-Reynolds-number conditions (20 < Re < 80), which was insufficient in practical applications. Under large fluid inertia, Matas et al. (Reference Matas, Morris and Élisabeth2009) pointed out another equilibrium position closer to the centreline of the circular channel. Miura, Itano & Sugihara-Seki (Reference Miura, Itano and Sugihara-Seki2014) and Shichi et al. (Reference Shichi, Yamashita, Seki, Itano and Sugihara-Seki2017) observed eight equilibrium positions located at four corners and wall bisectors under high Reynolds number (Re). Choi, Seo & Lee (Reference Choi, Seo and Lee2011) suggested that further works were needed to study the changes in the travel distance for particles to reach equilibrium positions in channel flow.
From a practical point of view, viscoelastic properties have been demonstrated to be the dominant behaviour in non-Newtonian fluids, and previous studies have mainly focused on the motion of spherical particles in non-Newtonian fluids (Pasquino et al. Reference Pasquino, D'Avino, Maffettone, Greco and Grizzuti2014; D'Avino, Del Greco & Maffettone Reference D'Avino, Del Greco and Maffettone2017). The migration of spherical particles in a circular channel flow of a viscoelastic fluid was first studied by Leshansky et al. (Reference Leshansky, Bransky, Korin and Dinnar2007). Del Giudice et al. (Reference Del Giudice, D'Avino, Greco, Maffettone and Shen2018) found that the gradients of normal stress differences in viscoelastic fluid yielded an elastic lift, which pushed particles towards the low-shear-rate regions when the fluid inertia was negligible. Villone et al. (Reference Villone, D'Avino, Hulsen, Greco and Maffettone2013) reported that the particle in a Giesekus fluid could be trapped by the secondary flow, the cross-streamline migration velocity overcame the velocity of secondary flow as the Weissenberg number (Wi) was reduced, and the particles migrated towards the channel centreline. Hu et al. (Reference Hu, Lin, Lin, Zhu and Yu2022) found that the polydisperse spherical particles in the square channel flow of viscoelastic fluid would migrate towards the channel centreline and form particle chains. Generally, the previous studies mainly considered the condition under low Reynolds number (D'Avino et al. Reference D'Avino, Del Greco and Maffettone2017; Liu et al. Reference Liu, Guo, Tian, Yang, Yan, Ding, Wei, Hu, Nie and Sun2017; Raoufi et al. Reference Raoufi, Mashhadian, Niazmand, Asadnia, Razmjou and Warkiani2019). The fluid inertial effect was found to be beneficial for particles separation in a viscoelastic fluid, called elasto-inertial particle focusing (Lim et al. Reference Lim, Ober, Edd, Desai, Neal, Bong, Doyle, Mckinley and Toner2014; Seo, Kang & Lee Reference Seo, Kang and Lee2014). When considering the effect of fluid inertia, Kim & Kim (Reference Kim and Kim2016) found that the elastic lift in viscoelastic fluid pushed the spherical particles towards the corners and centreline of the square channel, whereas the fluid inertia pulled the particles out of the corners. Yu et al. (Reference Yu, Wang, Lin and Hu2019) discovered that the spherical particle in a square channel would focus at corner, wall centreline, diagonal line and channel centreline equilibrium positions. Li, Mckinley & Ardekani (Reference Li, McKinley and Ardekani2015) found that the motion of spherical particles in an Oldroyd-B viscoelastic fluid depended on the interplay between the elastic and inertial effects, and the channel centreline (CC) equilibrium position occurred in flows with strong elasticity. For inelastic power-law fluids, Hu et al. (Reference Hu, Lin, Chen and Ku2020) found that the migration of spherical particles in a square channel was similar to that of Newtonian fluids, the number of equilibrium positions did not change when Re < 300.
The aforementioned proposed studies were, however, limited to spherical particles. Because there is no asymptotic analytical solution for non-spherical particles, the motion of ellipsoidal particles in viscous fluids is complicated. The rotational mode and its orientation distribution of spheroids are strongly coupled with the rotational behaviours, and the inhomogeneity of the particle orientation distribution also affects the flow field. One of the classic investigations was the motion of ellipsoidal particles in Couette shear flow of a Newtonian fluid (Jeffery Reference Jeffery1922; Qi & Luo Reference Qi and Luo2003). Rosén, Lundell & Aidun (Reference Rosén, Lundell and Aidun2014) found the log-rolling (LR), tumbling (TU), inclined-rolling (ILR), kayaking and inclined-kayaking rotational modes for the motion of a prolate spheroid in Couette flow. Huang et al. (Reference Huang, Yang, Krafczyk and Lu2012) found that an oblate spheroid would maintain the LR and ILR rotational modes, and that the rotational modes were insensitive to the initial orientation. In a non-Newtonian fluid of Couette shear flow, it was observed that the ellipsoidal cells orientated between the vorticity and the flow direction under a moderately elastic effect. In highly elastic fluids, the fibre-like particles aligned their symmetry axes along the flow direction (Iso, Cohen & Koch Reference Iso, Cohen and Koch1996a,Reference Iso, Koch and Cohenb; Lyon et al. Reference Lyon, Mead, Elliott and Leal2001). Recently, Li, Xu & Zhao (Reference Li, Xu and Zhao2023) discovered five rotational modes of a prolate spheroid in a Couette shear flow of viscoelastic fluids, where the prolate spheroid rotated with a newly asymmetric-kayaking mode as the fluid elasticity increased. Moreover, Dabade, Marath & Subramanian (Reference Dabade, Marath and Subramanian2015) proposed a torque model based on the generalized reciprocal theorem for both prolate and oblate spheroids in viscoelastic fluids, in which the first normal stress difference was found to be the dominant effect. D'Avino et al. (Reference D'Avino, Hulsen, Greco and Mafettone2014) observed four different motion modes of a prolate spheroid in an unconfined Couette shear flow of viscoelastic fluids as the Wi increased, i.e. log-rolling, metabistability, bistability and flow aligning.
Due to the slow process for the lateral migration of the particles in a channel flow, there are few related studies on the migration and control of the ellipsoidal particle in channel flow. In a Newtonian fluid, Huang & Lu (Reference Huang and Lu2017) reported that there were five types of rotational behaviours for prolate ellipsoidal particle migration in circular channel flow, and the rotational behaviours were mainly dependent on the fluid inertia, particle size and initial orientation. Hur, Tse & Di Carlo (Reference Hur, Tse and Di Carlo2010) experimentally found that the prolate spheroid in square channel flow predominantly maintained the tumbling rotational mode when Re was above a critical value (Re ≈ 50). Masaeli et al. (Reference Masaeli, Sollier, Amini, Mao, Camacho, Doshi, Mitragotri, Alexeev and Di Carlo2012) concluded that the cylindrical and ellipsoid particles in a microchannel flow preferred to tumble, while disk-like particles preferred to roll. Lashgari et al. (Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017) concluded that the final orientation of oblate particles exhibited chaotic behaviour when Re was beyond a critical value. Li, Xia & Wang (Reference Li, Xia and Wang2022) found that the oblate ellipsoidal particle in a square channel maintained a stable corner equilibrium position only when the Reynolds number was small. In viscoelastic fluid, D'Avino et al. (Reference D'Avino, Hulsen, Greco and Maffettone2019) found that a prolate ellipsoidal particle in a wide-slit microchannel of viscoelastic fluid would migrate to the closest wall or the channel centreline, depending on the initial position and fluid rheology, but the effects of fluid inertia and particle shape were not studied in their works. For the migration of spheroids in the channel flow of inelastic power-law fluid, Hu et al. (Reference Hu, Kang, Lin, Lin, Bao and Zhu2023a,Reference Hu, Lin, Lin and Zhub) found that the prolate spheroid exhibited TU and LR modes, and the LR mode was conditionally stable. For the oblate spheroid, only the LR mode existed when Re was large, while TU and ILR modes were observed when Re was decreased.
To summarize, although there have been a few studies conducted on the rotational dynamics of non-spherical particles in viscoelastic Couette shear flow, these studies are mostly concerned with fluid elasticity, where fluid inertia is absent. In addition, only a few works have been conducted concerning the inertial migration of non-spherical particles in the channel flow of power-law fluids, whose constitutive equation is relatively simple. The interplay between fluid elasticity and inertial effects on spheroid migration in channel flow has not yet been explored, and the corresponding rotational behaviour is still unknown. Therefore, the direct forcing/fictitious domain method is used in this study to determine how particle shape, fluid elastic effect and fluid inertia can together affect particle motion and rotational behaviour in channel flow. Moreover, the effect of the Wi and Re of the fluid, the effect of the particle aspect ratio on the elasto-inertial focusing and rotating characteristics of the spheroid, the number of stable equilibrium positions and the travel distance are systematically explored. These results are helpful for understanding the motion of non-spherical particles in channel flow and pave the way for designing efficient microfluidic devices.
The outline of this paper is as follows. The numerical methods adopted are described in § 2, with the validation results are reported in Appendix A. In § 3, the main results and discussions are presented. We study the mechanisms of elasto-inertial focusing and rotating characteristics of a spheroid in a square channel flow of Oldroyd-B viscoelastic fluids. A summary of the findings is given in § 4.
2. Numerical model
2.1. Simulation set-up
The elasto-inertial migration of neutrally buoyant particles in a straight square channel flow of Oldroyd-B viscoelastic fluids is shown in figure 1(a). The constant pressure gradient along the z direction is introduced to sustain the channel flow, the periodic boundary condition is adopted along the flow direction and the no-slip boundary condition is applied on the walls. The channel length, height and width are denoted as L, H and H, respectively. The computational domain spans [−L/2, L/2] × [−H/2, H/2] × [−H/2, H/2]. The Reynolds number is defined as Re = U 0ρfH/
$\mu$0, where U 0 is the maximum velocity of the fluid at the channel centreline, 
$\mu$0 is the total zero-shear-rate viscosity of the fluid and ρf is the fluid density. The Weissenberg number is defined as Wi = λtU 0/H (λt being the fluid relaxation time). The elasticity number, El, is the ratio of the inertial and elastic effects, i.e. El = Wi/Re = λt 
$\mu$0/ρfH 2.
Figure 1. Simulation set-up. (a) Schematic of particle migration in three-dimensional square channel flow of Oldroyd-B viscoelastic fluids; the space-fixed coordinate system and body-fixed coordinate system for the (b) prolate spheroid and (c) oblate spheroid. The sphere, prolate and oblate spheroids are represented by blue, grey and pink colours, respectively.
The space-fixed coordinate system and body-fixed coordinate system for the prolate spheroid and oblate spheroid are shown in figure 1(b,c). The spheroid is described by 
$x^{\prime 2}/{a^2} + y^{\prime 2}/{b^2} + z^{\prime 2}/{c^2} = 1$, where (x′, y′, z′) denotes a body-fixed coordinate system; a, b and c are the lengths of the three semi-axes of the particles. The aspect ratio is defined as α = b/a, where α = 1, α < 1 and α > 1 correspond to the cases of spheres (a = b = c), oblate ellipsoidal particles (a = c) and prolate ellipsoidal particles (a = b), respectively. The blockage ratio of the sphere, prolate and oblate ellipsoidal particles are defined as k = 2b/H. The blockage ratio is kept as k = 0.2 for all the ellipsoidal particles. The mesh convergence is compared in our previous work (Liu et al. Reference Liu, Lin, Ku and Yu2021). We further found that the effect of the domain length L on the mesh convergence is insensitive when L = 256 and 384 (the mesh size is h = H/128). The effect of discretization of the particles is not significant when the mesh size is h = H/128 and H/256 in the present work. Thus, the mesh numbers along the z, x and y directions with 256, 128 and 128 are adopted to balance the computational accuracy and efficiency. Taking the initial position of spherical particles at the channel cross-section (xin, yin) as an example, the closest initial position of the spheroid mass centre to the wall is uniformly set as 0.395 to avoid the initial overlap between the particle surface and the wall.
Due to the geometric symmetry of the square channel, the particles are released only in the upper right quadrant of the square channel cross-section (in the red triangle of figure 1a) with various initial positions. We identify four kinds of final equilibrium positions for the elasto-inertial migration of spheroids in the square channel flow of an Oldroyd-B viscoelastic fluid. The final equilibrium positions located on the cross-section midline (CSM), diagonal line (DL), corner (CO) and channel centreline (CC) are represented in figure 1(a) by circles with black, blue, red and green colours, respectively.
2.2. Direct forcing/fictitious domain method
The elasto-inertial focusing of neutrally buoyant particles in a square channel flow of the Oldroyd-B viscoelastic fluid is studied by the direct forcing/fictitious domain (DF/FD) method. For the motion of particles in both Newtonian and non-Newtonian fluids, the present method was described in detail and successfully verified in our previous publications (Yu & Shao Reference Yu and Shao2007; Wang, Yu & Lin Reference Wang, Yu and Lin2018; Yu et al. Reference Yu, Wang, Lin and Hu2019; Liu et al. Reference Liu, Lin, Ku and Yu2020; Yu et al. Reference Yu, Xia, Guo and Lin2021; Hu et al. Reference Hu, Lin, Lin, Zhu and Yu2022), and the detail of the method can be seen in Yu & Shao (Reference Yu and Shao2007) and Yu et al. (Reference Yu, Wang, Lin and Hu2019). Here, the dimensionless FD formulations for the incompressible Oldroyd-B fluid comprise the following three parts.
(1) Combined momentum equations:
(2.1)
\begin{gather}\frac{{\partial \boldsymbol{u}}}{{\partial t}} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = \frac{{{\mu _r}{\nabla ^2}\boldsymbol{u}}}{{Re}} - \boldsymbol{\nabla }p + \frac{{(1 - {\mu _r})\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B}}}{{Re\,Wi}} + \boldsymbol{\lambda }\quad \textrm{in}\;{\varOmega },\end{gather}(2.2)\begin{gather}\boldsymbol{u} = \boldsymbol{U} + {\omega _p} \times \boldsymbol{r}\quad \textrm{in}\;S(t),\end{gather}(2.3)\begin{gather}({\rho _r} - 1)V_p^\ast \frac{{\textrm{d}\boldsymbol{U}}}{{\textrm{d}t}} ={-} \int_P {\boldsymbol{\lambda }\,\textrm{d}\kern0.7pt\boldsymbol{x}} ,\end{gather}(2.4)\begin{gather}({\rho _r} - 1){\boldsymbol{J}^\ast }\frac{{\textrm{d}({\boldsymbol{\omega }_p})}}{{\textrm{d}t}} ={-} \int_P {\boldsymbol{r} \times \boldsymbol{\lambda }\,\textrm{d}\kern0.7pt\boldsymbol{x}} .\end{gather}(2) Continuity equation:
(2.5)\begin{equation}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0\quad \textrm{in}\;\varOmega ,\end{equation}where u and p are the fluid velocity and pressure, respectively; λ is the Lagrange multiplier that is defined in the solid domain S(t);
$\mu$r is the ratio of the solvent viscosity ($\mu$s) to the total zero-shear-rate viscosity of the fluid ($\mu$0); B is the configuration tensor; $\varOmega$ is the entire domain including the interior and exterior of the solid particle; U and ωp are the particle translational velocity and angular velocity; r is the position vector with respect to the mass centre of the particle; ρr is the particle–fluid density ratio, ρr = ρs/ρ (ρs is the particle density), here ρr = 1; $V_p^\ast= {V_p}/{H^3}$ is the dimensionless particle volume, Vp is the particle volume; J* = J/ρsH 5 is the dimensionless moment of inertia.(3) Oldroyd-B constitutive equation:
(2.6)\begin{equation}\frac{{\partial \boldsymbol{B}}}{{\partial t}} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B} - \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} - {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}}\boldsymbol{\cdot }\boldsymbol{B} + \frac{{\boldsymbol{B} - I}}{{Wi}} = 0\quad \textrm{in}\;\varOmega .\end{equation}Constitutive equation sub-problem for B is updated by
(2.7)\begin{align}\frac{{{\boldsymbol{B}^{n + 1}} - {\boldsymbol{B}^n}}}{{\mathrm{\Delta }t}} + {\boldsymbol{u}^{n + 1}}\boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol{B}^n} - {\boldsymbol{B}^n}\boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol{u}^{n + 1}} - {(\boldsymbol{\nabla }{\boldsymbol{u}^{n + 1}})^\textrm{T}}\boldsymbol{\cdot }{\boldsymbol{B}^n} + \frac{{{\boldsymbol{B}^{n + 1}} - \boldsymbol{I}}}{{Wi}} = 0,\end{align}where un +1 has been determined from the first two sub-problems. This equation is solved by the first-order time scheme, the central difference scheme for the velocity gradient and the third-order upwinding MUSCL scheme for the convective term.
To realize the transformation of torque and angular velocity between the body-fixed coordinate system and space-fixed frame, four quaternion variables are introduced:
\begin{equation}\left. {\begin{array}{@{}c@{}} {{q_1} = \textrm{sin}(0.5\theta )\,\textrm{cos}(0.5(\phi - \varphi )),}\\ {{q_2} = \textrm{sin}(0.5\theta )\,\textrm{sin}(0.5(\phi - \varphi )),}\\ {{q_3} = \textrm{cos}(0.5\theta )\,\textrm{sin}(0.5(\phi + \varphi )),}\\ {{q_4} = \textrm{cos}(0.5\theta )\,\textrm{cos}(0.5(\phi + \varphi )),} \end{array}} \right\}\end{equation}
where 
$\boldsymbol{n} = (\theta ,\phi ,\varphi )$ is the Euler angle, which is used to show the initial orientation angle in the body-fixed coordinate system. The orientation angle in a space-fixed frame can be obtained as 
$\boldsymbol{P} = ({P_x},{P_y},{P_z}) = (2({q_1}{q_3} + {q_2}{q_4}),\;2({q_2}{q_3} - {q_1}{q_4}),\;2(q_3^2 + q_4^2 - 0.5))$.
The transformation matrix M from the space-fixed frame to the body-fixed frame is
\begin{equation}M = 2\left[ {\begin{array}{@{}ccc@{}} {q_1^2 + q_4^2 - 0.5}&{{q_1}{q_2} + {q_3}{q_4}}&{{q_1}{q_3} - {q_2}{q_4}}\\ {{q_1}{q_2} - {q_3}{q_3}}&{q_2^2 + q_4^2 - 0.5}&{{q_2}{q_3} + {q_1}{q_4}}\\ {{q_1}{q_3} + {q_2}{q_4}}&{{q_2}{q_3} - {q_1}{q_4}}&{q_3^2 + q_4^2 - 0.5} \end{array}} \right].\end{equation}In the simulation, (2.8) and (2.9) are used to map the hydrodynamic moment vectors and the angular velocity vectors between the inertial and body-fixed frames. Then the fourth-order Runge–Kutta method is adopted to obtain each quaternion (q) and angular velocity.
3. Results and discussion
3.1. Mode classification
Within the present simulated parameters (1 ≤ Re ≤ 100, 0 ≤ Wi ≤ 2, 0.4 ≤ α ≤ 3), there are six (five) kinds of rotational behaviours for the elasto-inertial focusing of prolate (oblate) spheroids in a square channel flow of Oldroyd-B viscoelastic fluids. All the rotational modes are shown in figure 2.
(1) Tumbling (TU) mode: the prolate (oblate) spheroid rotates around the short (long) axis, and the spheroid rotates in the YOZ plane. The corresponding rotational behaviours are shown in figure 2(a,g).
(2) Tumbling mode in a diagonal plane (TUD): the spheroid migrates to the diagonal plane and keeps tumbling, which is a special case of the TU mode. The corresponding rotational behaviours are shown in figure 2(b,h).
(3) Log-rolling (LR) mode: the prolate (oblate) spheroid eventually rotates around its long (short) axis, which is almost perpendicular to the flow direction. The corresponding rotational behaviours are shown in figure 2(c,i).
(4) Inclined log-rolling (ILR) mode: the prolate (oblate) spheroid eventually rotates around its long (short) axis, and the orientation is parallel to the diagonal line of the square channel. The corresponding rotational behaviours are shown in figure 2(d,j).
(5) Vibrated inclined log-rolling (VILR) mode: here, the prolate spheroid also rotates like the ILR mode, but the spheroid vibrates along the long axis (see figure 2e), which is newly found in the present work as a special case of the ILR mode. This rotational mode is unique to the migration of the prolate spheroid.
(6) Kayaking mode: the prolate (oblate) spheroid rotates around the short (long) axis with precession and nutation (see figure 2f,k) between the tube axis and the wall. It looks like a kayak paddle moving around its equilibrium position with a small amplitude.
Figure 2. Rotational mode classification for the prolate (first line) and oblate (second line) spheroids with (a,g) tumbling (TU) mode; (b,h) tumbling mode in a diagonal plane (TUD); (c,i) log-rolling (LR) mode; (d,j) inclined log-rolling (ILR) mode; (e) vibrated inclined log-rolling (VILR) mode; ( f,k) kayaking mode.
To describe the rotational motion modes in more detail, the orientation in the space-fixed frame can be obtained as P = (Px, Py, Pz), and the long-time evolution of the orientation in the space-fixed frame is shown in the following sentences.
3.2. Spheroid migration under a low Reynolds number
The previous studies are mostly confined to inertia-free spherical particles immersed in viscoelastic flows (D'Avino et al. Reference D'Avino, Del Greco and Maffettone2017; Liu et al. Reference Liu, Guo, Tian, Yang, Yan, Ding, Wei, Hu, Nie and Sun2017), so we first study the migration and rotational behaviour of ellipsoidal particles at a low Reynolds number (Re = 1). In the general initial positions, the initial orientation of spheroid has negligible effect. Thus, the general initial Euler angle n = (
${\rm \pi}$/2, 
${\rm \pi}$/4, 0) is adopted in §§ 3.2, 3.3 and 3.4 to analyse the effect of general initial particle position. Only when the particles are initially located at special positions (corner and wall bisector) do some special initial orientations of the spheroids have an impact on the final equilibrium position and rotational mode. The effect of initial Euler angle on the rotational motion of spheroid at the corner and wall bisector is compared in § 3.5 in detail.
To intuitively obtain how fast the particle migrates in the channel, the trajectories are marked by circles for every 10 000 dimensionless time units. Thus, there is more than one circle on the trajectory when the migration velocity is extremely slow. Figure 3 shows the trajectory of particles on the cross-section at different Wi and initial positions (represented by the purple dashed circle). The background contour in figure 3 is the first normal stress difference (N 1) of the flow field at the particle mass centre in the XOY section when the particles migrate to the CC equilibrium position. Comparing the results at different Wi, it can be found that the regions of low first normal stress difference are generated in the channel centre and the corners, and the maximum first normal stress difference increases with increasing Wi. The orientation for particles initially placed at a corner (xin, yin) = (0.395, 0.395) is shown in figure 4.
Figure 3. Effect of the Weissenberg number on the trajectories of the prolate and oblate spheroids in the channel cross-section with Re = 1. (a) Wi = 0.01, α = 2 (top figure), α = 0.5 (bottom figure); (b) Wi = 1, α = 2; (c) Wi = 2, α = 2; (d) Wi = 0.5, α = 0.5; (e) Wi = 1, α = 0.5; ( f) Wi = 2, α = 0.5.
Figure 4. Long-time evolution of the orientation for the (a) prolate spheroid in the LR, VILR and ILR modes and (b) oblate spheroid in the LR, kayaking and TUD modes when the Weissenberg number is reduced from top to bottom, the particle is initially placed at a corner with (xin, yin) = (0.395, 0.395).
When the Wi is small (Wi = 0.01), as shown in figure 3(a), the first normal stress difference of the flow field is weak, and the elastic migration velocity of the spheroid is very slow. Therefore, the migration time to the CC equilibrium position with a stable rotational behaviour at the bottom of figure 4 is the largest. As shown at the bottom of figure 4(a), the prolate spheroid initially located at a corner maintains the kayaking mode for a long time, and it finally transforms into the ILR mode and migrates diagonally to the CC equilibrium position. At the bottom of figure 4(b), the oblate spheroid initially located at a corner maintains the TUD mode for a long time; consequently, the oblate spheroid finally maintains the LR mode when it migrates diagonally to the CC equilibrium position. Increasing Wi, as compared with the trend in figure 3(b,c,e,f), the value of the first normal stress difference near the wall increases, an obvious low first normal stress difference region is formed at a corner, and the particle shape and Weissenberg number affect the migration behaviours when the particles are initially placed near the corner. The particles whose initial positions are far from a corner migrate to the CC equilibrium position. The rotational motion of the spheroid results in the inability of the particles to stabilize in the regions of low first normal stress difference at a corner. For example, when Wi = 1, the prolate spheroid near a corner migrates diagonally towards the CO equilibrium position and maintains the VILR mode, as shown in figure 3(b) and the middle of figure 4(a). When Wi =0.5 (figure 3d and the middle of figure 4b), the oblate spheroid can stabilize at the corners and maintain the kayaking mode. When Wi = 1, the oblate spheroid near the corner migrates along the diagonal line with the unstable kayaking mode. Finally, it migrates to the wall near the corners with the LR mode, as shown in figure 3(e), figure 8(i) and the top of figure 4(b). Continuing to increase the Weissenberg number as Wi = 2, the CO equilibrium position disappears, and both the prolate spheroid and oblate spheroid in figure 3(c,f) migrate to the CC equilibrium position and maintain the LR mode. Overall, the spheroids with various initial positions migrate in the flow of the minimum first normal stress difference, i.e. towards the centreline or a corner of the square channel under a low inertial effect.
The travel distance of particles migrating to the final equilibrium position is an important parameter for particle separation in designing the length of microfluidic devices. We choose the initial position (xin, yin) = (0.2, 0.395) as a typical example to further analyse the changes in travel distance for the prolate and oblate spheroids to migrate to the channel centreline. Figure 5 shows that the larger the Wi, the faster the particles migrate to the CC equilibrium position and the shorter is the travel distance. The prolate spheroid maintains the VILR mode along the diagonal orientation with a small rotational amplitude, while the oblate spheroid maintains the TUD and kayaking modes with a large rotational amplitude. The oblate spheroid in figure 5(b) will migrate towards the corner along the x-direction, causing a non-monotonic change in the initial trajectory in the y-direction. Due to the low normal stress difference at the channel centre, the spheroid ultimately stabilizes at the CC equilibrium position. As a result, the travel distance of the oblate spheroid is slightly larger than that of the prolate spheroid. A greater first normal stress difference means the faster the particles migrate to the low shear rate direction and the shorter the channel needed for the particles to migrate to the CC equilibrium position. A larger Wi means the less obvious the effect of the particle shape. Therefore, the travel distance of particles to the CC equilibrium position is mainly affected by the elasticity of the fluid. When the particles gradually approach the diagonal line, the trajectory in the x and y directions remains unchanged, so the trajectory in the YOZ plane is presented in the following sections.
Figure 5. Effect of the Weissenberg number on the travel distance for the (a) prolate (α = 2) and (b) oblate (α = 0.5) spheroids when Re = 1 and (xin, yin) = (0.2, 0.395).
3.3. Effect of the particle aspect ratio
For ellipsoidal particles, the aspect ratio is an important parameter of the particle shape. We further analyse the migration behaviour and rotational mode of the prolate and oblate spheroids with different aspect ratios. The background contour in figure 6(a–c) is the first normal stress difference (N 1) of the flow field at the particle mass centre in the XOY section when the particles migrate to the channel centre. The changes in orientation when the particles migrate to the CO equilibrium position is compared in figure 6(e).
Figure 6. Effect of the particle aspect ratio on spheroid migration with Wi = 1 and Re = 1: (a) oblate spheroid with α = 0.6; (b) sphere with α = 1 and (c) prolate spheroid with α = 3. (d) Phase diagram of the final equilibrium position; here, the initial position in the x-axis is changed from 0.05 to 0.395, and the initial position in the y-axis is fixed as 0.395. (e) Long-time evolution of the orientation angle for the spheroid with (xin, yin) = (0.395, 0.395), and α = 0.6, 1.5, 3. The corresponding rotational modes are kayaking, TUD and ILR modes, respectively.
As shown in figure 6(b), the motion trajectory of the spherical particle is a smooth curve. The ellipsoidal particle maintains the rotational motion with a fluctuating trajectory on the cross-section. Spheroids with different aspect ratios all have CC and CO equilibrium positions, which are consistent with the results in figure 3(b,e). By comparing figure 3 with figure 6(a–c), it is found that a greater aspect ratio of the spheroids results in the smaller the region in the square channel cross-section where the spheroid can migrate to the CO equilibrium position. By changing xin from 0.05 to 0.395 and fixing yin = 0.395, we compare whether spheroids with different aspect ratios can migrate towards the centre or the corners. According to the final equilibrium position, a phase diagram regarding the aspect ratio and the initial particle position is drawn in figure 6(d). When 0.4 ≤ α ≤ 0.6 (1 ≤ α ≤ 3), a flatter (slenderer) spheroid means the larger (smaller) area required for the spheroid to stabilize at the corners. When 0.6 ≤ α ≤ 1, as shown at the top of figure 6(e), the oblate spheroid migrates diagonally in a stable kayaking mode. When 1 < α ≤ 1.75, as shown in the middle of figure 6(e), the prolate spheroid migrates diagonally and maintains the TUD mode at the CO equilibrium position. When 2 ≤ α ≤ 3, the prolate spheroid maintains the VILR and ILR modes, and the particle mass centre remains stable at the corners only when the particle is initially placed near the diagonal line. The rotational tumbling motion of prolate spheroid makes it easier to deviate from the corners and eventually migrate to the CC equilibrium position, resulting in a smaller region where the spheroids can migrate to the corners.
When particles migrate to the channel centreline (Wi = 1, Re = 1, (xin, yin) = (0.2, 0.395)), the effect of the particle aspect ratio on the travel distance (Z*) to the CC equilibrium position is shown in figure 7. Since there is no rotational motion for spherical particles, the travel distance of spherical particles to the CC equilibrium position is the shortest. When the aspect ratio is less (greater) than 1, the travel distance decreases (increases) as the aspect ratio increases; that is, a flatter and slenderer spheroid means the greater the travel distance to the CC equilibrium position.
Figure 7. Effect of the aspect ratio on the changes in the travel distance towards the CC equilibrium position.
We further analyse the interaction mechanism between the particles and the flow field. According to the Faxen law (Faxén Reference Faxén1922), a finite-sized particle in a parabolic Stokes flow exhibits only streamwise translational slip velocity with respect to the local unperturbed flow, u∞ (Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014). In addition, some of the fluid must flow past the particle through the space between the particle and channel walls. As a result, the incoming fluid tends to flow inwards from the channel walls to the shadow of the particle on the cross-section ahead of the particle, the particle-induced secondary flow is formed in the channel cross-section and the flow field away from the particles maintains zero velocity along the y-axis (v = 0). In figure 8, the velocity contour and the flow streamlines are plotted on the typical cross-section in front of the particle centre by one sphere diameter or long axis length of the spheroids. As shown in figure 8(c), spherical particles induce symmetrical streamlines that converge towards the four corners. When the prolate spheroid and oblate spheroid migrate to the CC equilibrium position, the prolate spheroid in figure 8(a,b) is oriented along the diagonal line; the oblate spheroid in figure 8( f,g) maintains the orientation perpendicular to the y-axis. Due to the shape of the spheroid, butterfly-shaped streamlines are induced in figure 8(a,b,f,g), which are different from those of the spherical particle. When the particle migrates to the CO equilibrium position, as shown in figure 8(d,e,h–j), a region of lower velocity is formed near the area between the particle and corner due to the interaction between the particles and the fluids near the walls, and secondary flow is induced near the corners. A region of high velocity is formed far from the corners, where the induced streamlines push the particles to the CO equilibrium position. This phenomenon is more obvious for the oblate spheroid with α = 0.6 and prolate spheroid with α = 1.5, as shown in figure 8(e,j). Therefore, under the interaction of positive and negative fluid velocities around the particles, the oblate (prolate) spheroid forms the stable kayaking (TUD) mode at the CO equilibrium position.
Figure 8. Velocity contours and the streamlines when the particles at the final equilibrium position with Re = 1 move into the paper. The particle maintains the CC equilibrium position with (Wi, α) = (a) (0.01, 2), (b) (1, 2), (c) (1, 1) and the CO equilibrium position with (Wi, α) = (d) (1, 2), (e) (1, 1.5). The oblate spheroid maintains the CC equilibrium position with (Wi, α) = ( f) (0.01, 0.5), (g) (0.5, 0.5) and the CO equilibrium position with (Wi, α) = (h) (0.5, 0.5), (i) (1, 0.5), ( j) (1, 0.6). The typical cross-section ahead of the particle centre by a distance of the major axis of the spheroid is chosen.
When particles migrate to the CC equilibrium position, the velocity contours in the XOY plane are shown in figure 8(a–c,f,g). The particle-induced velocity increases as Wi increases, so the travel distance for the particles migrating to the CC equilibrium position decreases as Wi increases. According to the Faxen law (Faxén Reference Faxén1922), the drag force, Fd, exerted by the fluid on a spherical particle is Fd = 6
${\rm \pi}$
$\mu$sR[u∞ − Us] + 
${\rm \pi}$
$\mu$sR 3∇2u∞, where Us is the streamwise translational slip velocity of the particle and R is the radius of sphere. After the particles migrate to the CC equilibrium position with Fd = 0, the velocity can be evaluated as u∞ − Us ∼ (−R 2∇2u∞/6), which is inversely proportional to the particle size. Thus, particles with small radius will move faster than the large particles. Due to the different rotational behaviour when (Re, Wi) = (1, 1), the rotational diameter of the spherical particles is the largest, then followed for the oblate spheroid and finally for the prolate spheroid. Therefore, the spherical particles migrate fastest to the CC equilibrium position, and rounder ellipsoidal particles result in a shorter travel distance.
3.4. Effect of elastic number
The elastic number of a viscoelastic fluid can explain the effect of fluid elasticity. This article further analyses the elasto-inertial migration behaviours of prolate and oblate spheroids at different elastic numbers. The background contours in figures 9 and 11 are the first normal stress difference (N 1) of the flow field at the particle mass centre in the XOY section when the particles are at the corresponding equilibrium position.
Figure 9. Effect of elastic number on the trajectories of the prolate (first line, α = 2) and oblate (second line, α = 0.5) spheroids in the channel cross-section when Re = 10 with (a,d) Wi = 0.01, El = 0.001; (b,e) Wi = 0.1, El = 0.01; (c,f) Wi = 1, El = 0.1.
When the elastic number is small, El = 0.001 (Re = 10, Wi = 0.01), as shown in figure 9(a,d), the particles migrate away from the corners and eventually reach the DL equilibrium position. The presence of a particle affects the distribution of the first normal stress difference. A negative first normal stress difference appears on the particle side close to the channel centreline, and a positive first normal stress difference appears on the opposite side of the particle. In addition, the first normal stress difference caused by the elastic effect is small, the migration of spheroid towards the final equilibrium position is extremely slow and a particle annulus is formed at approximately 0.28H away from the channel centreline (the shaded area in figure 9a,d), consistent with the experimental observation (Choi et al. Reference Choi, Seo and Lee2011; Abbas et al. Reference Abbas, Magaud, Gao and Geoffroy2014) in Newtonian fluid with Re = 10. As shown at the bottom of figure 10, the prolate spheroid maintains the TUD mode, and the oblate spheroid initially maintains the kayaking mode and eventually transforms into the stable ILR mode.
Figure 10. Long-time evolution of the orientation for the (a) prolate spheroid in the ILR, ILR and TUD modes and (b) oblate spheroid in the LR, kayaking and ILR modes, the elastic number is reduced from top to bottom. The particle is initial placed at a corner with (xin, yin) = (0.395, 0.395) and Re = 10.
When El = 0.01 (Re = 10, Wi = 0.1), the CO equilibrium position in figure 9(b,e) still does not exist, and the particle shape affects the final equilibrium position. As shown in the middle of figure 10(a), the prolate spheroid maintains the kayaking mode for a long time, the rotational period continues to increase as the particles approach the channel centreline and the ILR mode is finally formed. Since the prolate spheroids migrate extremely slowly to the CC equilibrium position, only the motion within 4000 dimensionless time steps is given in the middle of figure 10(a) to present the rotational mode more clearly. It is interesting to note in figure 9(e) and the middle of figure 10(b) that the oblate spheroid affects the flow field around the particles and enlarges the influence region of the first normal stress difference near the channel centre. Finally, the oblate spheroid migrates to the DL equilibrium position and maintains the kayaking mode, different from the prolate spheroid in figure 9(b).
When El = 0.1 (Re = 10, Wi = 1), as shown in figure 9(c,f), the contour of the first normal stress difference is similar as the results in figures 3 and 6(a–c). Thus, the elastic effect dominates, where the particles at a corner migrate to the CO equilibrium position again. As shown at the top of figure 10, the prolate spheroid finally maintains the ILR mode, and the oblate spheroid finally keeps the LR mode near a corner. Both the prolate spheroid and the oblate spheroid away from the corners migrate towards the CC equilibrium position. The final DL equilibrium position of the spheroids in figure 9 is close to the channel centreline as Wi increases.
When Re is further increased to 100, the elastic effect continues to decrease as the fluid inertia is enhanced, and the equilibrium position and rotational mode of the spheroids are changed significantly. In figure 11, all particles migrate away from the corners, indicating that the fluid inertia overcomes the influence of the first normal stress difference at the corners, and the CO equilibrium position disappears.
Figure 11. Effect of elastic number on the trajectories of the prolate (first line, α = 2) and oblate (second line, α = 0.5) spheroids in the channel cross-section when Re = 100 with (a) Wi = 0.01, El = 0.0001; (b) Wi = 0.5, El = 0.005; (c) Wi = 1, El = 0.01; (d) Wi = 0.01, El = 0.0001; (e) Wi = 1, El = 0.01; (c) Wi = 2, El = 0.02.
When El = 0.0001 (Re = 100, Wi = 0.01), as shown in figure 11(a,d), the first normal stress difference of the flow field is small, and the prolate and oblate spheroids migrate to the CSM equilibrium position. This shows that when the fluid inertia is the dominant effect, the elasto-inertial focusing behaviour is similar to the classic Segre–Silberberg effect in a Newtonian fluid. The first normal stress difference distribution of the flow field is bounded by the particle mass centre. A positive first normal stress difference appears on the particle side close to the wall and a negative first normal stress difference appears on the other particle side close to the channel centre. To clearly show the changes in orientation, only the orientation within t* = 700–1000 is shown at the bottom of figure 12(a), where the value of Py gradually approaches 0 and the prolate spheroid eventually migrate to the CSM equilibrium position with the TU mode. Under the same flow conditions, as shown at the bottom of figure 12(b), the oblate spheroid maintains the LR mode. This is consistent with the rotational mode of prolate and oblate spheroids in inelastic power-law fluids (Hu et al. Reference Hu, Kang, Lin, Lin, Bao and Zhu2023a,Reference Hu, Lin, Lin and Zhub).
Figure 12. Long-time evolution of the orientation for the (a) prolate spheroid in the LR, TUD and TU modes and (b) oblate spheroid in the kayaking, ILR and LR modes, the elastic number is reduced from top to bottom. The particle is placed at a corner with (xin, yin) = (0.395, 0.395) and Re = 100.
When El = 0.005 (Re = 100, Wi = 0.5), as shown in the middle of figure 12(a,b), the prolate spheroid and oblate spheroid maintain the DL equilibrium position, and keep the TUD mode and ILR mode, respectively. When El = 0.01 (Re = 100, Wi = 1), as shown in figure 11(c,e) and the top of figure 12(a,b), the prolate and oblate spheroids also show the same phenomenon as at (Re, Wi) = (10, 0.1) in figure 9(b,e), i.e. the prolate (oblate) spheroid can migrate to the CC (DL) equilibrium position. Due to the increase in Re and Wi, the time for the particles to transform into the stable mode is faster than that in the middle of figure 10(a,b) with the same elastic number. As shown in figure 11(e,f), the oblate spheroid does not migrate to the CC equilibrium position until the elastic number is enhanced as El = 0.02 (Re = 100, Wi = 2).
Here, we examine the dependence of the equilibrium position on the elastic number and Reynolds number in figure 13. For the prolate (α = 2) and oblate (α = 0.5) spheroids, the equilibrium positions when Re = 100, 50 and 10 depend strongly on the elastic number and weakly on the Reynolds number. In addition, the equilibrium position of the spheroid moves closer to the channel centreline as the elastic number increases. That is, a stronger elastic effect means the easier it is for the particles to migrate to the CC equilibrium position. These results are consistent with the elasto-inertial migration of spherical particles in a square channel (Li et al. Reference Li, McKinley and Ardekani2015; Yu et al. Reference Yu, Wang, Lin and Hu2019). Li et al. (Reference Li, McKinley and Ardekani2015) noted that the critical transition elastic number was dependent on the balance between the elastic force and the inertia-induced force, and a similar analysis is conducted here. The elastic force derived by Ho & Leal (Reference Ho and Leal1976) has the following form for two-dimensional Poiseuille flow of an Oldroyd-B fluid:
\begin{equation}{F_e} = \frac{{40}}{3}{\rm \pi} {\rho _f}U_0^2{d^2}\frac{d}{H}El(1 - {\mu _r}){y^\ast },\end{equation}where y* is the dimensionless lateral position in the range of −0.5 to 0.5. The inertia-induced force derived by Ho & Leal (Reference Ho and Leal1974) for the two-dimensional Poiseuille flow is
\begin{equation}{F_i} = C({y^\ast }){\rho _f}U_0^2{d^2}{\left( {\frac{d}{H}} \right)^2},\end{equation}
where C is a function of y* and has a maximum value of approximately 0.24 at y* ≈ 0.15. The balance between the elastic and inertial forces at y* = 0.15 yields the critical elastic number Elc ≈ 0.038(d/H)(1/(1 − 
$\mu$r)), which is approximately 0.0152 for the sphere with d/H = 0.2 and 
$\mu$r = 0.5. The theoretical results of spherical particles are compared in figure 13 with the green dashed line.
Figure 13. Changes in the equilibrium position for the migration of the prolate (α = 2) and oblate (α = 0.5) spheroids at different elastic numbers. The green dashed line shows the theoretically calculated critical elastic number.
As shown in figure 13, when El is less than (greater than) the critical elastic number, the particle shape does not affect the migration of the particles to the DL (CC) equilibrium position. The critical elastic numbers are all near the theoretical value (Elc = 0.0152), and the oblate spheroid remains the same as the theoretical value when Re = 100. At other Reynolds numbers, the critical elastic number changes complexly. As the Reynolds number increases, the critical elastic number in figure 13 shifts slightly to the left, i.e. the critical elastic number for the CC equilibrium position decreases. The curve of the oblate spheroid changes more gently than that of the prolate spheroid; that is, the oblate spheroid can migrate to the CC equilibrium position only under a larger elastic number than that for the prolate spheroid. Near the critical elastic number, the oblate spheroid is more likely to maintain the stable kayaking mode at the DL equilibrium position (in the top of figure 12b), and the effective cross-sectional radius on the diagonal plane is large. While the prolate spheroid transforms from the kayaking mode into the TUD mode (in the middle of figure 12a) at the DL equilibrium position, the effective cross-sectional radius along the diagonal line is smaller than that for the oblate spheroid. Therefore, the critical elastic number of an oblate spheroid is higher than that of a prolate spheroid.
To analyse the change in the equilibrium position of spheroids with different aspect ratios and elastic numbers, we summarize the phase diagram in figure 14 when the spheroid is initially located at a corner (xin, yin) = (0.395, 0.395). As El increases, the equilibrium position of the spheroid gradually changes from the CSM, DL, CC and CO. When the elastic number is greater than the critical value, there is a stable CO equilibrium position. Only when the elastic number is high can the spheroids with 0.4 ≤ α ≤ 0.5 and 2 ≤ α ≤ 3 at the corners keep the CC equilibrium position, as shown in figure 16(b). At high elastic numbers (El = 2), the first normal stress difference near the wall increases (figure 3c,f), and the rotational behaviour of ellipsoids with 0.4 ≤ α ≤ 0.5 and 2 ≤ α ≤ 3 becomes more pronounced, as shown in figure 16(b). The particles are more likely to leave the corner, and the high first normal stress difference near the wall will push the particles towards the low normal stress difference region at the channel centre. Thus, a non-monotonic CC–CO–CC transition of the equilibrium position is observed when 0.4 ≤ α ≤ 0.5 and 2 ≤ α ≤ 3. In addition, the aspect ratio of the spheroid has an impact on the equilibrium position near the critical elastic number. Prolate spheroids with α > 1.5 can stabilize at the CC equilibrium position, while the spheroid maintains the DL equilibrium position when α ≤ 1.5. Therefore, by controlling the elastic number near the critical value, spheroids with large aspect ratio migrate towards the channel centreline, while spheroids with small aspect ratio can stabilize at the corners, which is helpful for the separation of particles with different shapes.
Figure 14. Effect of the aspect ratio and elastic number on the phase diagram of the equilibrium position; here, the initial particle position is (xin, yin) = (0.395, 0.395).
We further compare the effects of the aspect ratio and elastic number on four kinds of equilibrium positions in figure 15.
Figure 15. Effect of aspect ratio and elastic number on the change in the final (a) CSM; (b) DL and CC; and (c) CO equilibrium positions.
The CSM equilibrium position, as shown in figure 15(a), only exists when El is small and Re is large. The CSM equilibrium position of the oblate spheroid increases rapidly, and it eventually stabilizes for the prolate spheroid, especially when α = 2 and 2.25. Since the prolate spheroid maintains the TU mode at the CSM equilibrium position (see figure 17b), it maintains a horizontal orientation for a long time as the aspect ratio increases, resulting in an equilibrium position close to the wall.
By increasing the elastic number, the particles migrate to the DL equilibrium position, which is shown in figure 15(b). Comparing different curves, it is found that a larger elastic number means the closer the particle is to the channel centre. When Re (Wi) maintains the same value, a larger Wi (Re) means the closer (further away) the particle is to the channel centre. When α < 1, the equilibrium position increases as the aspect ratio increases. When α > 1, the equilibrium position decreases as the aspect ratio increases. When the elastic number (El = 0.01) is near the critical value, the oblate spheroid and spheres are stabilized at the DL equilibrium position, while the prolate spheroid cannot keep the DL equilibrium position until the aspect ratio is close to 1. That is, a larger aspect ratio means the easier it is for the prolate spheroid to maintain the CC equilibrium position. Compared with the curves at (Wi, Re) = (0.5, 100), (0.5, 50) and (1, 100), it is found that the role of Re is smaller than that of Wi.
For the CO equilibrium position, as shown in figure 15(c), the change in the equilibrium position can be divided by the value of the spherical particle. The equilibrium position first lowers to the lowest point and then raises upwards, and the CO equilibrium position is close to the wall under a small elastic number. Due to the rotational motion of the spheroids, as shown in figure 8(i), the oblate spheroid with small aspect ratio (α = 0.4, 0.5) migrates near the wall with the LR mode, which is the closest to the wall. The oblate spheroid (0.5 < α < 1) maintains a stable CO equilibrium position with kayaking mode, the length corresponding to the rotational axis of the particles increases as the aspect ratio increases, and the equilibrium position begins to move away from the corners. The prolate spheroid begin to approach the corners with the TUD, VILR and ILR modes with increasing α. The cross-sectional length corresponding to the rotation axis of the prolate spheroid is decreased, the wall effect on the particles is weakened and thus the equilibrium position of the particles is close to the corners.
We further analyse the effect of the elastic number and aspect ratio on the travel distance when particles migrate to the final equilibrium position. When Wi = 1 and Re = 100, the travel distance of the oblate spheroid in figure 16(a) changes little and it cannot migrate to the CC equilibrium position. The fluctuating amplitude of the particle trajectory decreases as the aspect ratio decreases. The prolate spheroid cannot migrate to the CC equilibrium position until the aspect ratio is large enough, and a slenderer prolate spheroid means the faster it migrates to the CC equilibrium position. When El is large (Wi = 2, Re = 1), as shown in figure 16(b), the particles migrate to the CC (first row) or CO (second row) equilibrium positions, and the aspect ratio of the spheroids has an impact. When 0.6 ≤ α ≤ 1.75, the spheroid maintains the CO equilibrium position and the travel distance of the oblate spheroid does not change with the aspect ratio. However, a larger aspect ratio of the prolate spheroid means the greater the travel distance is to the CO equilibrium position. Spheroids with slender (α = 0.4, 0.5) or flatter (α = 2, 3) shapes can migrate towards the CC equilibrium position, and the particle shape does not influence the travel distance.
Figure 16. Effect of the aspect ratio on the change in the travel distance when (Wi, Re, El) = (a) (1, 100, 0.01); (b) (2, 1, 2).
After the particles migrated to the equilibrium position, we further analyse the interaction between the particles and the flow field. In figure 17, the velocity contours along the y-axis and the flow streamlines are plotted on the typical cross-section in front of the particle centre by one sphere diameter or major axis length of the spheroids. In a Newtonian fluid with Wi = 0 and Re = 100, both the prolate spheroid and oblate spheroid can migrate to the CSM equilibrium position, and maintain the TU mode and LR mode, respectively. When the orientation of the prolate spheroid is parallel to the cross-section, its shape is consistent with that of an oblate spheroid. Thus, two recirculation areas induced by the prolate spheroid and oblate spheroid in figure 17(a,f) are similar.
Figure 17. Effect of elastic number on the velocity contours and the streamlines on the cross-section that is ahead of the prolate (first line, α = 2) and oblate (second line, α = 0.5) spheroid centres by a distance of the major axis. The particle moves into the paper. The particle maintains the CSM equilibrium position with (Wi, Re, El) = (a,f) (0, 100, 0), (b,g) (0.01, 100, 0.0001); the CO equilibrium position with (Wi, Re, El) = (c,h) (0.01, 10, 0.001); (i) (0.1, 10, 0.01); ( j) (1, 100, 0.01); and the CC equilibrium position with (Wi, Re, El) = (d) (0.1, 10, 0.01), (e) (1, 100, 0.01). The prolate spheroid maintains the TU, TU, TUD, LR and LR modes from left to right. The oblate spheroid maintains the LR, LR, ILR, kayaking and kayaking modes from left to right.
For the Oldroyd-B viscoelastic fluid, as shown in figure 17(b,g) with (Wi, Re, El) = (0.01, 100, 0.0001), the particles migrate to the CSM equilibrium position. The particle-induced secondary flow is closer to the channel centre, which is quite different from the results in figure 17(a,f). When increasing the elastic number as (Wi, Re, El) = (0.01, 10, 0.001), the prolate spheroid and oblate spheroid in figure 17(c,h) maintain the stable TUD and ILR modes at the DL equilibrium position, respectively. At this time, the particles induce two diagonally symmetrical recirculation areas near the corners, while three corners far from the particles generate symmetrical streamlines along the diagonal line and point towards the particles, prompting the particles to migrate towards the corners. When the elastic number is further increased to the critical values, i.e. (Wi, Re, El) = (0.1, 10, 0.01) (figure 17d,i) or (Wi, Re, El) = (1, 100, 0.01) (figure 17e,j), the prolate spheroid can migrate to the CC equilibrium position and the velocity contour is smaller than that of the prolate spheroid at the DL or CSM equilibrium positions (figure 17a–c). The oblate spheroid migrates to the DL equilibrium position with kayaking mode and induces two recirculation areas near the corners, as shown in figure 17(i,j). Three corners away from the oblate spheroid generate flows directed towards the particles, the streamlines are no longer symmetrical along the diagonal line due to kayaking rotational motion.
By comprehensively comparing figure 17, it is found that the particle-induced migration velocity is mainly affected by Wi, and the Re has a small effect. The inertia effect on the oblate spheroid is greater than that of the prolate spheroid, and a higher particle-induced migration velocity is obtained in figure 17. Therefore, the equilibrium position of the oblate spheroid is more likely to change and the travel distance to the final equilibrium position is shorter than the prolate spheroid.
3.5. Effect of the initial orientation angle
In the above general initial positions, the initial orientation of the spheroid has a negligible effect. Due to the geometric symmetry of the square channel, there are two special initial positions in figure 1, i.e. the wall bisector (xin, yin) = (0, 0.395) and corner (xin, yin) = (0.395, 0.395). When the particles are initially located at the corner or wall bisector, some special initial orientations have an impact on the final equilibrium position and rotational mode. In the present work, the typical initial Euler angles are chosen as n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0), (0, 0, 0), (
${\rm \pi}$/2, 0, 0), (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0) and (
${\rm \pi}$/2, 
${\rm \pi}$/4, 0). The effect of the initial orientation on the rotational behaviour at two special positions (corner and wall bisector) and the general position (xin, yin) = (0.2, 0.395) are analysed in this section. Through numerous numerical simulations, the phase diagram of the rotational behaviours for the prolate spheroid with α = 2 (grey half-circle) and oblate spheroid with α = 0.5 (pink half-circle) is shown in figure 18.
Figure 18. Phase diagram of the rotational behaviours for the prolate spheroid with α = 2 (grey half-circle) and oblate spheroid with α = 0.5 (pink half-circle). The parameters are changed from the inner circle to the outer circle as (Wi, Re, El) = (0.1, 100, 0.001), (0.5, 100, 0.005), (1, 100, 0.01), (1, 10, 0.1) and (1, 1, 1).
Under different parameters, the prolate (oblate) spheroid has six (five) rotational modes. When (Wi, Re, El) = (1, 1, 1), as shown in figure 18, the spheroids mainly maintained the LR mode. Only when the spheroid is placed at the corners will it migrate to the CO equilibrium position, and the initial orientation and particle shape have an impact on the rotational mode. The prolate spheroids (oblate spheroids) maintain the VILR and TUD (LR and ILR) rotational modes at the corners. When the elastic number is near the critical value, such as (Wi, Re, El) = (1, 100, 0.01) and (0.5, 100, 0.005), the rotational modes of the spheroids are complex and depend on the particle shape, initial position and orientation. When the elastic number is further reduced, such as (Wi, Re, El) = (0.1, 100, 0.001), the fluid inertia plays a dominant role, and the types of rotational modes of the spheroids are reduced. The prolate spheroid located at a corner (wall bisector) maintains the TUD (TU) mode, and the initial orientation has no effect. The oblate spheroid at the wall bisector maintains the LR mode, and the initial orientation also has no impact. When the oblate spheroid is initially placed at a corner, it maintains the ILR and LR modes, and the effect of the orientation is small.
In the following section, we further analyse in detail the effect of the initial orientation on the migration of spheroids at two special positions and analyse whether the rotational modes of the spheroids change at four types of equilibrium positions.
3.5.1. CC equilibrium position
When the fluid inertia is low and the elastic number is large (Wi, Re, El) = (1, 1, 1), as shown in figure 19(a), regardless of the initial orientation, both the prolate spheroid and oblate spheroid at the initial position (xin, yin) = (0, 0.395) migrate towards the CC equilibrium position and finally maintain the LR mode, and the initial orientation does not affect the final CC equilibrium position. Thus, the conclusions on the particle equilibrium position in the above sections are applicable to most situations.
Figure 19. Effect of the initial orientation on the migration and rotational behaviours when Re = 1 and Wi = 1. (a) Changes in the trajectory in the ZOY plane with (xin, yin) = (0, 0.395). When the initial orientation is changed from top to bottom, the long-time evolution of the orientation angle for the (b) prolate spheroid and (c) oblate spheroid changes from other rotational modes into the LR mode.
The effect of the initial orientation on the trajectories and rotational modes of the prolate and oblate spheroids is slightly different. When the particles start to migrate, the first normal stress difference is shown in figures 19(b) and 20(a–c). The prolate and oblate spheroids affect the first normal stress difference of the flow field. When the initial Euler angle n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0), as shown at the top of figures 19(b) and 20(b), the prolate spheroid has a weak effect on the first normal stress difference of the flow field. The prolate spheroid always maintains the LR mode without tumbling and migrates the fastest to the CC equilibrium position in figure 19(a). When n = (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0), as shown in figure 20(c), the particles are initially in the kayaking mode, and the travel distance to the final CC equilibrium position in figure 19(a) is the second largest. For n = (0, 0, 0), as shown at the bottom of figures 19(b) and 20(a), since the prolate spheroid is initially in the TU mode, the rotational motion causes the largest travel distance to the CC equilibrium position in figure 19(a). In general, the initial orientation has a small influence on the travel distance of prolate spheroids.
Figure 20. Contour of the first normal stress difference on the cross-section. When the spheroids are released at (xin, yin) = (0, 0.395), the prolate and oblate spheroids maintain the (a,e) TU mode with n = (0, 0, 0), (b,f) LR mode with n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0) and (c,g) kayaking mode with n = (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0). When (xin, yin) = (0.395, 0.395), (d,h) the initial Euler angle n = (0, 0, 0).
For an oblate spheroid, as shown at the bottom of figure 19(a), the trajectories with different orientations almost overlap, and the initial orientations have no effect on the trajectory and the final rotational mode. Combining figure 19(c) with figure 20(e,f,g), the oblate spheroid turns from the TU mode and kayaking mode into LR mode when n = (0, 0, 0), (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0) and (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0). Thus, the particle trajectory and travel distance are not affected by the initial orientations.
3.5.2. CO equilibrium position
When the particles are initially placed at the corner with (Wi, Re, El) = (1, 1, 1), the initial orientation affects the equilibrium position and rotational modes. For the prolate spheroid with different initial orientations, although it migrates to the CO equilibrium position, the motion trajectories and rotational modes are quite different. When n = (0, 0, 0) or (
${\rm \pi}$/2, 0, 0), as shown in figure 20(d), figure 21(a,b) and the bottom of figure 21(c), the prolate spheroid always migrates diagonally with the TUD mode in the XOY plane. The final equilibrium position is closer to the corner than the other initial Euler angles (n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0), (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0), (
${\rm \pi}$/2, 
${\rm \pi}$/4, 0)). For other initial Euler angles (n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0), (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0), (
${\rm \pi}$/2, 
${\rm \pi}$/4, 0)), the prolate spheroid maintains the VILR mode, as shown in figure 21(c), which is consistent with the curve in the middle of figure 4(a).
Figure 21. Effect of the initial orientation on the migration and rotational behaviours when the prolate (first line, α = 2) and oblate (second line, α = 0.5) spheroids are released at (xin, yin) = (0.395, 0.395) with Re = 1, Wi = 1. The trajectories in the channel cross-section with (a) α = 2 and (d) α = 0.5; the corresponding trajectories in the XOZ and YOZ planes for the (b) prolate spheroid and (e) oblate spheroid in the TUD mode. When the initial orientation is changed from top to bottom, long-time evolution of the orientation for the (c) prolate spheroid in the VILR, VILR and TUD modes, ( f) oblate spheroid in the LR, LR and TUD modes.
In figure 21(d), oblate spheroids migrate to the CO equilibrium position, and spheroids with almost orientations eventually keep the LR mode (as shown in figures 6i and 21f). Only when n = (0, 0, 0) or (
${\rm \pi}$/2, 0, 0), as shown in figure 20(h) and at the bottom of figure 21( f), does the oblate spheroid rotate around the long axis and maintain the stable TUD mode on the diagonal line, which is similar to the prolate spheroid in figure 21(a,b). Although the prolate and oblate spheroids all maintain the TUD mode when n = (0, 0, 0) or (
${\rm \pi}$/2, 0, 0), as shown in figure 21(b,e), the oblate spheroid migrates to the CO equilibrium position with shorter travel distance than the prolate spheroid.
3.5.3. DL equilibrium position
As shown in § 3.4, the equilibrium position and rotational mode of the spheroids change when the fluid inertia is increased, especially at the critical elastic number. When the elastic number is near the critical value (Wi, Re, El) = (1, 100, 0.01), the prolate (oblate) spheroid at a corner with different initial orientations can migrate to the CC (DL) equilibrium position and maintain the LR (kayaking) mode, as shown in figures 11(c) and 12(a) (in figures 11e and 12b), and the initial orientation has no effect. When the elastic number is reduced as (Wi, Re, El) = (0.5, 100, 0.005), the prolate spheroid (oblate spheroid) at a corner maintains the TUD mode (ILR mode) at the DL equilibrium position, and the initial orientation also has no effect. Therefore, the results for these two cases are not presented. In this section, we further analyse whether the spheroids with (xin, yin) = (0, 0.395) can migrate to the DL equilibrium position and whether there is an initial orientation in which the spheroids can maintain a stable CSM equilibrium position when (Wi, Re, El) = (0.5, 100, 0.005).
When the particle is initially placed at the wall bisector with (xin, yin) = (0, 0.395), the prolate spheroids with different orientations migrate to the CSM equilibrium position. As shown in figure 22(a), when n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0), (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0) and (0, 0, 0), three rotational modes are found as the LR, kayaking and TU modes, respectively. The oblate spheroid with n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0) keeps the stable LR mode at the CSM equilibrium position. The oblate spheroid with n = (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0) and (0, 0, 0) initially maintains the TU mode, as shown in figures 23(g) and 22(b), it eventually transforms into the LR mode and stabilizes at the CSM equilibrium position. Thus, the effect of the initial orientation on the rotational mode of the oblate spheroid is negligible.
Figure 22. Long-time evolution of the orientation for the (a) prolate spheroid in the LR, kayaking and TU modes and (b) oblate spheroid in the LR mode when (Re, Wi) = (100, 0.5).
Figure 23. Effect of initial orientation on the velocity contours and the streamlines on the cross-section that is ahead of the spheroid centre by the distance of the major axis. The particle moves into the paper. When (Wi, Re, El, xin, yin) = (0.5, 100, 0.005, 0, 0.395), the prolate (first line, α = 2) and oblate (second line, α = 0.5) spheroids initially maintain the (a,f) LR mode, (b,g) TU mode and (c) the prolate spheroid maintains the kayaking mode. When (Wi, Re, El, xin, yin) = (0.5, 100, 0.005, 0.395, 0.395), the prolate (α = 2) and oblate (α = 0.5) spheroids at DL equilibrium position maintain the (d) TUD mode and (h) ILR mode, respectively. When (Wi, Re, El, xin, yin) = (0.1, 100, 0.001, 0.395, 0.395), (e) the prolate spheroid (α = 2) maintains the stable TUD mode, (i) the oblate spheroid (α = 0.5) maintains the stable ILR mode, ( j) the typical process of transitioning from the ILR to LR mode for the oblate spheroid.
As shown in figures 17 and 23, the particle-induced secondary flow on the cross-section is relatively strong as Re increases. When the spheroids are initially released at the wall bisector, as shown in figures 23(a,b) and 23( f,g), both the prolate and oblate spheroids induce obvious secondary flows that are symmetrical along the wall bisector (x = 0), in which the two recirculation areas close to the spheroid are small. Since the interaction between the oblate spheroid (in figure 23f,g with LR mode and TU mode) and the fluid near the wall is more obvious than that of the prolate spheroid, two recirculation areas away from the oblate spheroid are closer to the walls than that of the prolate spheroid. Therefore, the oblate spheroid deviates from the initial rotational mode and eventually maintains the LR mode. When the prolate spheroid maintains a stable kayaking mode at the CSM equilibrium position, as shown in figure 23(c), the particle-induced secondary flow is no longer symmetrical and two obvious recirculation areas are generated on the particle side. As shown in figure 23(d,h) with (xin, yin) = (0.395, 0.395), the obvious recirculation areas along the diagonal line are induced near the short axis of the spheroid, another two new recirculation areas are formed away from the spheroids, and the streamlines point towards the spheroids and push the prolate spheroid (in TUD mode) and oblate (in ILR mode) spheroid to stabilize at the DL equilibrium position.
Overall, for different initial orientations in figure 23(a–d,f–i), the particle-induced lateral migration velocity of the oblate spheroid on the cross-section is greater than that of the prolate spheroid. Therefore, the oblate spheroid migrates to the equilibrium position faster than the prolate spheroid.
3.5.4. CSM equilibrium position
Further reducing the elastic number as (Wi, Re, El) = (0.1, 100, 0.001), the particles initially released from the wall bisector stabilize at the CSM equilibrium position. Prolate (oblate) spheroids with different orientations maintain the TUD (LR) mode, and the initial orientation has no effect. Therefore, those results are not presented in this section. When the spheroid is initially placed at the corner, we mainly analyse the conditions under which spheroids can migrate to the CSM equilibrium position and stabilize at the CO equilibrium position.
The prolate spheroid, as shown in figures 24(a) and 23(e), eventually stabilizes at the DL equilibrium position and maintains the TUD mode when it is released at the corners, and the effect of initial orientation is negligible. Different from the streamlines in figure 23(d), the elastic number at this time is small, the two recirculation areas are positioned away from the diagonal line and the region in which the streamlines point to the corners becomes large. Thus, the effect of particles on the surrounding streamlines is weakened and the particles migrate slowly towards the DL equilibrium position when the elastic number is low.
Figure 24. Effect of the initial orientation on the long-time evolution of the orientation for the (a) prolate spheroid (α = 2) in the TUD mode and (b) oblate spheroid (α = 0.5) in the LR, ILR and ILR modes with (Wi, Re, El, xin, yin) = (0.1, 100, 0.001, 0.395, 0.395).
The oblate spheroid, as shown at the top of figure 24(b), cannot be stabilized at the DL equilibrium position when n = (
${\rm \pi}$/2, 
${\rm \pi}$/2, 0). The typical migration process is shown in figure 23( j), where an asymmetric secondary flow is formed. The particles slowly migrate to the CSM equilibrium position, and the final flow field is shown in figure 17(g), where a symmetrical secondary flow along the wall bisector (x = 0) is formed. For n = (
${\rm \pi}$/4, 
${\rm \pi}$/4, 0) and (0, 0, 0), as shown in the middle and bottom of figure 24(b), the oblate spheroid can maintain a stable DL equilibrium position with ILR mode, and the corresponding flow field is shown in figure 23(i). The oblate spheroid produces a secondary flow that is symmetrical along the diagonal line, which is similar to the prolate spheroid (figure 23e).
4. Summary
The elasto-inertial focusing and rotating characteristics of a spheroid in a square channel flow of Oldroyd-B viscoelastic fluids are studied by the direct forcing/fictitious domain (DF/FD) method. The effect of the fluid elastic number, the aspect ratio, initial orientation angle and position of the particle on the equilibrium position, rotational behaviour and travel distance are explored. Within the present simulated parameters (1 ≤ Re ≤ 100, 0 ≤ Wi ≤ 2, 0.4 ≤ α ≤3), the main conclusions are summarized below.
(1) The results show that there are four kinds of equilibrium positions and six (five) kinds of rotational behaviours for the elasto-inertial focusing of prolate (oblate) spheroids in a square channel flow of Oldroyd-B viscoelastic fluids. The prolate and oblate spheroids both undergo the TU, TUD, LR, ILR and kayaking modes. This work is the first to identify the VILR mode, which is a special case for the migration of the prolate spheroid.
(2) When the fluid inertia is small, the spheroids migrate towards the CC or CO equilibrium positions. While the particles whose initial position far from the corners migrate to the CC equilibrium position, the larger the Weissenberg number is, the faster the particles migrate to the CC equilibrium position. When 0.4 ≤ α ≤ 1 (1 ≤ α ≤ 3), the travel distance decreases (increases) as the aspect ratio increases. With increasing Wi, the maximum first normal stress difference is increased, and the particle shape and the Weissenberg number affect the CO equilibrium position and rotational mode when the particles are initially located near a corner.
(3) For a higher Reynolds number, the effect of the fluid elastic number on the migration and rotational behaviour is complex. The equilibrium position of the spheroid gradually changes from the CO, CC, DL to CSM equilibrium positions as El increases, depending on the particle shape, elastic number and rotational behaviour. The CSM equilibrium position increases nonlinearly with increasing α. Near the critical elastic number, DL and CC equilibrium positions arise depending on the particle shape. The prolate spheroid cannot migrate to the CC equilibrium position until the aspect ratio is large enough, and a slenderer the prolate spheroid means the faster it migrates to the CC equilibrium position. The CO equilibrium position first lowers to the lowest point and then raises with increasing particle aspect ratio.
(4) Only when the particles are initially located at two special locations (corner and wall bisector) do some initial orientations and particle shape have an impact on the final equilibrium position and rotational mode. When the elastic number is reduced to the critical value, the rotational modes of the spheroids are complex and depend on the particle shape, initial position and orientation. When the elastic number is further reduced, the fluid inertia plays a dominant role, the types of rotational modes of the spheroids are reduced.
Overall, abundant phenomena on the elasto-inertial focusing and the rotational behaviour of spheroids in the square channel flow of an Oldroyd-B viscoelastic fluid are found in this work. By controlling the elastic number near the critical value, spheroids with different aspect ratios can be efficiently separated. The results are helpful for quickly predicting the migration of spheroidal particles, which is useful for designing microfluidic devices with high efficiency.
Funding
This work was supported by the National Natural Science Foundation of China with Grant No. 12202392, Major Program of the National Natural Science Foundation of China with Grant No. 12132015, National Natural Science Foundation of China with Grant No. 12372251, China Postdoctoral Science Foundation founded project with Grant No. 2023 M741500, and the Joint Funds of the National Natural Science Foundation of China with Grant No. U2006221.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation
A.1. Motion of a spheroidal particle in Couette flow of Newtonian fluid
A comparison of the present orientation of spheroidal particles in Couette flow with the analytical results is shown in figure 25. The parameters used in the computation are Re = 0.1, Δt = 5×10−4, φ = 
${\rm \pi}$/2, θ = 
${\rm \pi}$/4, α = 2, k = 0.2 and L = H = W = 1. The numerical simulation results are in good agreement with the analytical results.
Figure 25. Comparison of particle orientation between the present result and analytical result (Jeffery Reference Jeffery1922).
A.2. Change in angular velocity of a sphere in an Oldroyd-B fluid
In this work, we further verify the feasibility of the present numerical method by comparing the rotation of a spherical particle in a three-dimensional shear flow of an Oldroyd-B viscoelastic fluid with the numerical and experimental results (Snijkers et al. Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2011). The parameters used in the computation are Re = 0.1, time step Δt = 5 × 10−4, radius of the particle is 0.1 and initial mass centre is (0, 0, 0). The computational domain spans [−L/2, L/2] × [−H/2, H/2] × [−H/2, H/2]. The angular velocities of the particles are shown in figure 26, where the present results are consistent with the previous experiments.
Figure 26. Comparison of angular velocities of spherical particles freely rotating in an Oldroyd-B viscoelastic fluid with different Weissenberg numbers.
