1. Introduction
In viscoelastic fluids, it is well known that suspended particles experience an elastic lift force due to the inhomogeneity of polymer deformation near the particle. Recent experiments evidenced that such disparity in polymer deformation can induce a one-million-fold enhanced lift force for a micrometre-sized spinning particle in viscoelastic fluids (Cao et al. Reference Cao, Das, Windbacher, Ginot, Krüger and Bechinger2023). This pronouncedly affects the particle travelling trajectories. With such elastic lift force in viscoelastic fluids, particles will exhibit fascinating dynamics, such as the cross-flow migration. This phenomenon provides a useful functional operation preceding downstream particulate analysis. For example, in biomedical fields, the viscoelasticity-induced cross-flow migration of particles has been widely utilized for cell sorting, separation and cytometry (D'Avino, Greco & Maffettone Reference D'Avino, Greco and Maffettone2017). Therefore, from the applied perspective, it is crucial to uncover the mechanism of the particle cross-flow migration, which paves the way for the advancement of microfluidic technologies for passive particle manipulation.
Overall, according to the types of lift forces acting on particles, the cross-flow migration of particles can be classified into inertial, viscoelastic and elasto-inertial migration. Different from the inertial migration induced by fluid shear-gradient lift force in Newtonian (Segré & Silberberg Reference Segré and Silberberg1961; Martel & Toner Reference Martel and Toner2014) or shear-thinning (Chrit, Bowie & Alexeev Reference Chrit, Bowie and Alexeev2020; Jannesari Ghomsheh, Jafari & Funfschilling Reference Jannesari Ghomsheh, Jafari and Funfschilling2023) fluids, the viscoelastic or elasto-inertial migration of particles is generated by the elastic lift force. The state of the art has elucidated the fact that viscoelastic and elasto-inertial migrations are indeed related to the non-uniform distribution of the first normal stress difference (
$N_1=\tau _{11}-\tau _{22}$, where 
$\tau _{11}$ and 
$\tau _{22}$ are the diagonal components of the stress tensor along the flow velocity and the velocity gradient directions, respectively) in the flow field. Previous studies have demonstrated that the particles preferentially migrate to the region with a small 
$N_1$ or fluid shear rate (Karnis & Mason Reference Karnis and Mason1966; Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007; Naillon et al. Reference Naillon, de Loubens, Chèvremont, Rouze, Leonetti and Bodiguel2019). For instance, in the microchannel flow of viscoelastic fluids, particles are principally attracted towards the channel centreline, i.e. the inverse Segré–Silberberg scenario, during the viscoelastic migration (Karnis & Mason Reference Karnis and Mason1966; Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007; Villone et al. Reference Villone, D'Avino, Hulsen, Greco and Maffettone2011b; Del Giudice et al. Reference Del Giudice, D'Avino, Greco, Netti and Maffettone2015; Del Giudice Reference Del Giudice2019; D'Avino Reference D'Avino2021) or elasto-inertial migration with highly fluid elasticity (Lim et al. Reference Lim, Ober, Edd, Desai, Neal, Bong, Doyle, McKinley and Toner2014; Seo et al. Reference Seo, Byeon, Huh and Lee2014; Lu & Xuan Reference Lu and Xuan2015). More interestingly, Chaparian et al. (Reference Chaparian, Ardekani, Brandt and Tammisola2020) found that the particle can also focus on the channel centre by entering into the unyielded core in the elastoviscoplastic fluids at the specific ranges of Weissenberg and Bingham numbers.
The viscoelastic lateral migration of particles is governed by several parameters, including channel geometry, suspension rheology and particle characteristics. The exploration of the effects of these parameters on particle cross-flow migration originated initially with the studies on spherical particles in viscoelastic fluids. Regarding the channel geometry, previous results indicate that it affects the particle viscoelastic migration through distinct distributions of 
$N_1$ in viscoelastic channel flow. For ease of fabrication, rectangular microchannels are widely utilized in microfluidic applications (Raoufi et al. Reference Raoufi, Mashhadian, Niazmand, Asadnia, Razmjou and Warkiani2019), which motivates the previous abundant research on particle viscoelastic migration in rectangular channel flow. Leshansky et al. (Reference Leshansky, Bransky, Korin and Dinnar2007) found experimentally that spherical particles suspended in viscoelastic fluids were attracted to the mid-plane of a rectangular microslit with a large aspect ratio, and finally formed a plane focusing, i.e. two-dimensional focusing (D'Avino et al. Reference D'Avino, Greco and Maffettone2017). Later, this two-dimensional focusing pattern has been verified by simulations (Villone et al. Reference Villone, D'Avino, Hulsen, Greco and Maffettone2011a). Furthermore, Seo et al. (Reference Seo, Byeon, Huh and Lee2014) and D'Avino et al. (Reference D'Avino, Romeo, Villone, Greco, Netti and Maffettone2012, Reference D'Avino, Greco and Maffettone2017) also found that the migration direction of the particle in viscoelastic fluids is dependent on the its initial positions. The aspect ratio of the channel obviously affects the elastic stress distribution in the channel. With the aspect ratio of the channel decreasing, the regions with small 
$N_1$ are located in the channel centre and the corner of the channel cross-section. This corresponds to the multi-focusing region of particles in the flow field (Del Giudice et al. Reference Del Giudice, D'Avino, Greco, Netti and Maffettone2015). In another widely used microchannel, i.e. the tube, particles have been observed to migrate towards the pipe centreline, forming a ‘line-focusing’ pattern (so-called three-dimensional focusing). Additionally, the phenomenon of particle focusing propelled by fluid elasticity has recently been extended to the viscoelastic channel flow with more complex geometries, including rhomb (Kwon et al. Reference Kwon, Kim, Kim and Cho2020), semi-ellipse (Tang et al. Reference Tang, Fan, Yang, Li, Zhu, Jiang, Shi and Xiang2019) and triangle (D'Avino Reference D'Avino2021). With these meticulously designed microchannels, the particles form a desired focusing pattern.
On the other hand, numerous polymeric fluids encountered in nature and realistic applications generally exhibit viscoelasticity as well as shear-thinning rheology. This fact prompted previous studies investigating the effect of suspension rheology on particle migration. For example, the viscoelastic fluid with a non-zero second normal stress difference (
$N_2$) can induce a secondary flow (Villone et al. Reference Villone, D'Avino, Hulsen, Greco and Maffettone2013) within the channel with a non-circle cross-section. Such secondary flow promotes particles to spirally approach the equilibrium positions in the square channel filled with Giesekus fluids (Villone et al. Reference Villone, D'Avino, Hulsen, Greco and Maffettone2013). Del Giudice et al. (Reference Del Giudice, D'Avino, Greco, Netti and Maffettone2015) also reported experimentally that particles can be driven towards the corners of the channel cross-section in highly shear-thinning viscoelastic fluids. Apart from particle distribution in the channel, the simulation results highlighted that the particle equilibrium position is also modulated by shear-thinning rheology. Specifically, the shear-thinning effect tends to push particles away from the channel centreline (Li, McKinley & Ardekani Reference Li, McKinley and Ardekani2015). In theoretical studies, Wang, Tai & Narsimhan (Reference Wang, Tai and Narsimhan2020) examined the effect of the normal stress on the elastic migration of spheroids in the second-order fluid. The above results inspired us to use suspension rheology and channel geometry for manipulating the focusing pattern of flowing particles in microfluidic applications.
Besides the channel geometry and suspension rheology, another prominent parameter governing the particle cross-flow migration is the particle shape. Several studies have already evidenced that the particle cross-flow migration can be greatly changed due to the effect of particle shape. For instance, Lu & Xuan (Reference Lu and Xuan2015) and Lu et al. (Reference Lu, Zhu, Hua and Xuan2015) reported a pinched flow fractionation for continuous shape-based particle separation. In their experiments, the equilibrium positions of spherical and peanut-shaped rigid particles are different. More importantly, for a non-spherical particle, an interesting phenomenon is that the spheroid can exhibit multi-orientation modes in viscoelastic shear flow with different fluid elasticities (D'Avino & Maffettone Reference D'Avino and Maffettone2015). Such orientation dynamics of spheroids will introduce additional complexities into the particle cross-flow migration process, thereby rendering the cross-flow migration of non-spherical particles more intricate. For this point, D'Avino et al. (Reference D'Avino, Hulsen, Greco and Maffettone2019) investigated numerically the viscoelastic migration of a single prolate spheroid in a wide-slit microchannel. Similar to spherical particles, most prolate spheroids released from different initial positions migrate to the channel centre plane and tend to orientate in the streamwise direction at the equilibrium position (Langella et al. Reference Langella, Franzino, Maffettone, Larobina and D'Avino2023). More recently, from the theoretical perspective, Tai & Narsimhan (Reference Tai and Narsimhan2022) derived a cross-flow migration model for spheroids in second-order fluid. Their prediction model suggests that the particle migration velocity is proportional to the length of particles projected in the velocity gradient direction (Tai, Wang & Narsimhan Reference Tai, Wang and Narsimhan2020a,Reference Tai, Wang and Narsimhanb). Consequently, the non-spherical particles exhibit slower migration compared to spheres.
Finally, as for the role of fluid inertia, unlike the inertial lateral migration of particles, the viscoelastic lateral migration of particles is attributed to the nonlinear nature of the constitutive equation of viscoelastic fluids. Thus the viscoelastic lateral migration does not rely solely on the fluid inertia, and can appear in viscoelastic flow even with the vanishing fluid inertia. Such facts motivate that most previous studies are often limited to small-Reynolds-number (
$Re\ll 1$) flow, obviously affecting the overall throughput of particle delivery in the flow system (D'Avino et al. Reference D'Avino, Greco and Maffettone2017). Increasing the fluid inertia in the flow system can achieve high-throughput particle focusing in viscoelastic fluid flows (Lim et al. Reference Lim, Ober, Edd, Desai, Neal, Bong, Doyle, McKinley and Toner2014), thus several researchers explored the combining effects of fluid inertia and elasticity on the particle focusing in elasto-inertial flows. For instance, Lim et al. (Reference Lim, Ober, Edd, Desai, Neal, Bong, Doyle, McKinley and Toner2014) reported experimentally on the deterministic focusing of particles in a microchannel at extremely high flow rates, corresponding to 
$Re \sim 10\,000$; they found that the particle viscoelastic focusing can also occur over a wide range of Reynolds numbers. Holzner, Stavrakis & de Mello (Reference Holzner, Stavrakis and de Mello2017) explored the elasto-inertial focusing of bio-particles within the straight microfluidic channel flow of viscoelastic PEO solutions with low viscosities. Li et al. (Reference Li, McKinley and Ardekani2015) simulated the elasto-inertial migration of a single spherical particle in the viscoelastic channel flow with finite fluid inertia. Raffiee, Ardekani & Dabiri (Reference Raffiee, Ardekani and Dabiri2019) investigated the particle elasto-inertial focusing patterns in viscoelastic channel flow. These numerical studies elaborated that the particle equilibrium position is indeed determined by the competition between fluid inertia and fluid elasticity. More importantly, from these results, tuning the finite fluid inertia and fluid elasticity appears to be a prevalent and feasible method for achieving an ‘optimal focusing’ for particles in microfluidic devices (D'Avino et al. Reference D'Avino, Greco and Maffettone2017; Raoufi et al. Reference Raoufi, Mashhadian, Niazmand, Asadnia, Razmjou and Warkiani2019). This also motivates the present study.
In summary, earlier studies on the elasto-inertial migration of particles are confined mainly to spherical particles. The dynamics of the non-spherical particle during its elasto-inertial migration is still not well understood. The main obstacle to further progress is the lack of understanding of the physical picture in the elasto-inertial migration of spheroids. How the spheroid orientation evolves, and how the spheroid elasto-inertial migration is influenced by its orientation, are the key questions to be answered. The present study is motivated by these questions and carried out by fully resolving the elasto-inertial migration of a single neutrally buoyant oblate spheroid suspended in viscoelastic duct flow. The focus on the oblate shape in the present study is motivated by the fact that the oblate spheroid is a canonical shape for bio-particles (such as red blood cells) in biomedical applications (Atwell et al. Reference Atwell, Badens, Charrier, Helfer and Viallat2022; Jiang, Liu & Tang Reference Jiang, Liu and Tang2022), where microfluidic technologies are employed extensively.
The paper is organized as follows. In § 2, we introduce the mathematical models and numerical methods of fully resolved simulations of the particle migration in the elasto-inertial duct flows. Then the elasto-inertial migration of the oblate particle is discussed in § 3, including the ‘oscillating’ migration (§§ 3.1 and 3.2) and rotational dynamics (§ 3.3) of the oblate particle. Finally, conclusions are drawn in § 4.
2. Methodology
In this section, a concise but complete summary of the simulation methodology is provided. The elasto-inertial flow of Oldroyd-B fluid is simulated within a monolithic projection framework (Li et al. Reference Li, Huang, Xu and Zhao2022). The particle migration is tracked in the Lagrangian frame. The interaction between fluid flow and particles is fully resolved by the immersed boundary method (IBM).
2.1. Governing equations
2.1.1. Elasto-inertial flow of viscoelastic fluid
The governing equations for incompressible and isothermal viscoelastic fluid flow, accounting for finite fluid inertia, are expressed as
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather} {\rho _f}\left( {\frac{{\partial \boldsymbol{u}}}{{\partial t}} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}} \right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}_s + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}_p + \boldsymbol{f}_{IB}, \end{gather} where 
$\boldsymbol {u}$ is the fluid velocity, 
${\rho _f}$ is the fluid density, 
$p$ is the pressure, 
$\boldsymbol {\tau }_s$ is the solvent stress, 
$\boldsymbol {\tau }_p$ is the polymer stress, and 
$\boldsymbol {f}_{IB}$ is the momentum forcing due to the fluid–particle interaction, representing the feedback forcing from the particle to fluid phase.
Generally, the solvent stress 
$\boldsymbol {\tau }_s$ satisfies the constitutive equation of Newtonian fluid:
\begin{equation} {\boldsymbol{\tau}_s} = \beta {\mu _0}\left( {\boldsymbol{\nabla} \boldsymbol{u} + {{\left( {\boldsymbol{\nabla} \boldsymbol{u}} \right)}^{\rm T}}} \right),\quad {\mu _0} \equiv {\mu _s} + {\mu _p}, \quad \beta \equiv \frac{{{\mu _s}}}{{{\mu _0}}}, \end{equation}
where 
${\mu _0}$ denotes the zero-shear-rate dynamic viscosity, 
${\mu _s}$ and 
${\mu _p}$ represent solvent and polymeric viscosity, respectively, and 
$\beta$ is the viscosity ratio.
Based on the kinetic theory, the polymer stress 
$\boldsymbol {\tau }_p$ in Oldroyd-B fluid can be determined by the conformation tensor 
$\boldsymbol{\mathsf{B}}$, which quantifies the deformation of the polymeric molecule in the flow field:
\begin{equation} {\boldsymbol{\tau}_p} = \frac{{{\mu _p}}}{\lambda }\,( {\boldsymbol{\mathsf{B}} - \boldsymbol{\mathsf{I}}} ), \end{equation}
where 
$\lambda$ is the polymer relaxation time. Then the evolution of the conformation tensor 
$\boldsymbol{\mathsf{B}}$ can be written as
\begin{equation} \frac{{\partial \boldsymbol{\mathsf{B}}}}{{\partial t}} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\mathsf{B}} - \boldsymbol{\mathsf{B}} \boldsymbol{\cdot} \left( {\boldsymbol{\nabla} \boldsymbol{u}} \right) - {\left( {\boldsymbol{\nabla} \boldsymbol{u}} \right)^{\rm T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{B}} = \frac{1}{\lambda }\,( {\boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{B}}} ). \end{equation}From the above governing equations, the nonlinearity in the elasto-inertial flow of viscoelastic fluids arises from the fluid convection in (2.2) and the nonlinear deformation of polymeric molecules described in (2.5).
2.1.2. Particle dynamics
To determine the intricate translation and rotation of an oblate spheroid induced by the interaction between the fluid and particle, the following Newton–Euler equation is employed to model particle motion:
\begin{gather} {m_p}\,\frac{{{{\rm d}}{{\boldsymbol{u}}_p}}}{{{\mathrm{d}}t}} = \oint_{{\varGamma _p}} {{\boldsymbol{\tau }} \boldsymbol{\cdot} {\boldsymbol{n}}\,{{\rm d}}s} + {{\mathcal{V}}_p}( {{\rho _p} - {\rho _f}} ){\boldsymbol{g}}, \end{gather}
\begin{gather} \frac{{{\rm d}( {{\boldsymbol{I}_p}{\boldsymbol{\omega}_p}} )}}{{{\rm d}t}} = \oint_{{\varGamma _p}} {\boldsymbol{r} \times \left( {\boldsymbol{\tau}\boldsymbol{\cdot} \boldsymbol{n}} \right)\,{{\rm d}}s}, \end{gather}
where 
$\boldsymbol {u}_p$ and 
$\boldsymbol {\omega }_p$ represent the particle translational and rotational velocities, respectively, 
$m_p$ is the particle mass, 
$\boldsymbol {I}_p$ is the particle moment of the inertial matrix, 
${\mathcal {V}}_p$ denotes the particle volume enclosed by particle surface 
$\varGamma _p$, 
$\rho _p$ is the particle density, 
$\boldsymbol {\tau }$ represents the hydrodynamic stress tensor, whose integration accounts for the fluid-particle interaction, 
$\boldsymbol {n}$ denotes the unit normal vector pointing outwards on the particle surface, 
$\boldsymbol {r}$ is the position vector on the particle surface from the particle's centre, and 
$\boldsymbol {g}$ is the acceleration due to gravity.
To accurately resolve the interplay between particle and viscoelastic flow, the IBM is adopted in the present study. Within the IBM framework, the hydrodynamic force and torque acting on the particle are rewritten as
\begin{gather} \oint_{{\varGamma _p}} {\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n} \,\mathrm{d}s} = \frac{\mathrm{d}}{{\mathrm{d}t}}\int_{{\varOmega _p}} {{\rho _f}\boldsymbol{u}}\, \mathrm{d}v - \int_{{\varOmega_p}} {{\boldsymbol{f}_{IB}}} \,\mathrm{d}v, \end{gather}
\begin{gather} \oint_{{\varGamma _p}} {\boldsymbol{r} \times \left( {\boldsymbol{\tau}\boldsymbol{\cdot} \boldsymbol{n}} \right)\,\mathrm{d}s} = \frac{\mathrm{d}}{{\mathrm{d}t}}\int_{{\varOmega _p}} {{\rho _f}\boldsymbol{r} \times \boldsymbol{u}} \,\mathrm{d}v - \int_{{\varOmega _p}} {\boldsymbol{r} \times {\boldsymbol{f}_{IB}}} \,\mathrm{d}v, \end{gather}
where 
${\varOmega _p}$ is the particle region bounded by surface 
${\varGamma _p}$. With (2.8) and (2.9), the governing equation of particle dynamics can be further approximated as
\begin{gather} {m_p}\,\frac{{{{\rm d}}{{\boldsymbol{u}}_p}}}{{{{\rm d}}t}} \approx{-} \sum_l^{{N_l}} {{\boldsymbol{F}_{IB,l}}}\,\Delta {s_l} + {\rho _f}\,\frac{{{\rm d}}}{{{{\rm d}}t}}\left( {\int_{{\varOmega _p}} \boldsymbol{u}\, {{\rm d}}v} \right) + {{\mathcal{V}}_p}( {{\rho _p} - {\rho _f}} ){\boldsymbol{g}}, \end{gather}
\begin{gather} \frac{{\mathrm{d}( {{\boldsymbol{I}_p}{\boldsymbol{\omega}_p}} )}}{{\mathrm{d}t}} \approx{-} \sum_l^{{N_l}} {\boldsymbol{r} \times {\boldsymbol{F}_{IB,l}}}\,\Delta {s_l} + {\rho _f}\,\frac{\mathrm{d}}{{\mathrm{d}t}}\left( {\int_{{\varOmega _p}} {\boldsymbol{r} \times \boldsymbol{u}} \,\mathrm{d}v} \right), \end{gather}
where 
$\Delta {s_l}$ is the surface area of each Lagrangian element on the particle surface, and 
${N_l}$ is the total number of Lagrangian points.
In the penalty IBM, the fluid–particle interaction forcing 
$\boldsymbol {F}_{IB}$ is determined by (Huang, Chang & Sung Reference Huang, Chang and Sung2011)
\begin{equation} {\boldsymbol{F}_{IB}} ={-} \kappa \left[ {( {{\boldsymbol{X}_{IB}} - \boldsymbol{X}} ) + \Delta t\,( {{\boldsymbol{U}_{IB}} - \boldsymbol{U}} )} \right], \end{equation}
where 
$\kappa$ is a large penalty constant in the penalty IBM (
$\kappa = 10^4$ in the present simulations), and 
${\boldsymbol {X}_{IB}}$ and 
${\boldsymbol {U}_{IB}}$, i.e. the position and velocity of the massless counterpart of the material Lagrangian points, respectively, are computed by
\begin{equation} {\boldsymbol{U}_{IB}} = \int_{{\varOmega _p}} \boldsymbol{u}\,\delta ( {\boldsymbol{X} - \boldsymbol{x}} ) \,\mathrm{d}v \end{equation}and
\begin{equation} {\boldsymbol{X}_{IB}}( t ) = {\boldsymbol{X}_{IB}}( {t = 0} ) + \int_0^t {{\boldsymbol{U}_{IB}}( t )\,\mathrm{d}t}, \end{equation}
where 
$\delta (\cdot )$ is the Dirac delta function.
Correspondingly, 
$\boldsymbol {X}$ and 
$\boldsymbol {U}$ are the position and velocity of the material Lagrangian points on the particle surface, respectively. They are calculated by (2.6) and (2.7), 
$\boldsymbol {U} = \boldsymbol {u}_p + \boldsymbol {\omega }_p \times \boldsymbol {r}$. Finally, the momentum forcing 
$\boldsymbol {f}_{IB}$ in (2.2) is spread from 
$\boldsymbol {F}_{IB}$ as
\begin{equation} {\boldsymbol{f}_{IB}} = \int_{{\varGamma _p}} {{\boldsymbol{F}_{IB}}}\,\delta ( {\boldsymbol{x} - \boldsymbol{X}} ) \,\mathrm{d}s. \end{equation}2.2. Numerical methods
The governing equations of viscoelastic fluid flow are discretized on the staggered grid using the finite difference method. For temporal discretization of the momentum and constitutive equations, all terms are integrated in time by the second-order Crank–Nicolson scheme. To retain the second-order accuracy in time, the incremental form is used to treat the pressure (
$\delta p$ form) and conformation tensor (
$\delta \boldsymbol {B}$ form) with the staggered time scheme. For spatial discretization, all terms are approximated by the second-order central difference scheme, except for the convective term in the constitutive equation, which is evaluated by the high-resolution MINMOD scheme (Alves, Pinho & Oliveira Reference Alves, Pinho and Oliveira2000). By rearranging the discretized governing equations into a monolithic matrix system
\begin{equation} \left[ {\begin{array}{ccc} {{{\mathcal{A}}_B}} & 0 & 0\\ { - {{\mathcal{D}}_B}} & {{{\mathcal{A}}_u}} & {\mathcal{G}}\\ 0 & {\mathcal{D}} & 0 \end{array}} \right]\left[ {\begin{array}{c} {\delta \boldsymbol{B}}\\ {\delta \boldsymbol{u}}\\ {\delta p} \end{array}} \right] = \left[ {\begin{array}{c} {\boldsymbol{r}_B^{n - {1/2}} - {{\mathcal{A}}_B}{\boldsymbol{B}^{n - {1/2}}}}\\ {{\boldsymbol{r}^n} - {{\mathcal{A}}_u}{\boldsymbol{u}^n}}\\ 0 \end{array}} \right] + \left[ {\begin{array}{c} 0\\ {\boldsymbol{mbc}{^{n + 1}}}\\ {cb{c^{n + 1}}} \end{array}} \right], \end{equation}the pressure, conformation tensor and velocity can be decoupled sequentially from the viscoelastic flow system based on the approximate factorization of the system coefficient matrix, as depicted in
\begin{align} &\left[
{\begin{array}{ccc} {{{\mathcal{A}}_B}} & 0 & 0\\ { -
{{\mathcal{D}}_B}} & {{{\mathcal{A}}_u}} & {\mathcal{G}}\\
0 & {\mathcal{D}} & { - \Delta
t\,{\mathcal{D}}{\mathcal{G}}}
\end{array}} \right]\left[ {\begin{array}{ccc} I & 0 & 0\\
0 & I & {\Delta t\,{\mathcal{G}}}\\ 0 &
0 & I \end{array}} \right]\left[ {\begin{array}{c} {\delta
\boldsymbol{B}}\\ {\delta \boldsymbol{u}}\\ {\delta p}
\end{array}} \right]\nonumber\\ &\quad = \left[ {\begin{array}{c}
{\boldsymbol{r}_B^{n - {1/2}} -
{{\mathcal{A}}_B}{\boldsymbol{B}^{n - {1/2}}}}\\
{{\boldsymbol{r}^n} -
{{\mathcal{A}}_u}{\boldsymbol{u}^n}}\\ 0 \end{array}}
\right] + \left[ {\begin{array}{c} 0\\
{\boldsymbol{mbc}{^{n + 1}}}\\ {cbc{^{n + 1}}} \end{array}}
\right],\end{align}
where 
${\mathcal {A}}_u$ and 
${\mathcal {A}}_B$ denote the coefficient matrix of the discretized governing equations of 
$\boldsymbol {u}^{n+1}$ and 
$\boldsymbol {B}^{n+1/2}$, respectively, 
$\mathcal {D}$ and 
$\mathcal {D}_B$ are the discrete divergence operators of 
$\boldsymbol {u}^{n+1}$ and 
$\boldsymbol {B}^{n+1/2}$, respectively, 
$\mathcal {G}$ represents the discrete gradient operator of pressure, 
$\boldsymbol {r}$ and 
$\boldsymbol {r}_B$ are the right-hand sides of discretized governing equations of 
$\boldsymbol {u}^{n+1}$ and 
$\boldsymbol {B}^{n+1/2}$, respectively, and 
$\boldsymbol {mbc}$ and 
$cbc$ are the discretizations of the boundary conditions in the momentum and continuity equations.
Then all flow quantities can be resolved without iteration within the present projection framework. The detailed numerical method can be found in Li et al. (Reference Li, Huang, Xu and Zhao2022).
For the particle phase, by employing the fourth-order Runge–Kutta scheme, the particle translation is tracked in the Eulerian inertial frame (
$x\unicode{x2013}y\unicode{x2013}z$), while particle rotation and orientation are simultaneously solved within the particle frame (
$x'\unicode{x2013}y'\unicode{x2013}z'$), in which the origin of the particle frame is in the particle centre, and the coordinate axes are aligned with the principal directions of the spheroid particle. Moreover, the rotation and orientation modes of the spheroids are presented in the particle co-moving frame (
$x''\unicode{x2013}y''\unicode{x2013}z''$) whose coordinate axes are parallel to those of the Eulerian inertial frame. The origin of the particle co-moving frame is attached at the mass centre of the particle. During particle migration, the particle orientation is represented by the quaternions (Goldstein Reference Goldstein1980), which are updated according to the particle angular velocity. The detailed validation of the present numerical methods can be found in Appendix A.
2.3. Simulation set-up
2.3.1. Computational domain and particle size
In the present study, an individual oblate spheroid is suspended in a duct flow of viscoelastic fluid, as illustrated in figure 1. The shear-thinning viscoelastic fluid (modelled by Giesekus fluid), whose 
$N_2$ is not null, can introduce a secondary flow in the duct cross-section (Villone et al. Reference Villone, D'Avino, Hulsen, Greco and Maffettone2013). Similar to the previous work (Yu et al. Reference Yu, Wang, Lin and Hu2019), the rheology of the viscoelastic fluid in this study is described by the Oldroyd-B model with 
$\beta = 0.1$ to eliminate the effect of the secondary flow on the spheroid migration. The size of the duct is set as 
$W \times H \times L = 4D \times 4D \times 8D$, where 
$L$ is the duct length in the streamwise direction (
$z$-axis), 
$H$ and 
$W$ are the height and width of the duct, respectively, and 
$D$ is the major diameter of the oblate spheroid, 
$D = 1$. Such a computational size in the streamwise direction of the microchannel 
$L/H = 2$ is also utilized in previous studies investigating the elasto-inertial migration of a sphere in Oldroyd-B fluid (Yu et al. Reference Yu, Wang, Lin and Hu2019). The effect of the computational size is examined in Appendix B. The time step is set as 
$\Delta t = 1\times 10^{-3}$, unless otherwise stated. The mesh resolution in this study is set as 
$\varDelta = 1/32D$, indicating that the major diameter is resolved by 
$32$ grids. The mesh convergence test is also shown in Appendix B.
Sketch of the flow configuration.
The confinement in the flow system of particle cross-flow migration is defined as 
$D/H$, which generally takes values 
$D/H = 0.1\unicode{x2013}0.5$ for inertial migration in Newtonian fluid (Hafemann & Fröhlich Reference Hafemann and Fröhlich2023) and elastic migration in viscoelastic fluid (Holzner et al. Reference Holzner, Stavrakis and de Mello2017; D'Avino et al. Reference D'Avino, Hulsen, Greco and Maffettone2019; Naillon et al. Reference Naillon, de Loubens, Chèvremont, Rouze, Leonetti and Bodiguel2019; Raffiee et al. Reference Raffiee, Ardekani and Dabiri2019; Raoufi et al. Reference Raoufi, Mashhadian, Niazmand, Asadnia, Razmjou and Warkiani2019; Zhou & Papautsky Reference Zhou and Papautsky2020) in previous studies. Considering the above typical confinement, the present study focuses on the confinement as 
$D/H = 0.25$.
2.3.2. Non-dimensionalization of the simulation results
In the present flow system, the characteristic velocity and length are set as the fluid velocity 
$U_c$ at the duct centreline and the duct height 
$H$, respectively. The characteristic density is fluid density 
${\rho _f}$. Thus the non-dimensionalization of the simulation results in the present study is realized as follows.
For fluid phase,
\begin{equation} {\boldsymbol{x}^*} = \frac{\boldsymbol{x}}{H},\quad {\boldsymbol{u}^*} = \frac{\boldsymbol{u}}{{{U_c}}},\quad {t^*} = \frac{t}{{{H/{{U_c}}}}},\quad {p^*} = \frac{p}{{\rho U_c^2}}. \end{equation}For particle phase,
\begin{equation} \boldsymbol{X}_p^* = \frac{{{\boldsymbol{X}_p}}}{H},\quad \boldsymbol{V}_p^* = \frac{{{\boldsymbol{V}_p}}}{{{U_c}}},\quad \boldsymbol{\omega} _p^* = \frac{{{\boldsymbol{\omega}_p}}}{{{{{U_c}} / H}}},\quad {D^*} = \frac{D}{H}, \end{equation}
where 
${\boldsymbol {X}_p}=({X}_p,{Y}_p,{Z}_p)$ and 
${\boldsymbol {V}_p}=({U}_p,{V}_p,{W}_p)$. The superscript 
$*$ represents the dimensionless quantities; for conciseness, the superscript 
$*$ is dropped unless explicitly stated otherwise.
Based on the above characteristic quantities, the dimensionless governing parameters in the present flow configuration are defined as follows.
The Reynolds number 
$Re$, which quantifies the fluid inertial effect, is defined as
\begin{equation} Re = U_cH/\nu_0, \end{equation}
in which 
$\nu _0$ is the zero-shear-rate kinematic viscosity 
$\nu _0 = \mu _0/{\rho _f}$.
The Weissenberg number 
$Wi$, which represents the fluid elastic effect, reads
\begin{equation} Wi = \lambda U_c/H. \end{equation} The elastic number 
$El$, representing the competition between fluid elastic and inertial effect, is defined as
\begin{equation} El = Wi /Re. \end{equation}For the particle phase, the particle shape is described by the particle aspect ratio
\begin{equation} AR = a/b, \end{equation}
where 
$a$ and 
$b$ are the polar and equatorial radii of the spheroid, respectively. As the state of the art highlighted in § 1, the features of particle viscoelastic and elasto-inertial migrations have been analysed extensively for both spheres (
$AR = 1$) and prolate spheroids (
$AR >1$). However, only a few works focus on the viscoelastic or inertial migrations of oblate spheroids (
$AR <1$) (Nizkaya et al. Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020; Tai & Narsimhan Reference Tai and Narsimhan2022). Thus to conveniently elucidate the unique behaviour of oblate spheroids in elasto-inertial flows, in this study, we focus on the oblate spheroid with the same aspect ratio (
$AR = 0.5$) used in earlier studies concerning the inertial migration of oblate spheroids (Nizkaya et al. Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020).
The inertia of a particle is characterized by the particle size and density ratio:
\begin{equation} \rho_r = \rho_p / \rho_f. \end{equation} In the present work, we focus on a neutrally buoyant particle, i.e. 
$\rho _r = 1$, which represents the typical density of bio-particles commonly manipulated by microfluidic technologies (Atwell et al. Reference Atwell, Badens, Charrier, Helfer and Viallat2022).
2.3.3. Initial and boundary conditions
For the initial conditions, a rest spheroid with zero translational and rotational velocities is initially released from 
$(X_{p0}, Y_{p0}, Z_{p0}) = (0, 0.25, 0)$. The initial orientation of the oblate particle is set by the Euler angle (
$\phi _0, \theta _0, \psi _0$) = (
$-1/4{\rm \pi}, 1/2{\rm \pi}, 0$). For the viscoelastic fluid phase, to eliminate the transient process during the start-up flow of viscoelastic fluid, the initial flow is a fully developed viscoelastic duct flow driven by a constant 
${\rm d}p/{\rm d}z = -U_c{\rm \pi} ^3\mu _0/(4kH^2)$, and the parameter 
$k$ takes value 
$k\approx 0.571$ for a square-shaped duct flow (Li et al. Reference Li, McKinley and Ardekani2015; Yu et al. Reference Yu, Wang, Lin and Hu2019).
For the boundary conditions, the periodic boundary is set in the streamwise (
$z$) direction of the duct flow, and the no-slip boundary is applied at the duct walls, which are normal to the 
$x$ and 
$y$ directions. The particle–fluid interface is set as a no-slip boundary realized by the IB forcing.
3. Results and discussion
In this section, we analyse the coupling effect of fluid inertia and fluid elasticity on the lateral migration (§ 3.1) and rotation modes (§ 3.2) of the neutrally buoyant oblate spheroid (
$\rho _r = 1$, 
$AR = 0.5$) in a viscoelastic duct flow driven by a constant pressure gradient 
${\rm d}p/{\rm d}z = -3.6201$. The above-mentioned governing parameters are set as 
$Re = 5 \unicode{x2013} 60$, 
$El = 0 \unicode{x2013} 0.1$.
3.1. Oscillating migration
First, the lateral migration of the oblate spheroid in the viscoelastic duct flow with different 
$El$ is depicted in figure 2. In figure 2(a), for the inertial migration of particles in Newtonian fluids (
$El = 0$), with the inertial lift force, the oblate spheroid gradually moves towards the wall, and ultimately approaches its equilibrium position close to the wall. However, the scenario becomes different in viscoelastic flows (
$El\ne 0$). The elastic lift force alters the particle's migration direction, causing it to be attracted to the duct centreline. The specific focusing position of the particle is governed by the elastic number 
$El$: the larger the fluid elasticity is, the closer the particle equilibrium position is to the duct centreline. This observed spheroid focusing behaviour is similar to that of spherical particles (Li et al. Reference Li, McKinley and Ardekani2015). Note that 
$Y_{peq}$ in the elasto-inertial flow at 
$El=0.01$ is a quasi-equilibrium position due to the perturbations induced by the particle rotation. The process of leaving out the duct symmetry plane for the spheroid is sufficiently long: the particle just moves 
$\sim O(10^{-3})\,H$ in the 
$x$-axis direction after it translates 
$\sim O(10^3)\,H$ along the streamwise direction. This indicates that the equilibrium position in the symmetry plane of the duct is unstable or meta-stable in the flow at 
$El=0.01$.
Lateral migration of the oblate spheroid in elasto-inertial duct flows. (a) Particle position 
$Y_p$ in the wall-normal direction. (b) Particle migration velocity 
$V_p$. The black dashed line represents the particle initial position. The dots of different colours in (b) denote the final positions of particles in the flow at the corresponding 
$El$.
However, different from spherical particles, figure 2(a) further reveals a prominent phenomenon: the time for the oblate particle approaching the duct centreline changes non-monotonically with 
$El$. For instance, as the fluid elasticity increases from 
$El = 0.03$ to 
$0.05$, the augmented elastic lift force exerted on the particle results in the faster particle focusing on the duct centreline. However, if the fluid elasticity is further increased to 
$El = 0.1$, then the migration of the oblate particle to its destination is retarded. Alternatively, the particle focusing period in highly elastic flows (
$El = 0.1$) is longer compared to that in moderately elastic flows, i.e. 
$El =0.03$ and 
$0.05$.
In summary, the above results suggest the existence of three critical elastic numbers within the present flow system: 
$El_{cr1} < El_{cr2} < El_{cr3}$, as illustrated in figure 3. Specifically, the first elastic number 
$El_{cr1}$ determines the direction of particle migration. When 
$El < El_{cr1}$, the fluid inertia is dominant in the flow system, causing the particle to move towards the duct wall. When 
$El > El_{cr1}$, the particle migration direction is reversed by the fluid elasticity, and the particle begins to migrate to the duct centreline. The second elastic number 
$El_{cr2}$ represents the critical fluid elasticity, at which the fluid elastic effect fully dominates in the flow system. The destination of particles occurs at the duct centreline in the flow at 
$El_{cr2}$. These two critical elastic numbers span a specific range of fluid elasticity, 
$El_{cr1} < El < El_{cr2}$, where the fluid inertia and elasticity are comparable. Within this range, the particle tends to focus on an equilibrium position between the duct centreline and the wall. The third critical elastic number 
$El_{cr3}$ governs the migration efficiency of oblate particles in highly elastic flows. If 
$El > El_{cr3}$, then the efficiency of particle migration to the duct centreline is reduced. Consequently, it is deduced that the oblate particle has the highest migration efficiency at 
$El = El_{cr3}$.
Sketch of the three critical elastic numbers in the present flow system.
Similar to the elasto-inertial migration of spherical particles, the first two critical elastic numbers (
$El_{cr1}$ and 
$El_{cr2}$) can be regarded as the consequence of the competition between fluid inertia and fluid elasticity. However, the mechanism responsible for 
$El_{cr3}$ has not been reported yet. The present results in figure 2(b) evidence that the generation of 
$El_{cr3}$ is directly associated with the particle migration velocity. By increasing 
$El$, the oscillating amplitude of particle migration velocity is also notably enlarged. The sign of the particle wall-normal velocity changes alternately during the particle migration in the highly elastic flow at 
$El = 0.1$. Such fluctuations of particle migration velocity result in the particle migrating to its final equilibrium position in an oscillating way, which delays the particle arrival at the duct centreline.
Based on the above analysis, the lateral migration of the oblate spheroid in figure 2 can be classified into two types, i.e. inertial migration and elasto-inertial migration, as shown in table 1. Moreover, table 1 also identifies the dominant factors in each migration type at different 
$El$.
Dominant factors in the lateral migrations of spheroids shown in figure 2.
It is well known that the migration of oblate particles is significantly affected by particle orientation. Tai & Narsimhan (Reference Tai and Narsimhan2022) have demonstrated that the particle migration velocity during its viscoelastic migration is proportional to the projection length of the particle on the shear-gradient (
$x$–
$y$) plane. This fact is also reflected in figure 4, which shows the relationship between the particle orientation and migration velocity. When the particle symmetry axis is approaching the streamwise direction (increasing 
$| {\cos ( {{\theta _z}} )} |$), the wall-normal projection length of the particle, 
$D\,| {\cos ( {{\theta _z}} )} |$, also increases, as shown in figure 5. Consequently, according to the analysis of Tai & Narsimhan (Reference Tai and Narsimhan2022), the particle migration velocity is augmented. However, from figure 4, the particle migration velocity is quickly reduced once the particle nearly aligns its major diameter along the wall-normal direction (a configuration shown in figure 5a).
Time evolution of the particle orientation and migration velocity at (a) 
$El = 0.05$ and (b) 
$El = 0.1$.
Two particle orientations during the particle migration in viscoelastic flows, 
$El = 0.1$. (a) Particle symmetry axis is nearly perpendicular to the 
$y''$-axis. 
$(\textit {b})$ Particle symmetry axis is nearly parallel to the streamwise 
$y''$-axis. The particle orientation is represented in the co-moving frame (
$x''$–
$y''$–
$z''$) with its origin in the particle centre. The coordinate axes of the co-moving frame are parallel to the corresponding axes of the inertial frame.
The above reduction of particle migration velocity observed when the particle symmetry axis approaches the streamwise direction is attributed to the shear-gradient lift force exerted on the particle. During the particle lateral migration in elasto-inertial flows, the particle mainly experiences three kinds of hydrodynamic forces, including the shear-gradient lift force, the wall-normal drag force and the elastic lift force. The magnitude of the shear-gradient lift force on an oblate particle is formulated as (Nizkaya et al. Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020)
\begin{equation} {F_{sl}} = \rho {a^3}bG_m^2{C_l}, \end{equation}
where 
$F_{sl}$ is the magnitude of the shear-gradient lift force, 
$G_m$ denotes the mean flow shear rate on the particle surface, and 
$C_l$ is the lift coefficient. Equation (3.1) suggests that the shear-gradient lift force is proportional to the flow shear rate on the particle surface 
$G_m$; when the particle symmetry axis approaches the streamwise direction, the projection length of the particle spanned in the shear plane is enlarged, thus the mean flow shear rate 
$G_m$ increases. This causes an increased shear-gradient lift force, which acts as a resisting force during the particle migration, finally resulting in the reduction of the particle migration velocity.
From the phenomenal point of view, the above features of particle migration velocity induce the ‘oscillating’ behaviour in the elasto-inertial migration of oblate particles, as shown in figure 6. It is found that the oscillation in particle position appears on the transition point of particle orientation, consistent with the particle migration velocity shown in figure 4. Compared to the spherical particle migration, the oscillating migration is the most prominent feature for oblate particles in the elasto-inertial flows.
Oscillating behaviour in the elasto-inertial migration of oblate particles, for (a) 
$El = 0.05$, (b) 
$El = 0.1$.
Finally, to clarify the effect of 
$X_{p0}$ on the oscillating migration of the oblate spheroid, figure 7 shows the spatial trajectory of the spheroids released from different initial positions, i.e. 
$X_{p0} = 0$, 
$-0.125$, 
$-0.25$ and 
$Y_{p0} = 0.25$. In this simulation, 
$\Delta t = 2\times 10^{-3}$. In figure 7, it is found that the initial 
$X_{p0}$ does not alter the oscillating characteristics in the lateral migration of the oblate spheroid. In contrast, the focusing pattern of the spheroid is determined by both the particle initial position and fluid elasticity. For example, in the viscoelastic flow at 
$El = 0.01$, the spheroids with different 
$X_{p0}$ are attracted to different equilibrium positions (figure 7b), including the diagonal equilibrium position 
$(X_{peq}, Y_{peq})=(-0.1638, 0.1635)$ and the centre equilibrium position 
$(X_{peq}, Y_{peq})=(0, 0.1180)$. This focusing pattern for the oblate spheroid with multiple equilibrium positions in the present study is similar to that of spherical particles (Yu et al. Reference Yu, Wang, Lin and Hu2019). This phenomenon is attributed to the non-uniform distribution of the polymer deformation (Del Giudice et al. Reference Del Giudice, D'Avino, Greco, Netti and Maffettone2015) and flow shear rate in the duct.
Lateral migration of the oblate spheroid with different 
$X_{p0}$ at 
$Re = 15$. (a,c) Time evolution of particle position 
$Y_{p0}$. (b, d) Particle trajectory projected on the cross-plane of the duct. Here, (a,b) 
$El = 0.01$, (c,d) 
$El = 0.1$. The dots of different colours in (b,d) denote the different initial positions of the oblate spheroid.
Moreover, figure 7(b) shows that the final equilibrium position of spheroids near the duct side-wall (
$X_{p0} = -0.125$ and 
$-0.25$) is higher than that of the oblate spheroid within the symmetry plane (
$X_{p0} = 0$). This result can be explained by the scaling analysis of the shear-induced and elastic lift forces. By substituting the definition of the first normal stress difference 
${N_1} = \tau _p^{zz} - \tau _p^{yy} \propto | \boldsymbol{\mathsf{B}} |$, and the power-law function of 
${N_1}$ in the dilute polymeric solution, where 
${N_1} = AG_m^n$ (
$1 < n \leqslant 2$) (Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007), and 
$A$ is the model parameter, the elastic lift force can be approximated as
\begin{equation} {F_{el}} = \oint_{{\varGamma _p}} {\boldsymbol{\tau}_p} \boldsymbol{\cdot} \boldsymbol{n}\,\mathrm{d}s = \oint_{{\varGamma _p}} {\frac{{{\mu _p}}}{\lambda }\,( {\boldsymbol{\mathsf{B}} - \boldsymbol{\mathsf{I}}} )} \boldsymbol{\cdot} \boldsymbol{n}\,\mathrm{d}s \propto G_m^n\quad ( {1 < n \leqslant 2} ). \end{equation} From the comparison between (3.1) and (3.2), it is concluded that the shear-induced lift force grows more quickly than the elastic lift force with increasing flow shear rate 
$G_m$. Therefore, the shear-induced lift force is enhanced more significantly than the elastic lift force from the symmetry plane to the side-wall of the duct, resulting in a higher equilibrium position for the spheroid released near the duct side-wall, i.e. 
$X_{p0} = -0.125$ and 
$-0.25$.
3.2. Scaling of particle equilibrium position
From the applied perspective, the scaling law of the particle equilibrium position plays an important role in the controlling strategy for non-spherical particles in microfluidic applications. To this end, this subsection discusses the characteristics of the particle equilibrium position in detail.
Figure 8 shows the phase diagram of the particle equilibrium position. As mentioned in § 3.1, there are three critical elastic numbers in the elasto-inertial migration process of oblate particles. We focus primarily on the second critical elastic number, 
$El_{cr2}$, which represents the critical fluid elasticity corresponding to 
$Y_{peq} = 0$. First, from figure 8(a), 
$El_{cr2}$ is affected by 
$Re$ in the flow system: for 
$Re > 15$, 
$El_{cr2}$ is almost unchanged with 
$Re$. However, in the low-
$Re$ flows (
$Re < 15$), 
$El_{cr2}$ is larger than that in the viscoelastic flow with higher 
$Re$. This phenomenon is similar to that observed in the elasto-inertial migration of spherical particles (Li et al. Reference Li, McKinley and Ardekani2015). More clearly, in figure 8(b), by the scaling law 
$Wi \sim 0.012\,Re$ in the 
$Wi$–
$Re$ space, the particle focusing patterns can be divided into two regimes. Such a scaling law, corresponding to 
$El_{cr2} = 0.012$, is close to that for the spherical particles (
$El_{cr2} = 0.01$). This result implies a similarity in the critical fluid elasticities required to achieve the duct-centring focusing for both spherical and spheroidal particles.
Phase diagrams of the particle equilibrium position. (a) Particle equilibrium position in 
$El$–
$Re$ space. (b) Particle equilibrium position in 
$Re$–
$Wi$ space. The symbol colour represents the 
$y$ coordinate of the particle equilibrium position 
$Y_{peq}$. The symbol shape represents the different particle focusing patterns: the square symbol denotes the particle focusing on the position between the duct wall and the centreline, while the circle symbol represents the particle being attracted to the duct centreline.
To elaborate further on the effects of fluid elasticity 
$El$ and fluid inertia 
$Re$ on the particle equilibrium position, figure 9 shows the modulation of the particle equilibrium position by 
$El$ and 
$Re$, respectively. Regarding the fluid elasticity 
$El$, figure 9(a) suggests that the particle equilibrium position is a monotonic function of 
$El$. With increasing 
$El$, the fluid elasticity is enhanced and promotes the particle equilibrium position, gradually covering to 
$Y_{peq} = 0$.
Effects of fluid elasticity and fluid inertia on the particle equilibrium position. (a) Effect of the fluid elasticity 
$El$. (b) Effect of the fluid inertia on the particle equilibrium position in elasto-inertial flows at 
$El =0.01$ and 
$0.05$. (c) Particle equilibrium position in Newtonian flows (
$El = 0$) at different 
$Re$. In (b), ED and FD denote the elasticity domination and fluid inertia domination regions, respectively.
In contrast to fluid elasticity, the effect of fluid inertia on the particle equilibrium position becomes more compliant. From figure 9(b), it can be observed that when keeping 
$El = 0.01$ in the flow system, the equilibrium position of the oblate particle varies non-monotonically with 
$Re$. Such a non-monotonic response of the particle equilibrium position to 
$Re$ is also evidenced by experiments (Holzner et al. Reference Holzner, Stavrakis and de Mello2017). More specifically, the variation of the particle equilibrium position to 
$Re$ in figure 9(b) can be approximately classified into two regions, including the fluid elasticity domination region (
$Re < 15$) and the fluid inertia domination region (
$Re > 15$).
First, within the fluid elasticity domination region (
$Re < 15$), different from the elastic migration of particles, the particle equilibrium position increases when decreasing 
$Re$ in the elasto-inertial flows at a constant 
$El = 0.01$. In such a flow condition, the equilibrium position of the spheroid is determined by the competition between fluid inertia and fluid elasticity. With 
$Re$ decreasing, 
$Wi = El\, Re$ will also be attenuated, causing the effect of fluid elasticity to gradually diminish in the flow system. Thus the elastic lift force on the spheroid is reduced, and the curve of the spheroid equilibrium position increases. The above discussion of the spheroid equilibrium position about 
$Re$ effect in the constant 
$El$ flow can also be found for the spherical particles in the elasto-inertial duct flow (Li et al. Reference Li, McKinley and Ardekani2015).
Moreover, although the fluid inertial effect is also suppressed as 
$Re$ reduces (the weaker fluid inertia helps the particle to focus towards the duct centreline), figure 9(b) suggests that the fluid elasticity is dominant in the flow system at 
$Re < 15$ and 
$El = 0.01$. Thus the equilibrium position of the spheroid is governed mainly by the variation of the fluid elasticity in the small-
$Re$ flow at 
$El = 0.01$.
More physically, the above peculiar phenomenon is associated with the polymer deformation near the particle. According to (2.5), the polymer conformation tensor 
$\boldsymbol{\mathsf{B}}$ is evolved by two effects: a polymer stretching-coiling term 
$\boldsymbol{\mathsf{B}} \boldsymbol {\cdot } ( {\boldsymbol {\nabla } \boldsymbol {u}} ) + {( {\boldsymbol {\nabla } \boldsymbol {u}} )^{\rm T}} \boldsymbol {\cdot } \boldsymbol {B}$, and a relaxation term 
$1/\lambda (\boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{B}})$. The first term represents the polymer deformation induced by the flow shear rate, while the second term denotes the relaxation effect of the polymer. Equation (2.5) indicates that the larger the flow shear rate (
${\boldsymbol {\nabla } \boldsymbol {u}}$), the stronger the polymeric molecule is stretched near the particle surface, as shown in figures 10(a,b). On the other hand, in the flow with a constant 
$El$, the relationship 
$Wi = El\, Re$ implies that the larger 
$Re$ corresponds to the larger 
$Wi$ or 
$\lambda$ in the flow system. Both of the above two aspects will generate the larger polymeric stress (figures 10c,d) near the particle in the flow at higher 
$Re$. Thus in the first region (
$Re <15$), the reducing particle equilibrium position occurs when increasing 
$Re$. In the second region (
$Re > 15$), the dominant effect in the flow system changes to the fluid inertia; the enhanced shear-gradient force induced by the larger 
$Re$ pushes the particle closer to the duct wall, corresponding to the larger particle equilibrium position.
Contours of the flow field near the particle: (a) shear rate at 
$Re = 5$; (b) shear rate at 
$Re = 15$; (c) conformation tensor 
$B_{zz}$ at 
$Re = 5$; (d) conformation tensor 
$B_{zz}$ at 
$Re = 15$. In (a,b), SR denotes the shear rate, which is scaled by 
$W_p/D$.
When 
$El$ increases further to 
$El = 0.05$, the fluid elasticity is fully dominant in the flow system, i.e. 
$El > El_{cr2} = 0.012$; as shown in figure 8, the spheroid will be attracted to the duct centreline at different 
$Re$. In Newtonian flows at 
$El = 0$, with 
$Re$ increasing, the enlarged fluid inertial effect causes the spheroid to get closer to the duct wall, thus the equilibrium position of spheroids is increased by 
$Re$, as shown in figure 9(c).
3.3. Particle rotational dynamics during the cross-flow migration
In addition to the oscillating migration of the oblate particle, another prominent feature of the non-spherical particles in elasto-inertial flow is the particle rotational dynamics, which is essential for assembling non-spherical bio-particles to realize a desired orientation with microfluidic technology. Thus this subsection analyses the rotation and orientation of the oblate particle during its elasto-inertial migration.
3.3.1. Particle rotation
First, figure 11 shows the three-dimensional (3-D) trajectories of particle translation and orientation. From figure 11, regardless of inertial (
$El=0$) or the elasto-inertial (
$El \ne 0$) migration, the oblate particle exhibits the oscillating migration. Meanwhile, the particle shows the kayaking-like rotation mode, revealing the coupling nature between the migration and rotation of non-spherical particles. Furthermore, figures 11(a–e) also demonstrate that the final particle orientation modes are determined by the coupling effect of fluid inertia and fluid elasticity. Note that although the final orientation of the oblate particle during inertial (
$El = 0$) and elasto-inertial (
$El= 0.1$) migration is the same, the specific particle rotation dynamics are different.
The 3-D trajectories of particle translation and rotation during its elasto-inertial migration, 
$Re = 15$: (a) 
$El = 0$, (b) 
$El = 0.01$, (c) 
$El = 0.03$, (d) 
$El = 0.05$, (e) 
$El = 0.1$.
For the particle rotation, figure 12 shows the particle rotation rate around its symmetry axis in the viscoelastic flow at different fluid elasticities. Overall, in the inertial migration (
$El = 0$), when the particle approaches its equilibrium position near the wall, the flow shear rate on the particle increases, thus the particle rotation rate is enlarged. However, the scenario becomes different in viscoelastic flows. By increasing fluid elasticity (
$El$), the average particle rotation rate is attenuated, causing a larger particle rotation period. This phenomenon can be explained from two aspects: (1) the tension along the streamline generated by fluid elasticity hinders the particle rotation rate; (2) as the particle moves towards the duct centreline, the flow shear rate acting on the particle diminishes, further contributing to the reduction of the particle rotation rate. Especially in the highly elastic flows at 
$El = 0.05$ and 
$0.1$, the particle rotation rate is reduced to zero, indicating a rotationless state for particles. This rotationless state for the oblate particle at the duct centreline is consistent with the previous experimental results on the elasto-inertial migration of cells (Holzner et al. Reference Holzner, Stavrakis and de Mello2017).
Rotation rate around the particle symmetry axis during its elasto-inertial migration with 
$Re = 15$.
Besides the above particle rotation rate, another important feature to characterize the particle rotational dynamics is the orbit drift of the particle symmetry axis. Figure 13 shows the temporal evolution of orbit parameter 
$C_b$, which is defined as (Yu, Phan-Thien & Tanner Reference Yu, Phan-Thien and Tanner2007)
\begin{gather} {C_b} = \frac{C}{{C + 1}}, \end{gather}
\begin{gather} {C = AR^{ - 1}\tan \theta \sqrt {(AR^2\sin^2\phi + \cos^2\phi)}} , \end{gather}
where 
$\theta$ and 
$\phi$ are the polar angle and azimuthal angle of the particle symmetry axis, respectively.
Orbit parameter of the oblate spheroid during its elasto-inertial migration at 
$Re = 15$.
Overall, similar to that in viscoelastic linear shear flows, 
$C_b$ is a nonlinear function of time. Compared to the inertial migration in Newtonian flow (
$El = 0$), the drifting rate of the particle orbit is suppressed by the effect of fluid elasticity. Moreover, figure 13 shows that the effect of fluid elasticity on the drift of orbit is different during particle migration. In the earlier stage of particle migration (
$t < 200$), the stronger the fluid elasticity, the larger the amplitude of 
$C_b$.
The convergent value of the particle orbit parameter (
$C_{beq}$) reveals the equilibrium rotation state of oblate particles. In the later stage of particle migration (
$t > 200$), 
$C_{beq}$ is almost unchanged with time in Newtonian (
$El = 0$) and highly elastic (
$El > 0.05$) flows. This behaviour of 
$C_{beq}$ suggests that the particle can achieve a stable rotation state or orientation at the particle equilibrium position. More importantly, the results in figure 13 inspire us to say that the orientation of complex-shaped particles can be controlled by tuning the rheology of suspension fluid in microfluidic applications.
Different from highly elastic flow (
$El = 0.05$ and 
$0.1$), the particle rotation becomes more complicated in the viscoelastic flow at 
$El = 0.03$. In figure 13, according to the trend of orbit parameter 
$C_b$, it is obvious that the evolution of 
$C_b$ at 
$El = 0.03$ can be divided into two stages. To conveniently explain the above physical picture in figure 13, figure 14 shows the 3-D rotational dynamics of the oblate particle at 
$El = 0.03$, including the 3-D trajectory of the particle tip, particle rotation rate and orbit parameter 
$C_b$. From figure 14(c), in the first stage (
$t < 240$), the particle orbit parameter 
$C_b$ gradually decreases, representing that the particle symmetry axis approaches the vorticity direction spirally (figure 14a). When 
$t > 240$, the particle rotation enters into the second stage, where 
$C_b$ begins to increase with time (figure 14c), indicating that the particle symmetry axis crosses the shear plane (figure 14a). In figure 14(c), the comparison of particle position 
$Y_p$ and orbit parameter 
$C_b$ clearly shows that the second stage starts when particle just reaches the duct centreline. Note that the second stage of particle rotation is a very slow process, in which the angular velocity is nearly zero, as shown in figure 14(b).
Rotational dynamics of the oblate particle in the viscoelastic flow at 
$El = 0.03$ and 
$Re = 15$. (a) A 3-D trajectory of the particle symmetry axis in the particle co-moving frame; the black and red orbits represent the particle orbits in two rotation stages, respectively. (b) Particle angular velocities in the particle frame. (c) Orbit parameter; the green dot denotes the particle orientation when 
$Y_p = 0$ at 
$t =240$.
In the duct centre region, the flow shear rate on the particle is weak, thus the particle rotation is mainly affected by the fluid inertial torque and fluid elastic torque induced by the particle slip velocity. To elaborate on the mechanism of the peculiar particle rotation mode presented in figure 14, we compare the streamwise velocity of the fluid surrounding and far away from the particle as shown in figure 15. Overall, figure 15 indicates that, similar to the spherical particle in elasto-inertial flows (Li et al. Reference Li, McKinley and Ardekani2015), the particle always slightly lags the fluid at the same lateral position of the particle (
$Y_{peq}$) due to the weak particle inertia. Moreover, the presence of the oblate particle just modulates the fluid velocity in a narrow region whose width is approximately one particle major diameter 
$D$. Beyond this region, the fluid velocity quickly recovers to that of the fluid in the far field.
Streamwise velocity of the fluid surrounding the particle and far undisturbed flow field with 
$Re = 15$: (a) 
$El = 0$, (b) 
$El = 0.01$, (c) 
$El = 0.05$. The filled coloured dots and the dashed coloured lines represent the equilibrium position of the particle centre. The coloured shaded area denotes the region with one particle major diameter 
$D$.
To quantify the relative importance of fluid inertial and fluid elastic effects, figure 15(c) approximately measures the slip velocity between the particle centre and the far-field fluid as 
$u_s = 0.0932$. Then the particle Reynolds number and Weissenberg number can be determined by 
$Re_p = u_sD/\nu _0 = 0.0874$ and 
$Wi_p = \lambda u_s/D= 0.0419$. The corresponding elastic number is 
$El_p = Wi_p/Re_p \sim 0.5$. Dabade, Marath & Subramanian (Reference Dabade, Marath and Subramanian2015) predicted theoretically the critical elastic number as 
$El_{pcr} = 1.337$, at which the fluid elastic torque fully governs the oblate particle (
$AR = 0.5$) orientation. From the present relationship 
$El_p < El_{pcr}$, it can be deduced that the fluid elasticity is not strong enough to be fully dominant in the flow system. Consequently, the fluid inertial and elastic effects appear to be comparable in the elasto-inertial duct flow at 
$El = 0.03$. From the previous pioneer studies, it is well known that fluid inertial torque and fluid elastic torque induce opposite particle orientation modes. Fluid inertial torque can induce the particle orientation maximizing the drag force (the so-called broadside-on mode Dabade et al. Reference Dabade, Marath and Subramanian2015). Conversely, the fluid elastic torque drives the particle to orientate with the minimum drag force, i.e. longside-flow alignment mode, which is termed as the longside-on mode in Dabade et al. (Reference Dabade, Marath and Subramanian2015). Thus the competition between the above two comparable hydrodynamic torques finally induces the unstable orientation: particle orientation gradually transitions between broadside-on mode and longside-flow alignment mode. This unstable particle orientation further generates the slight fluctuations in the particle equilibrium position shown in figure 11(c).
3.3.2. Steady particle orientation
The previous subsubsection mainly focuses on the 3-D rotational dynamics of the oblate particle during its migration. Figure 11 shows that the particle can also exhibit multi-steady orientation modes in the elasto-inertial flows. As mentioned before, these particle orientation modes are crucial to particle assembling in microfluidic applications. This motivates the following study on the steady orientation modes of the oblate particle in elasto-inertial duct flows with 
$El = 0$, 
$0.01$, 
$0.05$ and 
$0.1$.
Figure 16 demonstrates the steady orientation mode of the oblate particle during its elasto-inertial migration. With increasing the fluid elasticity from 
$El = 0$ to 
$0.1$, there are three kinds of steady orientation modes for the oblate particles.
(i) Log-rolling mode. The oblate spheroid is spinning around the vorticity direction in the shear plane, as shown in figure 16(a).
(ii) Kayaking-like mode. The oblate spheroid performs both precession and nutation around the vorticity direction (figure 16b), causing the spheroid to rotate in a kayaking orbit. The kayaking resembles the motion of a kayak paddle (Rosén et al. Reference Rosén, Do-Quang, Aidun and Lundell2015).
(iii) Longside-flow alignment mode. The major axis of the spheroid is along the streamwise direction, as shown in figures 16(c,d). According to the alignment of the spheroid symmetry axis in the channel cross-section, the longside-flow alignment mode further includes longside-flow alignment mode A and longside-flow alignment mode B.
Steady orientation mode of the oblate particle (
$AR = 0.5$) in elasto-inertial duct flows with 
$Re = 15$: (a) log-rolling mode, 
$El = 0$; (b) kayaking-like mode, 
$El = 0.01$; (c) longside-flow alignment mode A, 
$El = 0.05$; (d) longside-flow alignment mode B, 
$El = 0.1$. (a i,b i,c i,d i) The 3-D trajectories of the particle symmetry axis in the particle co-moving frame. The 3-D trajectory of the particle symmetry axis is also projected on the 
$y''\unicode{x2013}z''$ plane. The red and green dots represent the particle's initial and finally steady orientations, respectively. (a ii,b ii,c ii,d ii) The particle rotation rates (
${\omega '_{p = x,y,z}}$) within the particle frame. The corresponding flow velocity is sketched by the black arrows in each plot.
The mechanisms governing the above spheroid orientation modes are described as follows.
First, in the Newtonian duct flow (
$El = 0$), the oblate particle finally orientates in the log-rolling mode at its equilibrium position under the fluid shear effect, as shown in figure 16(a). According to the stability analysis (Einarsson et al. Reference Einarsson, Candelier, Lundell, Angilella and Mehlig2015), the log-rolling mode has been verified as a stable orientation mode for the oblate particle in the Newtonian shear flow.
Second, in the weakly elastic flow (
$El = 0.01$), when the fluid inertia and fluid elasticity are comparable, a kayaking-like mode is observed for the oblate particle. Similarly, our previous study has demonstrated a semblable orientation mode, i.e. asymmetry kayaking mode, for the prolate spheroid in elasto-inertial shear flows (Li, Xu & Zhao Reference Li, Xu and Zhao2023). The mechanisms behind this orientation mode for spheroids in elasto-inertial flows are similar. The fluid convective inertia prompts the spheroid to align its symmetry axis to the vorticity direction (Einarsson et al. Reference Einarsson, Candelier, Lundell, Angilella and Mehlig2015), whereas the fluid elastic effect strives to orient the spheroid symmetry axis perpendicular to the vorticity direction (Bartram, Goldsmith & Mason Reference Bartram, Goldsmith and Mason1975; D'Avino & Maffettone Reference D'Avino and Maffettone2015). From the above opposite particle orientation modes, it is clear that the present kayaking-like mode for the oblate particle during its elasto-inertial migration originates from the competition between fluid inertial and fluid elastic effects. Furthermore, comparing figures 16(a–c), it is found that fluid elasticity leads to the orbit shape of the particle symmetry axis being thinner, which is also evidenced by the orbit parameter 
$C_{b}$ shown in figure 13.
Finally, for the longside-on mode observed in the highly elastic flow (
$El = 0.1$), the particle equilibrium position becomes the duct centreline, where the particle elastic number is estimated as 
$El_p = Wi_p/Re_p \sim 1.598$. Compared to the critical elastic number 
$El_{pcr} = 1.337$ predicted by Dabade et al. (Reference Dabade, Marath and Subramanian2015), the fluid elastic torque is fully dominant (
$El_p > El_{pcr}$) in the flow system at 
$El = 0.1$. Thus the present longside-on mode for oblate spheroids is induced by the fluid elastic torque. It is worthwhile to note that although the configurations of particle steady orientation modes in Newtonian (
$El = 0$) and highly elastic (
$El = 0.1$) flow are similar, as shown in figures 15(a,d), the particle rotational behaviours in these two orientation modes are different. This can be distinguished by the particle rotation rates in figure 15. It is found that the particle rotation rate finally converged to 
$\omega _z = -1.021$ in the Newtonian flow, while 
$\omega _z = 0$ in the viscoelastic flow. This difference in particle rotation states evidences that the particle rotation process is greatly suppressed by the fluid elasticity.
Moreover, from the particle orientation modes in figures 16(c,d), the specific configuration of the particle orientation mode in a highly elastic flow is affected by the fluid elasticity. Although the oblate particle aligned its major diameter along the streamwise direction, i.e. longside-flow alignment mode, in figures 16(c,d), the angles between the particle symmetry axis and vorticity direction are different. In the flow at 
$El = 0.05$, the particle aligned its symmetry axis along the diagonal of the duct cross-section. However, with an increase in fluid elasticity (
$El = 0.1$), the particle symmetry axis is approximately attracted to the vorticity direction. These different configurations within the same orientation mode can be attributed to the change in the microstructure of polymer deformation in the flow field. The trace of the conformation tensor represents the polymer deformation, as shown in figure 17. When fluid elasticity is 
$El = 0.05$, figure 17(a) shows that the polymer is stretched along the diagonal of the duct cross-section. With further increasing 
$El$ to 
$El = 0.1$, the deformation of the polymer is enhanced and is nearly uniform in the space with the same distance to the particle. In other words, the polymer is stretched without a preferential direction around the particle at 
$El = 0.1$, resulting in the configuration of longside-flow alignment mode presented in figure 16(b).
Microstructure of polymer deformation around the oblate particle with 
$Re = 15$. (a) Longside-flow alignment mode A, 
$El =0.05$. (b) Longside-flow alignment mode B, 
$El = 0.1$. The grey contour surface represents 
${\rm trace}(\boldsymbol {B}) =10$.
4. Conclusions
The cross-flow migration of particles in viscoelastic flows is widely used in microfluidic applications. The focusing and orientational behaviours of spheroids in elasto-inertial flow can be used potentially for high-throughput particle separation and orientation control in microfluidic devices. To explore the physical picture in the elasto-inertial migration and orientation of spheroids, in this study, we numerically investigated the elasto-inertial migration of a single neutrally oblate particle (
$AR = 0.5$) in the viscoelastic duct flow with the immersed boundary method. The viscoelastic rheology is described by the Oldroyd-B model. The characteristics of the elasto-inertial migration of the oblate particle are first analysed. Then the mechanisms governing the particle rotational dynamics and orientation modes during the particle migration are also elaborated on in this study.
The main conclusions for the elasto-inertial migration of the oblate particle are as follows.
(i) Different from spherical particles, the lift force exerted on the oblate particle is altered by the particle orientation; the oblate particle exhibits ‘oscillating’ migration behaviour in elasto-inertial flows. The particle migration velocity is a nonlinear function of time. The transition of the particle migration velocity occurs when the particle symmetry axis approaches the streamwise direction.
(ii) The particle equilibrium position is determined by the elastic number (
$El$): the larger the fluid elasticity, the closer the particle equilibrium position is to the duct centreline. The duration for the oblate particle to approach the duct centreline changes non-monotonically with $El$, suggesting the existence of three critical elastic numbers in the present flow system. Within the present parameter range, the second critical elastic number $El_{cr2}$, which determines the critical fluid elasticity corresponding to $Y_{peq}= 0$, satisfies a scaling law $Wi\sim 0.012\,Re$.
The main conclusions for the rotational dynamics and steady orientation modes of the oblate particle are as follows.
(iii) The particle orientation modes during its elasto-inertial migration include stable and unstable orientation modes. Due to the competition between the fluid inertial and fluid elastic effects, the oblate particle shows an unstable rotation state at the duct centreline. Additionally, the particle rotation rate and the drifting rate of the particle orbit are both greatly hindered by the fluid elasticity.
(iv) We find three kinds of steady orientation modes for oblate particles: log-rolling mode, kayaking-like mode and longside-flow alignment mode. The log-rolling mode of the oblate particle in Newtonian flow (
$El = 0$) is induced by the fluid convective inertia. The kayaking-like mode is caused by the competition between fluid inertial and fluid elastic effects. Finally, the longside-flow alignment mode in the highly elastic flow originates from the fluid elastic torque. More interestingly, the specific configuration of the longside-flow alignment mode is attributed to the microstructure of polymer deformation in the flow field.
In summary, the main contributions of the present work are: (1) a first attempt to elaborate numerically on the elasto-inertial migration of non-spherical particles in viscoelastic duct flow; and (2) to give the link between the spheroid orientation and its elasto-inertial migration. The results of the present study could shed new insights into the cross-flow migration of complex particles in non-Newtonian fluids. Furthermore, from the applied perspective, the present results could also be used potentially to design the rheology-based controlling strategy for guiding particles to achieve optimal focusing and orientation in microfluidic applications without the need for external forcing fields.
Finally, several works are awaiting for future study: (1) exploring the effect of other rheological parameters of viscoelastic fluids (viscosity ratio, shear-thinning rheology, elastoviscoplastic fluids) on the migration of spheroids in complex fluid flows; (2) surface-stress-based analysis of oscillating migration of spheroids in the elasto-inertial flows; (3) identifying the exact critical elastic numbers, which are important for designing the rheology-based controlling strategy for manipulating the spheroids in microfluidic applications.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.572.
Funding
The authors are grateful for the support of the Natural Science Foundation of China through grant nos. 92252104, 12388101 and 92252204. Y.L. is also supported by the Postdoctoral Fellowship Program of CPSF under grant no. GZB20230360.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation of the numerical approach
The present study aims to investigate the elasto-inertial migration of oblate particles in viscoelastic duct flows. The particle dynamics is affected by the fluid inertial and elastic effects simultaneously. Thus it is important to check the ability of the present numerical approach to capture the effects of fluid inertia and fluid elasticity on the particle behaviour. To this end, we validated the numerical methods from three aspects: (i) the elastic migration of spheroids, (ii) the inertial migration of spheroids, and (iii) the elasto-inertial migration of the sphere in viscoelastic duct flow.
A.1. Elastic migration of prolate spheroids
For the elastic migration of spheroids, the present study focuses on the prolate spheroids. Although the elastic migration of the oblate spheroid has also been predicted theoretically, the suspension fluid is modelled by second-order fluid, which is generally utilized to model the polymeric fluid with weak fluid elasticity. Considering the significant fluid elasticity in the elasto-inertial channel flow, it is thus necessary to check the capability of the present numerical approach to capture the particle dynamics in highly elastic flows. D'Avino et al. (Reference D'Avino, Hulsen, Greco and Maffettone2019) simulated the lateral migration of a single prolate spheroid in a viscoelastic fluid in the wide-slit microchannel. Here, we first use it to validate the present numerical method.
In this validation case, the particle aspect ratio is set as 
$AR = 2$. Based on the channel height and average flow velocity, the bulk Reynolds number is 
$Re = 0.1$, thus the fluid inertial effect can be neglected in the flow system. The Weissenberg number is 
$Wi = 0.5$. Other simulation parameters are consistent with the reference study (D'Avino et al. Reference D'Avino, Hulsen, Greco and Maffettone2019). The comparison of results is shown in figure 18, in which the calculated particle position is consistent with the reference results from D'Avino et al. (Reference D'Avino, Hulsen, Greco and Maffettone2019). Similar to the spherical particle, the prolate spheroids released from different initial positions can be attracted to the channel centreline. Figure 18 indicates that the present numerical method can readily estimate the elastic migration of the non-spherical spheroids in viscoelastic fluids.
Comparison of calculated particle position with the reference results (D'Avino et al. Reference D'Avino, Hulsen, Greco and Maffettone2019), denoted by dots. The line represents the calculated results.
Moreover, figure 19 further shows the distribution of the first normal stress difference 
$N_1$ in the flow field during the particle elastic migration. The results suggest that 
$N_1$ above the particle region is always larger than that in the region below the particle. Such a gradient of 
$N_1$ in the flow field can generate an elastic lift, which pushes the prolate particle to the channel centreline (Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007; Raoufi et al. Reference Raoufi, Mashhadian, Niazmand, Asadnia, Razmjou and Warkiani2019). Finally, for the orientation of the prolate particle in figure 19, it is found that the prolate particle exhibits the flow alignment mode at the channel centreline, which is consistent with the recent experimental results (Langella et al. Reference Langella, Franzino, Maffettone, Larobina and D'Avino2023). All the above results agree with the simulation from D'Avino et al. (Reference D'Avino, Hulsen, Greco and Maffettone2019).
Distribution of the first normal stress difference 
$N_1$ in the flow field: (a) 
$t=0$, (b) 
$t=0.596$, (c) 
$t=4.1969$, (d) 
$t=7.1947$, (e) 
$t=10.1925$, (f) 
$t=69.5485$.
A.2. Inertial migration of oblate spheroids
The inertial migration of the oblate spheroid in the Newtonian channel flow is used to check the capability of the present method to capture the fluid inertial effect on particle migration. In this test case, the Reynolds number is 
$Re = 22$, representing the finite fluid inertia in the flow system. The particle aspect ratio is 
$AR = 0.5$, and the particle blockage ratio is 
$D/H = 0.3$. For the simulation parameters, the mesh resolution is set as 
$\varDelta = 1/32D$, and the time step is 
$\Delta t = 10^{-3}$. The comparison of results is shown in figure 20. First, from figure 20(a), the evolution of the particle position calculated by the present method agrees with the reference results (Nizkaya et al. Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020), indicating that the present method is capable of determining the inertial migration of oblate spheroids. Note that the small deviation between the two results might be caused by the different numerical methods.
Inertial migration of oblate (
$AR = 0.5$) and spherical (
$AR = 1$) particles: (a) particle equilibrium position; (b) particle migration velocity. The open and filled coloured dots denote the initial and final particle migration velocities, respectively.
Figure 20(a) indicates that the sphere particle is faster than the oblate particle, but the equilibrium positions of the two types of particles are almost the same. This suggests that it is difficult to separate the spherical and oblate particles with the same major diameter by the particle inertial migration. Moreover, figure 20(b) also compares the particle migration velocity between oblate and spherical particles. Different from the spherical particle, the particle migration velocity of the oblate particle shows ‘oscillating’ characteristics, which reduce its migration efficiency.
Finally, figure 21 gives the 3-D rotational and migration behaviour of the oblate particle during its inertial migration. It is clearly seen that the oblate particle exhibits a ‘kayaking-like’ rotation state during the inertial migration, and orientates in a ‘log-rolling’ mode at the its equilibrium position. These rotational dynamics are consistent with the results reported by Nizkaya et al. (Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020).
The 3-D rotation and migration of an oblate particle in the Newtonian channel flow at 
$Re = 22$.
A.3. Elasto-inertial migration of spherical particles
As mentioned before, in the elasto-inertial migration of particles, particle behaviour is determined by the competition between fluid inertia and fluid elasticity. Thus it is necessary to validate the performance of the present method in simulating the coupling effect on the particle dynamics. Considering that the elasto-inertial migration of non-spherical particles is still developing, here we used the elasto-inertial migration of spherical particles to further validate our numerical approach. In this test case, to save the computational cost, the present computation size is 
$L = 8D$ in the streamwise direction, which is much smaller than that (
$L = 16D$) in the reference study (Li et al. Reference Li, McKinley and Ardekani2015). Different from Li et al. (Reference Li, McKinley and Ardekani2015), all the present simulations are conducted in the inertial frame with the periodic condition in the streamwise direction. The mesh resolution is 
$\varDelta =1/32D$, and the time step is 
$\Delta t=10^{-3}$. Other rheological parameters of viscoelastic fluid are consistent with the reference study (Li et al. Reference Li, McKinley and Ardekani2015).
The comparison results are demonstrated in figure 22. The particle position is generally consistent with the reference results, especially for the time of the sphere reaching its equilibrium position. The deviation is caused by the smaller computational size in the present study. Overall, the comparison of results shows that the present numerical method could estimate the dependence of particle migration on the coupling effect (
$El$) of fluid inertia and elasticity in elasto-inertial flows. More importantly, figure 22 shows that at the same 
$El$, the shear-thinning effect of polymeric fluids (represented by 
$\alpha = 0.1$ in the Giesekus model) weakens the fluid elasticity, and increases the fluid inertial effect in the flow system. Thus the particle equilibrium position is closer to the channel wall than that in the Oldroyd-B fluid (
$\alpha = 0$). These results further evidence that it is feasible to control the particle positions in the flow field by tuning the rheology of suspension fluids.
Comparison of the calculated elasto-inertial migration of spheres and reference results (Li et al. Reference Li, McKinley and Ardekani2015). Here, 
$\alpha = 0.1$ represents the Giesekus-type fluid with the shear-thinning rheology.
Appendix B. Effect of the computational domain size and mesh convergence
As mentioned in § 3.2, the computational size of the microchannel is set as 
$W\times H\times L = 4D\times 4D\times 8D$. Considering the periodic boundary in the streamwise direction of the channel flow, this appendix checks the effect of the computational size of 
$L$ on the particle elasto-inertial migration. In this test case, the computational size 
$L$ in the streamwise direction is set as 
$L = 8D$, 
$16D$, 
$20D$. The comparison of results simulated in three different microchannels with different 
$L$ is shown in figure 23. From figure 23, it is found that the effect of 
$L$ does not qualitatively change the present results calculated with 
$L = 2H$. The present channel length is similar to that (
$L = H$) in the previous study on the elasto-inertial migration of spheres conducted by Yu et al. (Reference Yu, Wang, Lin and Hu2019). More importantly, Yu et al. (Reference Yu, Wang, Lin and Hu2019) indicated that the shorter channel length (
$L = H$) is important in simulating the particle elasto-inertial migration, because the shorter channel length can reveal the particle hydrodynamic interactions in practical particle focusing or separation applications.
Effect of computational size in the streamwise direction on the particle migration, with 
$Re =15$, 
$El = 0.1$.
Compared to the particle migration in the fluid inertia or fluid elasticity dominant flows, the elasto-inertial migration of the particle is more complex. Considering that the fluid inertia and elasticity are comparable for the spheroid migration in the flow at 
$El = 0.01$ and 
$Re =15$, we focused mainly on the effect of mesh resolution on the case of spheroid migration in the elasto-inertial flow at 
$El = 0.01$ and 
$Re =15$.
In the present test cases, we simulated the migration of the oblate spheroid in the elasto-inertial duct flow on four types of meshes with different spatial resolution, and the results are shown in figure 24. From figure 24(a), the time evolution of the particle mass centre calculated on two finer meshes (
$\varDelta = 1/32D$ and 
$1/64D$) agree well with each other. Furthermore, figures 24(b,c) show the particle lateral velocities at different positions and orientations during the particle migration. It is found that the effect of the mesh resolution on the particle migration is negligible when 
$\varDelta = 1/32D$. These results suggest that the present mesh resolution (
$\varDelta = 1/32D$) is reliable in capturing the oscillating characteristics of the oblate spheroid in the elasto-inertial migration.
Effect of mesh resolution on the elasto-inertial migration of the oblate spheroid at 
$El = 0.1$ and 
$Re =15$. (a) Time evolution of particle position 
$Y_p$. (b) Particle migration velocity at different position 
$Y_p$. (c) Particle migration velocity at different orientation.
Open access
