Oscillation migration of a neutrally buoyant oblate spheroid in the elasto-inertial duct flow

Abstract

The peculiar migration and rotational dynamics of non-spherical particles in non-Newtonian flows stem from the interplay between fluid rheology and fluid inertial effects. In this paper, the cross-flow migration of a neutrally buoyant oblate spheroid (aspect ratio $AR = 0.5$) immersed in the elasto-inertial duct flow is investigated by particle-resolved simulations with the immersed boundary method. Different from spherical particles, due to the orientation-dependent lift force, the oblate spheroid migrates in an oscillating manner in the duct. The travelling period for particles reaching the duct centreline undergoes a non-monotonic change with elastic number, revealing the existence of a critical elastic number governing the migrating efficiency of oblate particles within the present flow system. For the particle rotation and orientation, the present results indicate that the particle rotation rate and orbit drifting rate are both hindered by the fluid elasticity. With increasing the fluid elasticity, three different orientation modes – log-rolling mode, kayaking-like mode and longside-flow alignment mode – are observed successively during the elasto-inertial migration of the oblate spheroid. Potentially, the present results could be used 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.

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. 2023). 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 2017). 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 1961; Martel & Toner 2014) or shear-thinning (Chrit, Bowie & Alexeev 2020; Jannesari Ghomsheh, Jafari & Funfschilling 2023) 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 1966; Leshansky et al. 2007; Naillon et al. 2019). 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 1966; Leshansky et al. 2007; Villone et al. 2011b; Del Giudice et al. 2015; Del Giudice 2019; D'Avino 2021) or elasto-inertial migration with highly fluid elasticity (Lim et al. 2014; Seo et al. 2014; Lu & Xuan 2015). More interestingly, Chaparian et al. (2020) 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. 2019), which motivates the previous abundant research on particle viscoelastic migration in rectangular channel flow. Leshansky et al. (2007) 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. 2017). Later, this two-dimensional focusing pattern has been verified by simulations (Villone et al. 2011a). Furthermore, Seo et al. (2014) and D'Avino et al. (2012, 2017) 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. 2015). 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. 2020), semi-ellipse (Tang et al. 2019) and triangle (D'Avino 2021). 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. 2013) 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. 2013). Del Giudice et al. (2015) 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 2015). In theoretical studies, Wang, Tai & Narsimhan (2020) 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 (2015) and Lu et al. (2015) 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 2015). 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. (2019) 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. 2023). More recently, from the theoretical perspective, Tai & Narsimhan (2022) 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 2020a,b). 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. 2017). Increasing the fluid inertia in the flow system can achieve high-throughput particle focusing in viscoelastic fluid flows (Lim et al. 2014), 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. (2014) 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 (2017) explored the elasto-inertial focusing of bio-particles within the straight microfluidic channel flow of viscoelastic PEO solutions with low viscosities. Li et al. (2015) simulated the elasto-inertial migration of a single spherical particle in the viscoelastic channel flow with finite fluid inertia. Raffiee, Ardekani & Dabiri (2019) 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. 2017; Raoufi et al. 2019). 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. 2022; Jiang, Liu & Tang 2022), 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. 2022). 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

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(2.2)\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:

(2.3)\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:

(2.4)\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

(2.5)\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:

(2.6)\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}

(2.7)\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

(2.8)\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}

(2.9)\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

(2.10)\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}

(2.11)\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 2011)

(2.12)\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

(2.13)\begin{equation} {\boldsymbol{U}_{IB}} = \int_{{\varOmega _p}} \boldsymbol{u}\,\delta ( {\boldsymbol{X} - \boldsymbol{x}} ) \,\mathrm{d}v \end{equation}

and

(2.14)\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

(2.15)\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 2000). By rearranging the discretized governing equations into a monolithic matrix system

(2.16)\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

(2.17) \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. (2022).

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 1980), 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. 2013). Similar to the previous work (Yu et al. 2019), 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. 2019). 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.

Figure 1. 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 2023) and elastic migration in viscoelastic fluid (Holzner et al. 2017; D'Avino et al. 2019; Naillon et al. 2019; Raffiee et al. 2019; Raoufi et al. 2019; Zhou & Papautsky 2020) 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,

(2.18a–d)\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,

(2.19a–d)\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

(2.20)\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

(2.21)\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

(2.22)\begin{equation} El = Wi /Re. \end{equation}

For the particle phase, the particle shape is described by the particle aspect ratio

(2.23)\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. 2020; Tai & Narsimhan 2022). 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. 2020).

The inertia of a particle is characterized by the particle size and density ratio:

(2.24)\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. 2022).

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. 2015; Yu et al. 2019).

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. 2015). 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$.

Figure 2. 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}$.

Figure 3. 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$.

Table 1. 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 (2022) 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 (2022), 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).

Figure 4. Time evolution of the particle orientation and migration velocity at (a) $El = 0.05$ and (b) $El = 0.1$.

Figure 5. 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. 2020)

(3.1)\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.

Figure 6. 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. 2019). This phenomenon is attributed to the non-uniform distribution of the polymer deformation (Del Giudice et al. 2015) and flow shear rate in the duct.

Figure 7. 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. 2007), and $A$ is the model parameter, the elastic lift force can be approximated as

(3.2)\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. 2015). 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.

Figure 8. 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$.

Figure 9. 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. 2017). 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. 2015).

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.

Figure 10. 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.

Figure 11. 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. 2017).

Figure 12. 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 2007)

(3.3)\begin{gather} {C_b} = \frac{C}{{C + 1}}, \end{gather}

(3.4)\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.

Figure 13. 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).

Figure 14. 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. 2015), 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.

Figure 15. 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 (2015) 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. 2015). 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. (2015). 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.

1. (i) Log-rolling mode. The oblate spheroid is spinning around the vorticity direction in the shear plane, as shown in figure 16(a).

2. (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. 2015).

3. (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.

Figure 16. 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. 2015), 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 2023). 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. 2015), whereas the fluid elastic effect strives to orient the spheroid symmetry axis perpendicular to the vorticity direction (Bartram, Goldsmith & Mason 1975; D'Avino & Maffettone 2015). 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. (2015), 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).

Figure 17. 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.

1. (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.

2. (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.

1. (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.

2. (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.