1. Introduction
Elastoviscoplastic (EVP) fluids are a category of non-Newtonian materials that exhibit a dual behaviour: they act like viscoelastic solids when subjected to shear stresses below a critical yield stress, but transform into a fluid-like state once this stress threshold is exceeded (Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014; Saramito Reference Saramito2016; Ewoldt & McKinley Reference Ewoldt and McKinley2017). These fluids are frequently combined with particles and are utilised in various industrial and scientific contexts, including hydraulic fracturing in the oil industry (Barbati et al. Reference Barbati, Desroches, Robisson and McKinley2016), the pumping of concrete in the construction industry (Fataei Reference Fataei2022), transportation of slurries (Li et al. Reference Li, He, Liu, Liu, Huang and Jiang2020), 3-D printing (Eom et al. Reference Eom2025) and the development of biomaterials for drug delivery (Gratton et al. Reference Gratton, Ropp, Pohlhaus, Luft, Madden, Napier and DeSimone2008; Mitragotri & Lahann Reference Mitragotri and Lahann2009). In many applications, the suspended particles are non-spherical, adding significant complexity to the suspension dynamics. Even for simple spheroidal shapes, particle orientation introduces an additional degree of freedom, leading to intricate particle–fluid and particle–particle interactions (Gunes et al. Reference Gunes, Scirocco, Mewis and Vermant2008; Dabade, Marath & Subramanian Reference Dabade, Marath and Subramanian2015; Rosén Reference Rosén2016; Voth & Soldati Reference Voth and Soldati2017). Oblate (flattened, disk-shaped) and prolate (elongated, rod-shaped) spheroids are widely used as canonical representations of non-spherical particles, relevant to a broad range of biomedical and industrial applications, from red blood cells and platelets to synthetic fibres and engineered micro-/nano-particles (Amini et al. Reference Amini, Lee and Di Carlo2014).
A major breakthrough in the dynamics of non-spherical particles was achieved through Jeffery’s pioneering study (Jeffery Reference Jeffery1922) on the angular dynamics of ellipsoidal particles in unbounded shear flow at low Reynolds numbers. This work led to the identification of the Jeffery orbit, where an ellipsoid undergoes a combination of rotation and tumbling, following closed trajectories determined by its aspect ratio and initial orientation. Subsequent studies explored how particle inertia (Saffman Reference Saffman1956; Rosén et al. Reference Rosén, Lundell and Aidun2014; Einarsson et al. Reference Einarsson, Candelier, Lundell, Angilella and Mehlig2015; Di Giusto et al. Reference Di Giusto, Bergougnoux, Marchioli and Guazzelli2024), and the rheology of the carrier medium (e.g. shear-thinning viscosity (Abtahi & Elfring Reference Abtahi and Elfring2019; Li et al. Reference Li, Xu and Zhao2023a ), or viscoelasticity (Leal Reference Leal1975; D’Avino et al. Reference D’Avino, Hulsen, Greco and Maffettone2014; Wang, Yu & Lin Reference Wang, Yu and Lin2019; Li et al. Reference Li, Xu and Zhao2023b ) influence Jeffery orbits. While these studies provided valuable insights into the angular dynamics of spheroids in Couette flows, many practical applications – such as cell focusing techniques (Di Carlo et al. Reference Di Carlo, Irimia, Tompkins and Toner2007) – typically involve Poiseuille duct flows, where particles experience both rotation and lateral migration.
The seminal work of Segré & Silberberg (Reference Segré and Silberberg1961) established that in a dilute suspension of spherical particles with inertia, initially randomly distributed, the particles gradually migrate across streamlines to reach a stable configuration, forming a narrow ring approximately 0.6 times the tube radius from the centre (Segré & Silberberg Reference Segré and Silberberg1962; Matas, Morris & Guazzelli Reference Matas, Morris and Guazzelli2004). Similarly, in square or rectangular duct flows of Newtonian fluids, spheres tend to accumulate along the centre of the duct walls, at a distance of approximately 0.3 times the channel width from the cross-sectional centre (Chun & Ladd Reference Chun and Ladd2006; Di Carlo Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Choi, Seo & Lee Reference Choi, Seo and Lee2011; Miura, Itano & Sugihara-Seki Reference Miura, Itano and Sugihara-Seki2014; Kazerooni et al. Reference Kazerooni, Fornari, Hussong and Brandt2017). In a laminar Newtonian flow, the particle equilibrium position is determined by the competition between two opposing inertial forces: the shear-gradient lift, which drives particles away from the centreline (Ho & Leal Reference Ho and Leal1974; Asmolov Reference Asmolov1999), and the wall-induced lift, which pushes them towards the channel centre (Zeng, Balachandar & Fischer Reference Zeng, Balachandar and Fischer2005; Martel & Toner Reference Martel and Toner2014). This type of migration was later analysed in the context of non-spherical particles (Hur et al. Reference Hur, Choi, Kwon and Carlo2011; Masaeli et al. Reference Masaeli, Sollier, Amini, Mao, Camacho, Doshi, Mitragotri, Alexeev and Di Carlo2012; Uspal, Burak Eral & Doyle Reference Uspal, Burak Eral and Doyle2013; Takeishi & Imai Reference Takeishi and Imai2017). The cross-streamline migration of oblate particles in square ducts resembles that of spheres, exhibiting similar equilibrium positions and manifolds. However, unlike spheres, spheroidal particles tumble during migration, resulting in non-uniform lateral trajectories and longer downstream focusing lengths (Lashgari et al. Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017; Hafemann & Fröhlich Reference Hafemann and Fröhlich2023). Additionally, spheroids with higher aspect ratios (more elongated) tend to have their centres closer to the duct mid-plane at equilibrium (Lashgari et al. Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017).
The rheological properties of the carrier fluid represent another critical factor influencing the dynamics of particle migration (Bartram, Goldsmith & Mason Reference Bartram, Goldsmith and Mason1975; Won & Kim Reference Won and Kim2004; Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007; Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2018; Shaqfeh Reference Shaqfeh2019; Tanriverdi et al. Reference Tanriverdi2025). Even in flows with negligible inertia (
${\textit{Re}} \simeq 0$
), the presence of viscoelasticity can induce cross-streamline migration of particles, unlike in Newtonian fluids where no migration occurs in creeping flows (Gauthier, Goldsmith & Mason Reference Gauthier, Goldsmith and Mason1971; Lormand & Phillips Reference Lormand and Phillips2004; D’Avino et al. Reference D’Avino, Greco and Maffettone2017). Note that the Reynolds number,
$Re=\rho VL/\mu$
, represents the ratio of inertial to viscous forces, where
$\rho$
is the fluid density,
$V$
is a characteristic velocity,
$L$
is a characteristic length, and
$\mu$
is the dynamic viscosity. Indeed, the interplay between inertia and elasticity complicates the particle migration dynamics significantly, as discussed by Li, McKinley & Ardekani (Reference Li, McKinley and Ardekani2015) and Tanriverdi et al. (Reference Tanriverdi, Cruz, Habibi, Amini, Costa, Lundell, Mårtensson, Brandt, Tammisola and Russom2024). In the context of non-spherical particles in viscoelastic media, Tai et al. (Reference Tai, Wang and Narsimhan2020a
,
Reference Tai, Wang and Narsimhanb
) derived approximate analytical expressions for the viscoelastic force and torque acting on arbitrarily shaped particles suspended in a general quadratic flow. These studies, based on a second-order fluid model, were limited to weakly viscoelastic regimes characterised by slow, gradual deformations. The authors found that spheroidal particles tend to align their long axes with the vorticity direction, although finite aspect ratios lead to deviations from predictions for slender bodies. Furthermore, scaling analyses reveal that the lift experienced by prolate particles initially exceeds that of spheres, but diminishes as the particles transition to a log-rolling state near the channel centre. This behaviour ultimately results in a reduced average migration velocity for more elongated particles. Subsequent experimental work by the same authors (Tai & Narsimhan Reference Tai and Narsimhan2022) confirmed these theoretical predictions, notably showing that spherical particles migrate laterally faster than spheroids.
Langella et al. (Reference Langella, Franzino, Maffettone, Larobina and D’Avino2023) experimentally studied the lateral migration of prolate particles in a viscoelastic fluid flowing through a microchannel under creeping flow conditions. At relatively low flow rates, where the particles experience the constant-viscosity regime of the fluid, viscoelasticity drives the spheroids towards the channel centre without inducing a preferential orientation. More recently, Hu et al. (Reference Hu, Lin, Zhu, Yu, Lin and Li2024) investigated numerically the elasto-inertial migration of a neutrally buoyant spheroid in Oldroyd-B viscoelastic duct flows. The authors identified distinct rotational modes of spheroidal particles and demonstrated that, by tuning the elastic number near a critical threshold, spheroids with varying aspect ratios can be effectively separated. Furthermore, Li, Xu & Zhao (Reference Li, Xu and Zhao2024) showed that, unlike spherical particles, oblate spheroids undergo oscillatory cross-stream migration within the duct, driven by orientation-dependent lift forces. They also identified a critical elastic number that controls the migration efficiency of oblate particles in such flows. The elasto-inertial migration of spheroidal particles has been extensively utilised in various microfluidics and biological applications (Liu et al. Reference Liu, Xue, Chen, Shan, Tian and Hu2015; Lu et al. Reference Lu, Zhu, Hua and Xuan2015; Liu et al. Reference Liu, Liu, Semenec, Yuan, Yan, Cain and Li2022).
While much of the existing literature has focused on flows dominated by either inertia or elasticity, the combined influence of elasticity and yield stress remains less explored, despite its potential to produce unique particle behaviours (Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Cheddadi et al. Reference Cheddadi, Saramito, Dollet, Raufaste and Graner2011; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Fraggedakis, Dimakopoulos & Tsamopoulos Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016; Lopez et al. Reference Lopez, Naccache, de Souza and Paulo2018; Izbassarov & Tammisola Reference Izbassarov and Tammisola2020; Pourzahedi, Zare & Frigaard Reference Pourzahedi, Zare and Frigaard2021; Habibi Reference Habibi2026). For example, Chaparian et al. (Reference Chaparian, Ardekani, Brandt and Tammisola2020a ) showed that, in EVP channel flows, a single sphere can penetrate the central plug and reach the centreline at specific Weissenberg and Reynolds numbers, in contrast to purely viscoelastic cases where particles tend to remain off centre. Complementing this, Zade et al. (Reference Zade, Shamu, Lundell and Brandt2020) experimentally demonstrated that inertia modulates migration patterns in Carbopol duct flows: at low inertia, particles cluster near the corners, while at higher inertia, they distribute in a circular pattern between the duct centre and the walls. More recently, Habibi et al. (Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025) investigated the dynamics of spherical particle suspensions in duct flows of EVP fluids via direct numerical simulations. They discovered that increasing the yield stress of the fluid promotes particle migration towards the duct corners in EVP fluids with shear-thinning plastic viscosity, as a result of the yield stress amplifying the first normal stress difference. They also found that particle-induced stresses can disrupt the central plug region, leading to the formation of separated islands of unyielded areas in the EVP flow.
While a limited number of studies have investigated the migration of spherical particles in EVP flows, the complex interplay between the particle dynamics and EVP fluid rheology remains poorly understood. This knowledge gap is even more pronounced for non-spherical particles, which more accurately represent those encountered in practical applications spanning biomedical, industrial and geophysical systems. The anisotropic geometry of these particles introduces additional complexity, giving rise to distinct migration patterns, orientation dynamics and clustering behaviour influenced by both particle shape and the rheological properties of the carrier fluid. To address this, we perform fully resolved, three-dimensional direct numerical simulations (DNS) of spheroidal particles in pressure-driven EVP duct flows. Our simulations systematically investigate the influence of particle shape, inertia, volume fraction and key rheological parameters – including yield stress, elasticity and shear-thinning viscosity – on both the transient evolution and the final spatial distribution and orientation of particles. Additionally, we compare the behaviour of spheroidal and spherical particles to explore potential opportunities for shape-based particle separation in EVP carrier fluids.
The remainder of this paper is organised as follows. The problem formulation and corresponding non-dimensional numbers are outlined in § 2. Next, the governing equations for the fluid and particles, along with the numerical methods, are presented in § 2.1 and § 2.2. The simulation results concerning the migration and orientation of spheroidal particles within EVP fluids are discussed in § 3. Finally, the key conclusions of this study are presented in § 4.
2. Problem statement
We investigate numerically the behaviour of non-colloidal rigid spheroidal particles with an equivalent diameter
$D$
, suspended in a laminar EVP fluid flowing through a straight duct with a square cross-section. The equivalent diameter is defined as
$D =2(ab^2)^{1/3}$
, where
$a$
and
$b$
are the semi-major and semi-minor axes of the spheroids. The particles are either spherical, oblate spheroids (lentil-like), or prolate spheroids (rugby ball-shaped). Simulations are performed in a Cartesian domain of dimensions
$L_x= 2H$
,
$L_y= 6H$
and
$L_z= 2H$
, where
$H$
is the half-height of the square-duct cross-section, and
$y$
,
$z$
and
$x$
denote the streamwise, vertical and spanwise directions, respectively (see figure 1). Particles are initially at rest and randomly distributed and oriented throughout the domain. The blockage ratio, defined as
$\kappa = D/(2H)$
, is set to 0.25 in most cases; additional simulations with
$\kappa = 0.2$
and
$\kappa = 0.125$
are conducted to examine particle size effects. Brownian motion is neglected, as the particles are assumed large enough for thermal fluctuations to be insignificant. In addition, all particles are neutrally buoyant, with densities equal to that of the surrounding fluid. In what follows, the velocity is non-dimensionalised by the bulk flow velocity
$V_b$
, length by the half-height of the channel
$H$
and time by
$H/V_b$
. Pressure and the extra stress tensor are scaled by
$\rho {V_b}^2$
. Note that the total viscosity of the material (
$\mu$
) is the sum of the solvent and polymer viscosities (
$\mu = \mu _s + \mu _p$
). The key non-dimensional parameters governing the system are summarised in table 1, where
$\tau _y$
is the yield stress,
$\lambda$
is the relaxation time of the EVP fluid,
$N$
is the number of particles and
$L_x$
,
$L_y$
and
$L_z$
are the duct dimensions in each direction.
Table 1. Non-dimensional numbers in the problem of spheroidal particle suspension in EVP duct flows.
Figure 1. Three-dimensional view of the computational domain with solid volume fraction
$\phi = 6\,\%$
. Trains of oblate particles are formed at the duct core due to the particle migration towards the centre. Velocity distributions are provided for the streamwise direction (on the
$yz$
-plane) and the
$z$
-direction (on the
$xy$
-plane), with
$N_1$
showing the first normal stress difference in the
$xz$
-plane. Particle colours are for illustrative purposes only.
The computational domain (see figure 1) is discretised using 160
$\times$
480
$\times$
160 Eulerian grid points (for the
$\kappa = 0.2$
) along the spanwise (
$x$
), streamwise (
$y$
) and vertical (
$z$
) directions, respectively. Each spheroidal particle is resolved using 32 Eulerian grid points across its equivalent diameter, with surface discretisation comprising 3720 (oblate), 3610 (prolate) and 3219 (spherical) uniformly distributed Lagrangian points. This level of resolution is chosen in order to properly resolve the fluid–particle interactions. All simulations are conducted at a constant bulk velocity
$V_b$
, which in most cases corresponds to a fixed Reynolds number
${\textit{Re}}=20$
. Periodic boundary conditions are applied to both the fluid and particles in the streamwise direction, while no-slip and no-penetration conditions are enforced on the remaining domain boundaries. While wall slip is a characteristic feature of yield-stress fluid flows (Péméja et al. Reference Péméja, Géraud, Barentin and Le Merrer2019; Younes et al. Reference Younes, Bertola, Castelain and Burghelea2020), the present study assumes a no-slip boundary condition. This choice is motivated by our previous investigations (Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025), which demonstrated excellent agreement with experimental data for particle suspensions in Carbopol under the no-slip assumption (Zade et al. Reference Zade, Shamu, Lundell and Brandt2020). Furthermore, the no-slip condition isolates the influence of bulk rheology and precludes the introduction of additional unconstrained empirical parameters associated with slip models.
2.1. Governing equations
The dynamics of the suspending fluid is described by the incompressible Navier–Stokes equations, subject to mass and momentum conservation
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u} &= -\boldsymbol{\nabla }{p} + \frac {\beta _s}{\textit{Re}}{\nabla} ^2 \boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }^p + \boldsymbol{f}. \\[0pt] \nonumber \end{align}
Here,
$\boldsymbol{u}$
and
$p$
denote the velocity and pressure fields, respectively, while
$\boldsymbol{\tau }^p$
represents the polymeric stress tensor and
${\textit{Re}}$
is the Reynolds number. The term
$\boldsymbol{f}$
on the right-hand side of the momentum conservation accounts for the immersed boundary force, which captures particle-fluid interactions (detailed in § 2.2). The total stress tensor is composed of Newtonian solvent and polymeric contributions:
$\boldsymbol{\tau } = \boldsymbol{\tau }^s + \boldsymbol{\tau }^p$
, where
$\boldsymbol{\tau }^s = \beta _s (\boldsymbol{\nabla }\boldsymbol{u} + {\boldsymbol{\nabla }\boldsymbol{u}}^T$
), and
$\beta _s = \mu _s/\mu$
denotes the viscosity ratio. In all simulations,
$\beta _s$
is fixed at 0.1. To model the non-Newtonian behaviour of the EVP fluid, we employ the Giesekus (Giesekus Reference Giesekus1982) modification of the Saramito equation (Saramito Reference Saramito2007) to capture the second normal stress difference and the shear-thinning viscosity of the EVP carrier fluid
\begin{align} {\textit{Wi}}\, \overset {\triangledown }{\boldsymbol{\tau }^p} + F\left (\boldsymbol{\tau }^p + \frac {\textit{Wi} \,\,\alpha }{\bigl(1-\beta _s\bigr)}\bigl(\boldsymbol{\tau }^p\boldsymbol{\cdot }\boldsymbol{\tau }^p\bigr)\right ) = \frac {2(1 - \beta _s)}{\textit{Re}} \boldsymbol{D}, \end{align}
Here,
$Wi$
denotes the Weissenberg number, defined in table 1, and the function
$F = \max \left(0, ({|\boldsymbol{\tau }^p_d| - Bi/Re})/{|\boldsymbol{\tau }^p_d|} \right)$
accounts for the presence of yield stress in the Saramito model where
${\textit{Bi}}$
denotes the Bingham number. The term
$|\boldsymbol{\tau }^p_d|\equiv \sqrt {\boldsymbol{\tau }^p_d : \boldsymbol{\tau }^p_d /2}$
, where
$\boldsymbol{\tau }^p_d = \boldsymbol{\tau }^p - (\mathrm{tr}\, \boldsymbol{\tau }^p/\mathrm{tr}{\boldsymbol{I}})\, \boldsymbol{I}$
, defines the deviatoric part of
$\boldsymbol{\tau }^p$
, and
$\boldsymbol{I}$
denotes the identity tensor. The shear-thinning behaviour is modelled through the mobility parameter
$\alpha$
with
$0 \lt \alpha \lt 1$
. When
$\alpha = 0$
, the model reduces to the Saramito EVP model, and if both
$\alpha = 0$
and
${\textit{Bi}} = 0$
, it recovers the Oldroyd-B viscoelastic model (Oldroyd Reference Oldroyd1950). Note that
$\overset {\triangledown }{\boldsymbol{\tau }^p}$
represents the upper-convected time derivative of the polymer stress tensor, defined as (Oldroyd Reference Oldroyd1950)
\begin{align} \overset {\triangledown }{\boldsymbol{\tau }^p} \equiv \frac {\partial \boldsymbol{\tau }^p}{\partial t} + u\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }^p -{\boldsymbol{\nabla }\boldsymbol{u}}^T\boldsymbol{\cdot }\boldsymbol{\tau }^p - \boldsymbol{\tau }^p\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}. \end{align}
The translational and rotational dynamics of each particle is governed by the Newton–Euler equations
\begin{align} \rho _p \, V_p \, \frac {\mathrm{d} \boldsymbol u_p}{\mathrm{d} t} &= \oint _{\partial V} \boldsymbol \sigma \boldsymbol{\cdot }\boldsymbol n \,\, \mathrm{d}A + \boldsymbol F_c, \end{align}
\begin{align} \frac {\mathrm{d} \boldsymbol (I_p\boldsymbol \omega _p)}{\mathrm{d} t} &= \oint _{\partial V} \boldsymbol r\times (\boldsymbol \sigma \boldsymbol{\cdot }\boldsymbol n) \,\, \mathrm{d}A + \boldsymbol T_c, \end{align}
where
$\boldsymbol u_p$
and
$\boldsymbol \omega _p$
are the linear and angular velocities of the particle,
$\rho _p$
,
$V_p$
and
$I_p$
are its density, volume and moment of inertia, respectively. The particle domain is denoted by
$\partial V$
and the position vector relative to the particle centre by
$\boldsymbol{r}$
. Note that the Cauchy stress tensor is expressed as
$\boldsymbol{\sigma } = -p\boldsymbol{I} + \boldsymbol{\tau }^p + \beta _s (\boldsymbol{\nabla }\boldsymbol{u} + {\boldsymbol{\nabla }\boldsymbol{u}}^T$
). Collision forces and torques due to particle–particle and particle–wall interactions are represented by
$\boldsymbol{F}_c$
and
$\boldsymbol{T}_c$
. These short-range interactions are modelled using a soft-sphere collision framework with a lubrication correction, following Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015) for spherical and Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016) for spheroidal particles. Further details regarding the lubrication correction and the collision model are provided in Appendix C.
2.2. Numerical method
The numerical framework employs a direct-forcing immersed boundary method (IBM), in which the particles are represented by moving Lagrangian grids interacting with the fluid field defined on a stationary Eulerian mesh. The governing momentum and constitutive equations for the carrier fluid are discretised using a finite-difference scheme. To solve the fluid momentum equations, a projection method is used to decouple pressure and velocity, with the pressure Poisson equation efficiently handled by a scalable Fast Fourier Transform-based solver (Costa Reference Costa2018). Spatial derivatives are primarily computed using central finite differences, except for the advective term in the constitutive equation, which is treated with a fifth-order weighted essentially non-oscillatory scheme (Shu Reference Shu2009). Temporal integration of all governing equations is performed using a third-order Runge–Kutta method. The code is parallelised in two spatial directions using pencil-based domain decomposition and the Message Passing Interface framework, with a manager–worker communication strategy optimising particle data exchange and load balancing (Costa Reference Costa2018).
The fluid–solid interaction is modelled using the direct-forcing IBM algorithm proposed by Breugem (Reference Breugem2012) and later extended to ellipsoidal particles by Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016), incorporating lubrication, friction and collision models to capture short-range particle interactions. In the present study, we adopt this numerical framework to simulate suspensions of oblate and prolate spheroidal particles in EVP duct flows. The simulations are more computationally demanding than those with spherical particles (Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025) due to the need to update each particle moment of inertia at every time step. In the IBM, the fluid domain is discretised using a uniform staggered Cartesian grid (
$\Delta x = \Delta y = \Delta z$
), while the particle surfaces are represented by a set of uniformly spaced Lagrangian markers that move with the interface. The fluid momentum equations are solved over the entire Eulerian grid, and the influence of immersed particles is incorporated by introducing a localised forcing term,
$\boldsymbol{f}$
, on the right-hand side of the momentum conservation equation (2.2). This force satisfies the no-slip and no-penetration boundary conditions at the particle surface by acting near the fluid–solid interface (Mittal & Iaccarino Reference Mittal and Iaccarino2005). The computed immersed boundary force is applied to both the fluid and the particles to advance their motion in time. Further details of the numerical framework can be found in our earlier studies (Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016; Izbassarov et al. Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018; Niazi Ardekani Reference Niazi Ardekani2019). Note that we performed box size and resolution studies to ensure the reported results are free from numerical artefacts (see Appendix A). Additionally, both the IBM and EVP fluid solvers were validated through separate test cases (see Appendix A). Due to the fine resolution and slow particle migration dynamics, simulations are computationally intensive, with some cases requiring up to six weeks on 32 processor cores, depending on the focusing distance needed to reach steady state.
3. Results
This section presents numerical simulations of the spheroidal particle dynamics in EVP duct flows. We first examine results for particle migration and orientation across a range of carrier fluids, including Newtonian, viscoplastic (VP), viscoelastic (VE) and EVP models. These cases establish a baseline within our numerical framework, allowing the specific effects of EVP rheology to be distinguished from simpler constitutive models. In § 3.2, we discuss particle migration in VE duct flows in greater detail. We then extend the analysis to EVP fluids and shear-thinning yield-stress fluids in § 3.3 and § 3.4, respectively. Finally, in § 3.5, we examine the collective behaviour of spheroidal particles in EVP fluids for dilute and semi-dilute suspensions. A summary of the simulation cases performed in this study is provided in table 2.
Table 2. Summary of the non-dimensional parameters explored in our simulations. Domain dimensions
$L_x$
,
$L_y$
and
$L_z$
are expressed in units of the particle diameter for
$\kappa = 0.25$
. For cases with lower
$\kappa$
, the domain dimensions increase proportionally.
3.1. Influence of fluid rheology on lateral migration
We first examine the migration behaviour of oblate particles in Newtonian (
$El = 0,\ Bi = 0$
), VE (
$El = 0.05, \ Bi = 0$
), VP (
$El = 0, \ Bi = 1$
) and EVP (
$El = 0.05, \ Bi = 1$
) duct flows (see table 1 for the definition of elasticity number
$El$
and Bingham number
$Bi$
). The VE fluid is modelled by the Oldroyd-B equation (with
${\textit{Bi}}=0$
and
$\alpha =0$
in (2.3)), while the EVP fluid is modelled using the Saramito constitutive equation (
$\alpha =0$
in (2.3)). The VP fluid is modelled by considering a negligible elasticity in the Saramito model, thus approximating the Bingham plastic model. Figure 2(a) depicts a cross-section of the square duct, with each quadrant corresponding to one of the carrier fluids. Dashed ellipses denote the initial positions and orientations of the particles projected onto the cross-section, and solid ellipses indicate their final states. The lines connecting these two states represent particle trajectories in the duct cross-section (
$xz\hbox{-}$
plane). Insufficient shear stress at the duct centre in yield-stress fluids leads to an unyielded plug, shown as circles in the VP and EVP quadrants. The simulations are conducted at a fixed Reynolds number of
$20$
and a blockage ratio of
$0.25$
. The aspect ratio of the oblate particles is set to 1/3. In addition, the orientation of the oblate particle is described by the unit vector
$\boldsymbol{n} = [n_x, n_y, n_z]$
, aligned with its symmetry axis. Similarly, the rotation of the particle is characterised by the vector
$\boldsymbol{\omega } = [\omega _x, \omega _y, \omega _z]$
. In all simulations, the particle is initially oriented with
$\boldsymbol{n} = [0, 0, 1]$
, so that its symmetry axis is aligned with the vertical direction (
$z$
), while both its initial velocity and rotation are set to zero.
Figure 2. (a) Lateral trajectories of an oblate particle (with an aspect ratio of
$AR = 1/3$
) in a square duct, with each quadrant corresponding to a different carrier fluid: Newtonian, VP, VE and EVP. Circles and triangles indicate initial and equilibrium positions, while dashed and solid ellipses denote initial and final orientations, respectively. (b) Orientation of the oblate symmetry axis along the
$z$
-axis as a function of non-dimensional time. (c) Temporal evolution of the oblate angular velocity along the spanwise direction (
$\omega _x$
). (d) Three-dimensional visualisation of the motion of a single oblate particle in a VE duct flow. (e) Instantaneous snapshot of the flow field with the oblate particle in EVP duct flow. The first normal stress difference,
$N_1$
, is shown in the
$xz$
-plane, while the contours in the
$xy$
- and
$yz$
-planes represent the streamwise velocity component (
$V$
). The central plug region is illustrated in purple.
In the Newtonian case, the particle tends to stay away from the walls due to wall-repulsion forces and migrates slightly towards the duct centre while aligning along the diagonal. This behaviour is consistent with the observations of Lashgari et al. (Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017) regarding the migration of oblate particles with high aspect ratios (
$\kappa \geqslant 0.25$
) in Newtonian duct flows. In the VP case, the particle exhibits negligible migration towards the centre but still aligns along the diagonal. This limited migration results from repulsion from both the duct walls and the central plug. The role of the plug in hindering particle migration towards the centre is discussed in detail in § 3.3. In contrast to the Newtonian and VP fluids, the migration behaviour in the VE and EVP cases is markedly different: the particle migrates towards and ultimately reaches the duct centre.
Figure 2(b) shows the time evolution of the
$z$
-component of the orientation vector,
$n_z$
, for selected cases. The colours correspond to those in panel (a). In the Newtonian case, after an initial transient during which the particle undergoes tumbling motion, it rapidly reaches a steady orientation of
$\boldsymbol{n} = [0.7, 0, 0.7]$
. In the VP case, the particle exhibits pronounced tumbling before its rotation gradually damps, ultimately achieving a steady diagonal orientation similar to the Newtonian case. In contrast, in the EVP case, the orientation dynamics exhibits a long transient and strong tumbling due to interactions with the central plug region. The particle then yields the plug and enters the central region, ultimately becoming fully trapped around
$t \approx 280$
, with a steady orientation of
$n_z = 0.9$
.
Furthermore, figure 2(c) shows the time evolution of the angular velocity of the oblate particle around the spanwise (
$x$
) axis. Since the particles settle in a diagonal orientation in the Newtonian, VP and VE cases,
$\omega _x$
and
$\omega _z$
exhibit similar behaviour, indicating rotation about the shortest axis (spin). After the initial tumbling motion in the VE and EVP cases,
$\omega _x$
decays to zero, indicating that spinning ceases at equilibrium. In contrast, for the Newtonian and VP fluids, steady spinning persists, even though the particle orientation about the
$y$
-axis remains aligned along the duct diagonals (i.e.
$\omega _y = 0$
).
To further examine the spatial distribution of particles within the duct, panels (c) and (d) show three-dimensional (3-D) instantaneous snapshots of the flow field and the oblate particle in the VE and EVP fluids, captured after the particle dynamics reached a steady state. In both cases, the cross-section in the
$xz$
-plane displays the distribution of the first normal stress difference (
$N_1 = \tau _{\textit{yy}} - \tau _{zz}$
), while the
$yz$
- and
$xy$
-planes illustrate contours of the streamwise velocity. In the EVP flow, the central plug region is highlighted in purple, indicating that the oblate particle is trapped within the plug and is transported together with it.
The presence of VE normal stresses is the primary reason for the distinct particle migration observed in VE and EVP fluids compared with non-elastic cases such as Newtonian and VP flows. As depicted in the 3-D snapshots for both VE and EVP flows,
$N_1$
exhibits maxima near the walls (
$N_1 = 2$
) and minima at the centre and corners (
$N_1 = 0$
). This non-uniform distribution generates an asymmetric elastic force that drives particles towards the duct centre (Ho & Leal Reference Ho and Leal1976; Zhou & Papautsky Reference Zhou and Papautsky2020). As a result, particles can overcome inertial forces and migrate towards the duct centre, depending on the local
$N_1$
gradient.
In the following subsections, we present the details of the particle dynamics in VE (§ 3.2) and EVP fluids (§ 3.3 and § 3.4) and elaborate on the role of the governing parameters in particle migration and orientation.
3.2. Viscoelastic carrier fluid
Figure 3. (a) Migration trajectories of oblate particles in VE duct flows with varying fluid elasticity (
$El = 0.012$
to
$0.1$
) shown in the
$xz$
-plane. Dashed ellipses indicate the initial position and orientation, while solid lines mark the final state. (b) Focusing length versus particle spanwise position, where vertical dashed lines indicate the equilibrium positions. (c) Vertical particle position over time for different elasticity values. For ease of comparison particles start from the same vertical position.
Figure 4. (a) Cross-sectional trajectories of oblates with aspect ratios ranging from 1/4 to1, confinement ratios from 0.125 to 0.25 and different initial positions. The inset of the figure shows the vertical particle position over time for various particle shapes. (b) Time evolution of the angular velocity of the oblate with
$AR = 1/3$
along the streamwise (
$y$
) and spanwise (
$x$
) directions at
$El = 0.01{-}0.1$
.
Here, we examine the migration behaviour of oblate particles in a VE duct flow, modelled by the Oldroyd–B equation (
${\textit{Bi}} = 0$
and
$\alpha = 0$
). Figure 3(a) shows a cross-section of the square duct, including the initial and final particle positions and orientations, as well as the particle trajectories in the duct cross-section (
$xz$
-plane). Simulations are conducted at a fixed Reynolds number of
${\textit{Re}}=20$
and a blockage ratio of 0.25, employing oblate particles with an aspect ratio of 1/3.
To investigate the effect of fluid elasticity on particle migration and orientation, the elasticity number (
$El$
) is varied systematically between
$0.01$
and
$0.1$
, with each quadrant in figure 3(a) corresponding to one of these elasticity values. In all cases, the particles migrate towards the duct core and align with the diagonal symmetry line within their respective quadrants. In carrier fluids with low elasticity (
$El\lt 0.05$
), particles migrate only partially towards the duct centre, settling at intermediate positions between the centre and the corners. In contrast, at higher elasticities (
$El \geqslant 0.05$
), particles complete their migration to the duct core. Figure 3(b) shows the focusing length (
$y/H$
) – defined as the downstream distance required for oblates to reach their equilibrium position – plotted against their initial lateral offset from the duct centreline. Trajectories are colour coded consistently with panel (a). For the initial positions considered here, lower elasticity appears to yield shorter focusing lengths (e.g. the particle in
$El = 0.0125$
focuses within approximately 180 channel heights, compared with roughly 800 and 300 for
$El = 0.05$
and
$El = 0.1$
, at the same bulk Reynolds number). This trend is primarily due to the proximity of the release points to the respective equilibrium positions; the ordering of the focusing lengths may reverse if particles are instead introduced near the centreline.
To investigate this further, we plot the vertical position of oblate particles versus time in figure 3(c). The results reveal that particles in fluids with lower elasticity (
$El = 0.0125$
and
$0.025$
) reach equilibrium more quickly, as they have a partial migration and travel a shorter distance to their equilibrium position. In contrast, at higher elasticity, particles fully migrate to the centreline, with faster migration observed in the
$El = 0.1$
case compared with
$El = 0.05$
.
This migration behaviour appears robust across various initial positions, particle sizes and aspect ratios, as illustrated in figure 4(a). The figure shows oblates with aspect ratios ranging from 1/4 (highly elongated) to 1 (spherical), all having the same equivalent diameter of 1 and suspended in a fluid with fixed elasticity of
$El=0.05$
. Owing to the geometric symmetry of the square duct, initial particle positions are sampled across only one eighth of the cross-sectional area. Regardless of initial conditions, all particles ultimately migrate towards the duct centre. The spherical particle follows a linear path (black line in figure 4
a) because its initial position lies on the diagonal, and the geometric symmetry of the configuration preserves a straight trajectory. In contrast, non-spherical particles break this symmetry through their rotational and tumbling motions, resulting in the oscillatory lateral trajectories observed in figure 4(a) (Li et al. Reference Li, Xu and Zhao2024). This behaviour is further quantified in the inset of figure 4(a), which plots the vertical position of various particle shapes over time. Among them, the sphere demonstrates the fastest lateral migration, whereas the prolate particle with
$AR=3$
migrates the slowest. This behaviour is due to the ongoing reorientation of oblate and prolate particles along Jeffery-like orbits (Jeffery Reference Jeffery1922), causing trajectory and velocity oscillations that delay their migration to the centre. Note that at the initial stage of motion, the oblate particle with
$AR = 1/2$
and the prolate particle experience faster migration (higher lift) compared with the spherical particle, although overall, the spherical particle exhibits the highest migration velocity among them. Irrespective of shape, the particles ultimately align to minimise their projected surface area in the streamwise direction. For example, the oblates (
$AR = 1/2$
or 1/3) adopt a diagonal orientation to reduce drag.
The angular velocity of the oblate particles with
$AR = 1/3$
around the streamwise axis (
$\omega _y$
) for different elasticity numbers is shown in figure 4(b). The particles exhibit Jeffery-like orbits, with oscillations in angular velocity that gradually decay to zero, indicating convergence to a stable orientation aligned with the duct diagonal. While any orientation with the minor axis lying in the
$xz$
-plane results in the same projected area in the streamwise direction, the diagonal alignment is preferred because it maximises the distance between the particle surface and the duct walls, reducing wall-induced hydrodynamic interactions (Lashgari et al. Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017). In this way, the particle simultaneously minimises drag and avoids near-wall confinement, leading to the observed equilibrium orientation. Comparing
$\omega _y$
for different elasticity numbers shows that the oscillation amplitude decreases as fluid elasticity increases. For example, the peak
$\omega _y$
drops from 0.078 at
$El = 0.012$
to 0.053 at
$El = 0.1$
. This reduction is due to VE stresses that dampen particle rotation, consistent with the findings of Snijkers et al. (Reference Snijkers, D’Avino, Maffettone, Greco, Hulsen and Vermant2011) for spherical particles. Furthermore, figure 4(b) shows the angular velocity of oblate particles around the spanwise axis (
$x$
). Since particles settle in a diagonal orientation,
$\omega _x$
and
$\omega _z$
show similar behaviour, indicating rotation around the shortest axis (spin). For
$El \geqslant 0.05$
,
$\omega _x$
decays to zero, meaning spinning stops at equilibrium. However, for
$El \lt 0.05$
, steady spinning continues, even though the orientation about the
$y$
-axis remains stable. Similar to the focusing length in figure 3(b), oblates in low-elasticity fluids (
$El \lt 0.05$
) reach a steady orientation faster than in high-elasticity cases (
$El \geqslant 0.05$
). Among the latter, alignment occurs notably quicker at
$El = 0.1$
than at
$El = 0.05$
.
In the following subsections, we investigate how yield stress, shear-thinning viscosity and volume fraction of particles change the migration dynamics in EVP fluids.
3.3. Elastoviscoplastic fluid
Figure 5. Particle migration in an EVP fluid with
$El = 0.025$
. (a) Trajectories of oblate particles (
$AR = 1/3$
) in Saramito duct flows for varying Bingham numbers (
${\textit{Bi}} = 0$
–1). (b) Streamwise angular velocity (
$\omega _y$
) and spanwise angular velocity (
$\omega _x$
) of oblate particles as a function of time. (c) Vertical position of particles over time.
Figure 6. Particle migration in EVP fluids with (a)
${\textit{Bi}} = 1$
and varying elasticities
$El = 0.01-0.05$
, (b) trajectories of oblate particles in Saramito duct flows for varying Bingham numbers (
${\textit{Bi}} = 0$
–1), at
$El=0.05$
, (c) streamwise velocity profiles for varying
${\textit{Bi}}$
and
$El = 0.05$
. Far-field profiles (symbols) are compared with profiles across the particle centre at
$z/(2H)=0.5$
(lines); filled circles mark the particle position. The inset provides a magnified view of the central region for clarity. (d) Particle final equilibrium positions as a function of
${\textit{Bi}}$
and
$El$
. Cases marked with star symbols indicate instances of unsteady angular dynamics.
Now, we shift our focus to EVP carrier fluids, which, in addition to VE stresses, exhibit a yield stress. Figure 5(a) illustrates the final position and orientation of oblate particles in the cross-section of an EVP duct flow, simulated using the Saramito model with elasticity
$El = 0.025$
, and Bingham numbers varying from 0 to 1. Owing to the lack of sufficient shear stress to yield the material at the duct centre, an unyielded region (plug) forms at the duct core. This plug region is illustrated in the figure for various Bingham numbers, demonstrating its expansion with increasing yield stress. Additionally, figure 5(b) shows the angular velocities of the particles about the
$x$
and
$y$
axes, and figure 5(c) presents the particle vertical position over time for each Bingham number.
In the case with
${\textit{Bi}} = 0$
, the oblate particle migrates towards the duct centre, reaching an equilibrium position situated between the duct walls and the centre, but biased towards the centre. The oblate aligns diagonally within the duct due to the configuration of minimal energy dissipation in shear flows (see § 3.2 for more discussion), and the angular velocity along the streamwise axis (
$\omega _y$
) approaches zero. Crucially, the values of
$\omega _x$
remain non-zero, indicating that the oblate undergoes a spinning motion while maintaining its diagonal alignment. Increasing the Bingham number of the suspending fluid up to
${\textit{Bi}}=0.5$
preserves similar steady-state particle behaviour, with final position and orientation consistent with that in VE fluids. However, the particle angular velocities in the
$y$
- and
$x$
-directions exhibit significantly more pronounced oscillations and tumbling during its migration towards the steady state in the EVP fluid.
At
${\textit{Bi}} = 1$
, however, the particle angular dynamics undergoes a fundamental change. The oblate particle now exhibits time-dependent rolling and tumbling motions around an equilibrium position situated on the diagonal symmetry line. Consequently, a stable, steady-state orientation or rotation rate cannot be defined at
${\textit{Bi}} = 1$
and
$El=0.025$
. One representative orientation of the oblate is illustrated with a dotted line in the figure 5(a). The particle angular velocity also displays a distinct behaviour: the angular velocity along the
$y$
-axis no longer attains a steady-state value but exhibits quasi-periodic oscillations (see panel (b)). Similarly, the angular velocities along the
$x$
-axis display quasi-periodic behaviour. This gyroscopic-like motion can be attributed to the substantial size of the unyielded region (plug) at the duct centre at
${\textit{Bi}} = 1$
where the plug comes in contact with the oblate. To clarify, the particle experiences continuous collisions with this plug, which exerts a varying force on the oblate. Analogous to a gyroscope, the oblate responds by precessing and nutating, thereby resisting changes to its angular momentum in the flow. The magnitude of the angular velocities in this case does not decay over time. Additionally, figure 5(c) shows the vertical position of the particle as a function of time for varying Bingham numbers, indicating that, for
${\textit{Bi}}=0{-}1$
, the particle centre of mass will eventually reach the same equilibrium position. However, as the Bingham number increases, the migration process becomes slower, and the particle trajectory exhibits more pronounced oscillations which delay the migration. Notably, the vertical motion appears to remain periodic at
${\textit{Bi}}=1$
.
Turning to the influence of the elasticity number on the observed peculiar dynamics at
${\textit{Bi}}=1$
, figure 6(a) presents the final positions and orientations of oblate particles suspended in an EVP fluid with
${\textit{Bi}} = 1$
, for
$El = 0.0125$
to
$0.05$
. The corresponding particle trajectories, projected onto the duct
$xz\hbox{-}$
cross-sectional plane, are also shown. At
$El = 0.0125$
and
${\textit{Bi}} = 1$
, the particle approaches the duct core but remains outside the plug, settling into a diagonal orientation. For
$El = 0.025$
, the particle interacts with the plug region, exhibiting an unsteady, quasi-periodic angular dynamics (as in figure 5
a). However, at the highest elasticity (
$El=0.05$
), the particle migrates further towards the duct centre, penetrating the plug region to reach the duct core. It is noteworthy to note that, although the plug region does not fully encompass the particles, it can nonetheless trap them, preventing rotation and causing the particles to travel at the same velocity as the central plug. In this case, the final orientation of the particle within the plug region is determined by its initial position and orientation, leading to a distribution of potential angles due to variations in the particle entry trajectory. Chaparian et al. (2020a) observed centreline migration of a small sphere (
$\kappa = 0.1$
) under moderate elasticity and
${\textit{Re}} = 200$
, despite the presence of a plug region. They attributed this to particle rotational inertia inducing local yielding in the surrounding material, enabling deeper plug penetration. Conversely, in our case, the particle remains partially outside the central plug, so the
$N_1$
gradient from the yielded material still actively drives it towards the duct core, resulting in immediate trapping within the plug.
To further examine this effect, figure 6(b) illustrates the steady-state particle configuration in a Saramito fluid, maintained at the highest elasticity value of
$El = 0.05$
and for Bingham numbers ranging from
${\textit{Bi}} = 0$
to 1. For each Bingham number, the central plug zone is also depicted. In all cases, the particles attain an equilibrium position at the duct core. At a moderate Bingham number (
${\textit{Bi}} \lt 0.5$
), the small plug zone yields a dynamics similar to the VE case (
${\textit{Bi}} = 0$
), with diagonal alignment. In contrast, at higher Bingham numbers (
${\textit{Bi}} \geqslant 0.5$
), the larger relative plug size alters the particle final orientation, which now depends on its initial position and orientation at entry to the plug zone. As observed in the figure 6(b), increasing the yield stress intensifies oscillations in the particle
$xz\hbox{-}$
trajectory compared with the VE migration. As we will show, this enhanced oscillatory motion delays particle migration towards the duct centre in EVP fluids. Figure 6(c) illustrates the steady streamwise velocity profiles (
$V$
) for EVP flows at varying
${\textit{Bi}}$
. The unperturbed far-field profiles are shown as solid lines, while symbols represent the perturbed profiles at the particle’s cross-section. The final particle positions are indicated by filled circles at the duct centre. In the VE baseline (
${\textit{Bi}}=0$
), the far-field profile exhibits a Poiseuille duct flow distribution. In the presence of the particle, the local velocity profile is modified, exhibiting a plateau at the centreline corresponding to the particle’s translational velocity. Since the particle equilibrates at the centreline without rotating, the flow disturbance remains localised to its immediate vicinity, and the velocity field recovers rapidly to its far-field state. Notably, the particle lags the flow, as its translational velocity is lower than the unperturbed centreline fluid velocity. This observation is consistent with Li et al. (Reference Li, McKinley and Ardekani2015), who reported that, at relatively large blockage ratios, central particle migration in Oldroyd-B flows results in a negative slip velocity. In cases with a non-zero Bingham number, a central plug region forms even in the far-field profiles where the shear stress is insufficient to yield the material. Notably, in these cases, the velocity profile across the particle centre identically matches the far-field profile. This occurs because the particle is fully trapped within the central plug; since the relative velocity between the particle and the surrounding medium is zero, the particle induces no additional velocity disturbances in the flow field (see the 3-D visualisation in figure 2
e).
Figure 6(d) illustrates the equilibrium position
$Z_{eq}$
of oblate particles (
$AR =$
1/3) in an EVP duct flow with a confinement ratio of 0.25, as a function of elasticity (
$El$
) and Bingham number (
${\textit{Bi}}$
). All particles are initially at rest, positioned at
$Z_{eq} = 3$
and oriented along the
$x$
-axis. The colour of each point denotes the final equilibrium height within the duct cross-section. For
$El \geqslant 0.05$
, particles migrate fully to the duct centre (even in the presence of the central plug region), while for
$ 0.012 \leqslant El \lt 0.05$
, only partial migration to the duct centre occurs, with equilibrium positions between the centre and the walls (see figure 6
a). An interesting observation is that, at low elasticity numbers, particles remain further away from the centre with increasing
${\textit{Bi}}$
, suggesting that yield stress hinders particle migration towards the duct centre. This behaviour can be attributed to two underlying factors. First, the emergence of a central plug region reduces the effective flow area, thereby increasing shear rates and enhancing inertial forces that drive particles towards the walls and corners (Chaparian et al. Reference Chaparian, Ardekani, Brandt and Tammisola2020a
). Additionally, increasing yield stress in EVP flows promotes particle migration towards the corners by affecting the distribution and gradients of the
$N_1$
, even in the absence of a central plug region as shown by Habibi et al. (Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025) for spherical particles. As a result, increasing yield stress slightly displaces particles away from the duct centre, even at fixed elasticity. In figure 6(d), cases exhibiting unsteady motion – characterised by persistent tumbling and wobbling – are marked with a hexagonal symbol. In each case, the particle comes in contact with the central plug region, leading to an unsteady rotational dynamics (see figure 5 for details).
To further understand particle migration in EVP fluids, the competing forces on the particle surface must be considered. In these fluids, as in VE media, cross-streamline migration is mainly driven by the gradient of
$N_1$
(Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025). Following Ho & Leal (Reference Ho and Leal1976), the elastic lift force on a particle of diameter
$D$
scales with the first normal stress difference as:
$ F^p_e = -({40}/{3}) \pi \rho U_c^2 D_{eq}^2 \kappa El \,(1 - \beta _s) Z_p$
, where
$Z_p$
is the dimensionless vertical offset and
$U_c$
is the centreline velocity. In contrast, Nizkaya et al. (Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020) introduced a general form for the inertial lift force acting on an oblate particle as:
$ F^p_i = \rho \, C_l \, a^3 b \, \dot {\gamma }^2$
, where
$C_l$
is the lift coefficient. By balancing inertial and elastic forces in EVP fluids, we estimate a critical elasticity number
$El_c \approx 0.01$
, above which oblate particles with
$ AR = 1/3$
migrate towards the centreline (see figure 6
a). This threshold includes both partial migration – where the particle is biased towards the centre but remains between the centreline and wall – and complete migration to the duct core. This value aligns with previous findings for spherical particles in VE (Li et al. Reference Li, McKinley and Ardekani2015) and EVP fluids (Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025) and for oblates with
$AR=1/2$
in VE media (Li et al. Reference Li, Xu and Zhao2024). As the elasticity of the fluid increases, the resulting lift force strengthens, allowing the particles to overcome the yield-stress barrier and enter the plug region. Our simulations for
${\textit{Bi}} = 0$
–1 show that an elasticity number of
$El = 0.05$
is sufficient to drive particles fully into the duct core (complete migration). Additional effects of particle shape in EVP fluids are discussed in Appendix B.
3.4. Effect of shear thinning in elastoviscoplastic fluids
Up to this point, we have examined the effects of viscoelasticity and yield stress on the oblate particle migration. However, many VE and EVP materials in nature and industries exhibit shear-thinning behaviour. To this end, we employ the Saramito–Giesekus model to simulate the carrier fluid, as it captures viscoelasticity, yield stress and shear-thinning effects within a unified framework (see § 2.1). It also accounts for weak secondary flows in a duct flow (Li et al. Reference Li, McKinley and Ardekani2015; Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025).
Figure 7. Effect of shear-thinning viscosity in oblate migration in an EVP fluid with
$El = 0.05$
and varying mobility parameters (
$\alpha$
). (a) Trajectories of oblate particles suspended in fluids with
$\alpha = 0-0.3$
. In the case of
$\alpha = 0.1$
, particles are released from various initial positions. (b) Temporal evolution of the spanwise (
$\omega _x$
) and streamwise (
$\omega _y$
) angular velocities of oblates. (c) Time evolution of the
$z$
component of the orientation axis of the oblate particle.
Figure 7(a) illustrates particle migration and orientation in an EVP fluid across a range of shear-thinning plastic viscosities, from
$\alpha = 0$
(no shear-thinning viscosity) to 0.3 (highest shear-thinning viscosity), with an elasticity of 0.05 and a Bingham number of 1. The central plug region for each fluid is also depicted. Notably, the introduction of shear thinning significantly expands the plug size compared with the case with constant plastic viscosity (
$\alpha = 0$
). In the non-shear-thinning EVP case (
$\alpha = 0$
), as described before, the particles migrate to the centre of the duct. However, in the shear-thinning cases
$\alpha \gt 0$
, the particles exhibit a surprising behaviour, moving away from the duct centre (plug zone) towards the walls and settling in an equilibrium position between the centre and the corners, while orienting perpendicular to the diagonal symmetry lines. We further simulated two additional cases with varying initial positions within the quadrant of
$\alpha = 0.1$
. The results indicate that even when released near the central plug, the particle migrates away from the duct centre and equilibrates at the same lateral position between the centre and the walls. Due to symmetry, a particle released exactly on the symmetry line (
$x=2$
) follows a straight trajectory toward the duct walls while orienting itself perpendicular to the
$z\hbox{-}$
axis. Migration of oblate particles in shear-thinning fluids is notably rapid; within fewer than 100 non-dimensional time units for
$\alpha = 0.1$
and 50 time units for
$\alpha = 0.2$
and 0.3, the particle angular velocity decays to zero (see figure 7
b), and it adopts a stable, perpendicular orientation. In contrast to Newtonian, VE and Saramito-type fluids – where particles undergo sustained Jeffery-like oscillations before settling – the shear-thinning case features only a single oscillation prior to reaching equilibrium. Notably, unlike in Saramito fluids, where particles tend to align along the duct diagonals while spinning, in Saramito–Giesekus fluids they remain perpendicular to them, exhibiting neither spin nor tumbling motion. Panel (c) shows the time evolution of the
$z-$
component of the orientation axis for the oblate particle. We compare the case of
$\alpha = 0$
with cases involving shear thinning. The colour coding in this plot is consistent with that used in panel (a). As observed, the particle reaches a stable orientation significantly faster (
$t \approx 50$
) in the presence of shear-thinning plastic viscosity. In contrast, for the
$\alpha = 0$
, the particle migrates towards the duct centre and penetrates the plug region, requiring a much longer time to achieve a steady-state orientation (
$t \approx 280$
).
Habibi et al. (Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025) analysed the physical mechanism governing the behaviour of spherical particles in Saramito–Giesekus fluids. Their findings show that the shear-thinning effect – introduced by the nonlinear term
$ {\textit{Wi} \alpha }/{(1-\beta _s)}(\boldsymbol{\tau }^p \boldsymbol{\cdot }\boldsymbol{\tau }^p)$
– reduces elastic stresses and diminishes centre-directed migration. Consequently, particles preferentially migrate towards the walls and corners when
$\alpha \gt 0.1$
. Additionally, shear thinning lowers the apparent viscosity at high shear rates, steepening velocity gradients near the walls and flattening the velocity profile. This amplifies the shear-gradient lift force – scaling as
$F_L \sim \rho U^2 a^6 / H^4$
(Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009) – further promoting wall-directed migration (Villone et al. Reference Villone, D’avino, Hulsen, Greco and Maffettone2013; Li et al. Reference Li, McKinley and Ardekani2015; Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025). As a result of these two mechanisms, particles in Saramito–Giesekus fluids migrate towards the duct walls rather than the centre, in contrast to the Saramito case.
A novel observation in the present study is the stable orientation of oblate particles perpendicular to the duct diagonals in the Saramito–Giesekus suspension. Orienting the side with the smaller surface area in the streamwise direction reduces the surface exposed to the flow and is therefore energetically favourable; however, this does not explain why the particle avoids diagonal alignment seen in VE and Saramito EVP cases. However, strong elastic and repulsive forces from adjacent walls (see figure 3 d) inhibit diagonal orientations with particle tips directed into corners. Consequently, particles stabilise in orientations perpendicular to the diagonals, balancing minimal flow resistance with wall avoidance.
To examine the influence of yield stress on the orientation of oblate particles in shear-thinning fluids, we vary
${\textit{Bi}}$
from 0 (purely VE) to 1, as shown in figure 8(a). In the VE and weak yield-stress cases (
${\textit{Bi}} \lt 1$
), particles partially migrate towards the duct centre, while at
${\textit{Bi}} = 1$
, they slightly shift towards the duct corners. This indicates that although shear-thinning effects tend to drive particles towards the corners (through weakening of the elastic effects), the distinct orientation perpendicular to the duct diagonals requires a sufficiently large central plug region. Figure 8(b) shows the vertical position of the particle over time, revealing that increasing yield stress pushes the equilibrium position further from the centre, suggesting a stronger influence of inertia. While particles reach positional equilibrium quickly in both VE and EVP fluids, at
${\textit{Bi}} \lt 1$
their angular motion exhibits long oscillations, gradually damping to zero. In contrast, the
${\textit{Bi}} = 1$
case rapidly reaches zero angular velocity within 50 time units. The perpendicular orientation of the oblate relative to the duct diagonals remains stable, as it aligns with the axis of maximum moment of inertia towards the walls, and the shear-stress asymmetry at its tips is insufficient to induce rotation. Conversely, when the oblate is oriented along a diagonal (as in cases with
${\textit{Bi}}\lt 1$
), a pronounced stress imbalance develops across its sides, leading to sustained particle spin at equilibrium. To further investigate this phenomenon, panel (c) presents the streamwise velocity profiles of the Saramito–Giesekus fluids across a range of
${\textit{Bi}}$
. As observed, the spatial extent of the plug region expands significantly with increasing yield stress. Since the bulk velocity
$V_b$
is held constant, the reduction of the centreline velocity at higher
${\textit{Bi}}$
must be compensated by a lateral expansion of the 3-D velocity paraboloid. This redistribution of the flow field steepens the velocity gradients in the near-wall regions for flows with higher
${\textit{Bi}}$
. This effect is captured in the inset, where a magnified view of the velocity profile at the duct corner is shown; here, the dashed lines representing the wall velocity tangents for
${\textit{Bi}}=0$
and
${\textit{Bi}}=1$
clearly illustrate the sharpening of the gradient. As discussed earlier, these intensified velocity gradients enhance the inertial lift forces – which scale with the local shear rate (Ho & Leal Reference Ho and Leal1974; Li et al. Reference Li, McKinley and Ardekani2015) – thereby driving the particles further away from the duct centre towards the walls.
Figure 8. Effect of Bingham number on oblate migration in an EVP fluid with
$El=0.05$
and
$\alpha =0.2$
. (a) Trajectories of oblate particles suspended in fluids with
${\textit{Bi}} = 0-1$
, and (b) time evolution of the particle vertical position and the spanwise angular velocity (
$\omega _x$
). (c) Streamwise mean velocity profiles of Saramito–Giesekus fluids for various
${\textit{Bi}}$
. The inset provides a magnified view of the near-wall region, where dashed lines represent the velocity tangents. (d) Mapping of particle equilibrium positions and orientations in the
${\textit{Bi}}$
$El$
-plane for Saramito–Giesekus fluids with
$\alpha = 0.2$
. The vertical equilibrium position is represented by the colour bar (right). Orientation modes are distinguished by symbols: squares for normal-to-diagonal alignment, circles for diagonal alignment and stars for cases exhibiting an unsteady angular dynamics.
To conclude this subsection, figure 8(d) illustrates the equilibrium position (
$Z_{eq}$
) of oblate particles (
$AR = 1/3$
) in Saramito–Giesekus duct flows with a confinement ratio of 0.25, plotted as a function of
$El$
and
${\textit{Bi}}$
for
$\alpha = 0.2$
. In these simulations, all particles were released from rest at
$Z=3$
with an initial orientation normal to the
$z-$
axis. The data points are colour coded to represent the final equilibrium height within the duct cross-section. Across the entire range of elasticity values investigated, the particles migrate away from the walls towards the duct corners. This lateral migration is further enhanced by increasing elasticity, which drives the particles towards the corner regions – the areas of minimum shear rate and elastic stress. This observation aligns with the findings of Habibi et al. (Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025), who demonstrated that increased elasticity in Saramito–Giesekus fluids laden with spherical particles promotes migration towards the duct corners. Notably, the particles maintain a greater lateral distance from the duct centre with increasing
${\textit{Bi}}$
. This suggests that in Saramito–Giesekus fluids – analogous to the behaviour observed in Saramito fluids – the yield stress effectively promotes migration towards the corners by augmenting the inertial lift forces. The angular dynamics of the particles is also mapped in the same figure. Depending on the specific combination of
$El$
and
${\textit{Bi}}$
, the particles either align with the duct diagonals (denoted by circular symbols) or orient themselves normal to the diagonals (square symbols). At low
${\textit{Bi}}$
, the particles favour a diagonal orientation; this represents the most energetically favourable state as it minimises drag and maximises the distance from the confining duct walls (Lashgari et al. Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017). At intermediate
${\textit{Bi}}$
, the particles exhibit unsteady angular behaviour, characterised by persistent tumbling and wobbling, similar to regimes observed in Saramito fluids (cf. figure 5
a). In these cases, the particle is situated in close proximity to, or in direct contact with, the central plug region; the resulting contact torque with the rigid plug leads to this unsteady rotational dynamics. With a further increase in yield stress, this transient unsteady regime eventually transitions into a stable state where the particle remains oriented normal to the diagonal.
Figure 9. Instantaneous snapshots of the flow field (left column) and mean particle concentration
$\varPhi (x,z)$
(right column) across the duct cross-section are shown for (a) the Saramito EVP fluid with
$El = 0.05$
,
${\textit{Bi}} = 1$
, and (b) the Saramito–Giesekus fluid with
$El = 0.05$
,
${\textit{Bi}} = 1$
and
$\alpha = 0.2$
. In all cases, the solid volume fraction is
$\phi = 6\,\%$
and the Reynolds number is
${\textit{Re}} = 20$
. The 3-D contours in the flow field snapshots depict the first normal stress difference (
$N_1$
) and the streamwise velocity (
$v$
).
3.5. Collective behaviour of oblate particles
In this subsection, we examine the migration and orientation behaviour of suspensions of oblate particles and compare this dynamics with those of spherical particles. We first investigate a suspension of oblate particles with an aspect ratio of 1/3, at a solid volume fraction of
$\phi =6\,\%$
, flowing through a duct characterised by a blockage ratio
$\kappa =0.25$
. Figure 9 displays instantaneous three-dimensional snapshots of the flow field and particle distribution at statistically steady states for two carrier fluids: a Saramito fluid (figure 9
a) and a shear-thinning Saramito–Giesekus fluid (figure 9
b). Each snapshot illustrates the distributions of the first normal stress difference (
$N_1$
) and the streamwise velocity (
$v$
) on their respective planes. Complementing these visualisations, we present the probability distribution function of the concentration,
$\varPhi (x,z)$
, representing the time and spatially averaged solid phase distribution across the duct cross-section. All simulations are conducted at a Reynolds number
${\textit{Re}}=20$
. For the EVP model, the elasticity number is fixed at
$El=0.05$
and the Bingham number at
${\textit{Bi}}=1$
. It is worth emphasising that the Saramito model assumes a constant polymer viscosity (
$\alpha = 0$
), while the Saramito–Giesekus model incorporates shear-thinning behaviour through a non-zero mobility factor (
$\alpha = 0.2$
); recall the governing equations in § 2.1.
Figure 10. Time and spatially averaged distribution of yielded and unyielded regions (
${\textit{YR}}$
) in a suspension of oblate particles within an EVP fluid modelled using the Saramito–Giesekus formulation, together with an instantaneous 3-D view of unyielded regions (
${\textit{YR}}=0$
) and particle positions. The
$YR$
map shows the yielding probability at each location, where
$YR=0$
indicates a solid-like behaviour and
$YR=1$
indicates a fully yielded region.
In the Saramito carrier fluid, particles migrate towards the duct centre, penetrating the central plug region and forming a concentrated cluster of oblate particles. Within this cluster, the particles exhibit a broad range of angular orientations without clear preferential alignment. This contrasts with the single-particle simulations presented previously in this work, where an isolated particle adopts a stable diagonal orientation at the centreline. The disparity is attributed to the collective dynamics, where particle–particle interactions disrupt the alignment seen in isolated cases and give rise to a wider distribution of orientations.
In contrast to the Saramito suspension, the Saramito–Giesekus fluid drives particles towards the walls and corners of the duct, leading to the formation of distinct chains of oblate particles with preferred orientations nearly perpendicular to the diagonal, as in the single-particle case. While the majority of particles migrate from the duct centre towards the walls, a subset remains near the centre. To explain this, we examine the unyielded regions of the flow shown in figure 10, which includes a 3-D snapshot of the unyielded zones together with the time-averaged yielded region (
${\textit{YR}}$
) across the duct cross-section. The
${\textit{YR}}$
map represents the probability of the material being yielded at each location, with
${\textit{YR}} = 0$
indicating fully unyielded, solid-like behaviour while with
${\textit{YR}} = 1$
represents fully yielded region. At the duct centre, where shear stress is minimal, the material remains unyielded, forming a central plug that moves at the maximum flow velocity. As illustrated in the figure, while the majority of particles are driven towards the corners, those initially located near the centre can become entrained within the plug. These particles remain effectively stationary relative to the plug and are passively advected at the plug velocity. Trapping particles were observed also in the single-particle case in previous sections, but only for the Saramito fluid.
Another notable phenomenon observed in the Saramito–Giesekus simulations is the emergence of secondary flows, attributed to the combined effects of the second normal stress difference (
$N_2$
) and the non-circular duct geometry (Yue, Dooley & Feng Reference Yue, Dooley and Feng2008). The presence of particles further amplifies these flows by introducing disturbances into the velocity field. As shown in figure 11(a), the secondary flow field in the EVP suspension at a particle volume fraction of
$\phi = 6\,\%$
consists of eight vortical structures. The peak magnitude of these secondary flows reaches approximately
$4 \times 10^{-3}$
. No secondary flows are observed at the duct centre, as this region is occupied by the unyielded plug, where the EVP material behaves as a solid and thus cannot support shear-driven vortical motion. In Saramito–Giesekus flows, the secondary flow magnitude exceeds the lateral migration velocity of the particles (
$O(10^{-3})$
), thereby driving the particles further towards the duct corners.
Figure 11. (a) Secondary flows in EVP suspension with
${\textit{Bi}} = 1$
,
$El = 0.05$
and
$\alpha = 0.2$
. For both cases
$\phi = 6\,\%$
and
${\textit{Re}} = 20$
. (b) The difference between the time-averaged velocities of the fluid and particles, normalised by the bulk flow velocity. White regions indicate locations where no particles are present.
In addition to secondary flows, figure 11(b) presents the time-averaged streamwise velocity profiles of the percentage difference between the fluid and oblate particle velocities (
$V_f {-} V_p$
), relative to the bulk velocity in the EVP suspension, plotted across the duct section. As shown, the particle velocity in the duct core matches that of the plug (zero velocity difference), confirming that particles enclosed within the plug are carried along with it. Further away from the centreline, however, particle average velocities slightly lag behind the fluid (approximately 11 % of the bulk velocity), primarily due to particle inertia. However, there exists a narrow region where the fluid velocity falls below that of the particles.
Figure 12. Effect of particle shape on migration in an EVP medium: (a) instantaneous snapshots of particle distributions and flow fields in Saramito–Giesekus fluids for oblates with
$AR = 1/3$
and spheres, with planes indicating
$N_1$
values; (b) mean particle concentration profiles for aspect ratios
$AR = 1/3$
,
$1/2$
and 1.
Finally, we compare the migration behaviour of oblate and spherical particles. Figure 12 depicts 3-D snapshots of suspensions containing 6 % volume fraction of either oblate or spherical particles in a Saramito–Giesekus duct flow. All rheological parameters are identical between the two cases, with particle shape being the only difference. To further illustrate particle distribution, the time and spatially averaged concentration
$\varPhi (y,z)$
is also presented for the oblate particles with
$AR =$
1/3 and 1/2 and the spherical suspensions. This quantity represents the mean solid phase concentration across the duct cross-section. In all cases, particles migrate from the duct centre towards the corners. Spherical particles settle at the corners, while oblate particles reach equilibrium slightly away from them, collectively forcing a small number of particles to be trapped in the centre. Among the oblate cases, particles with an aspect ratio of 1/3 remain farther from the wall compared with those with
$AR=1/2$
. This difference arises from the higher shear gradients at the oblate tips, which enhance wall repulsion and the first normal stress difference (note the green circle in figure 12
a). These increased forces, particularly near the sharper tips of oblates, prevent them from reaching the corners. Due to their anisotropic shape and larger effective surface area when aligned perpendicular to the diagonals, oblates experience stronger wall-induced repulsion. This interaction, especially in confined geometries, outweighs migration forces and causes oblates to stabilise near – but not at – the corners. Conversely, isotropic spheres experience symmetric wall forces and, under strong VE and shear-thinning effects in the Saramito–Giesekus fluid, can settle stably at the corners where shear and elastic force are minimal (Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025). This shape-dependent behaviour is particularly relevant for microfluidics applications, as it enables particle separation based on shape and aspect ratio. Such an approach can complement existing sorting techniques, including inertial microfluidics (Masaeli et al. Reference Masaeli, Sollier, Amini, Mao, Camacho, Doshi, Mitragotri, Alexeev and Di Carlo2012).
4. Concluding remarks
A flowing suspension of particles in EVP fluids presents a complex dynamics due to the interplay of multiple mechanisms influenced by both particle geometry and the rheological properties of the fluid. Given the importance of EVP fluids in both natural and industrial contexts, their detailed study is essential. To this end, we performed 3-D, interface-resolved DNS to investigate the dynamics of rigid spheroidal particles in EVP duct flows. In all cases, inertial effects were included. We employed the Saramito and Saramito–Giesekus models to capture the EVP behaviour of the fluid, resolving particle–fluid interactions using the IBM, along with lubrication corrections and a soft-sphere collision model to account for inter-particle and particle–wall interactions.
We initially examined the migration of single spheroidal particles in VE fluids, then extended our study to incorporate yield-stress effects characteristic of EVP fluids. We further analysed suspensions of oblate particles and identified key differences between individual and collective behaviours. The main findings are summarised below.
-
(i) In Saramito-modelled EVP fluids, particles migrate towards the duct centre when the elasticity number exceeds 0.01, consistent with VE cases, due to the gradient of the first normal stress difference (
$N_1$
). For elasticity numbers above 0.05, particles are able to fully enter the plug region, moving at the same velocity as the plug and without rotation relative to it. We also demonstrated that increasing yield stress at fixed elasticity causes oblates to shift towards the duct walls and corners, in contrast to previous findings about spherical particles. -
(ii) The angular dynamics of particles is strongly influenced by both elasticity and yield stress. At moderate elasticity (
$El \simeq 0.025$
) and high yield stress (
${\textit{Bi}} \geqslant 1$
), the plug region becomes comparable in size to the particle, resulting in partial contact and inducing gyroscopic motion to the particles. As a result, particles interacting with the plug exhibit unsteady, periodic spinning and tumbling. In contrast, particles that fully enter the plug (
$El \geqslant 0.05$
) adopt its velocity and cease rotation, with their final orientation strongly influenced by their entry angle, determined by their initial position and orientation. At lower elasticity (
$El \lt 0.02$
), particles align with the duct diagonals between the centre and corners without interacting with the plug, resembling VE behaviour. -
(iii) In Saramito–Giesekus fluids, particles tend to migrate away from the centre towards the corners due to weakened elastic forces and enhanced wall-directed inertial effects from the shear-thinning viscosity. Nevertheless, for sufficiently large plug regions (
${\textit{Bi}} \geqslant 1$
), oblate particles stay away from the corners and align perpendicularly to the duct diagonals. In such cases, particles rapidly reach a steady state with negligible angular velocities, deviating from classical Jeffery orbits. When the plug is smaller, particles move towards the wall along a diagonal trajectory and spin, but do not exhibit tumbling. -
(iv) The dynamics of dilute and semi-dilute suspensions largely mirrors that of single particles. In Saramito fluids, particles accumulate in the plug region and display a wide range of orientations due to collision-induced variability. In Saramito–Giesekus fluids, particles form trains near duct corners, aligned perpendicularly to the diagonals and moving slightly faster than the local fluid. Secondary flows also emerge in Saramito–Giesekus cases, with magnitudes comparable to the particles lateral migration velocity, further promoting their migration towards the corners.
-
(v) Finally, we compared suspensions of oblate particles (aspect ratios 1/3 and 1/2) with spherical particles (
$AR = 1$
). Spherical particles accumulate at the duct corners, overcoming wall-repulsion forces. In contrast, oblate particles experience stronger wall-repulsion and elastic forces, causing them to settle slightly away from the corners. This distance increases with decreasing aspect ratio (i.e. for more flattened particles). These shape-dependent interactions suggest a promising route for particle separation in microfluidics applications, where the complex flow behaviour of EVP fluids can be harnessed to sort particles by shape.
Acknowledgements
The authors acknowledge the computer time provided by SNIC and NAISS (National Academic Infrastructure for Super-computing in Sweden).
Funding
This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 955605 YIELDGAP. We gratefully acknowledge the support of European Research Council through Starting Grant MUCUS (Grant No. ERC-StG-2019-852529).
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Code validation
A.1. Resolution and box size studies
The numerical solver used in this study has been successfully applied and validated in various configurations, including the sedimentation of rigid particles in EVP media (Sarabian et al. Reference Sarabian, Rosti, Brandt and Hormozi2020), the dynamics of soft particles and droplets in EVP flows (Izbassarov et al. Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018; Izbassarov & Tammisola Reference Izbassarov and Tammisola2020), EVP flow through porous structures (De Vita et al. Reference De Vita, Rosti, Izbassarov, Duffo, Outi and Brandt2018; Chaparian et al. Reference Chaparian, Izbassarov, De Vita, Brandt and Tammisola2020b ) and turbulent EVP flows (Rosti et al. Reference Rosti, Izbassarov, Tammisola, Hormozi and Brandt2018; Izbassarov et al. Reference Izbassarov, Rosti, Brandt and Tammisola2021). Additionally, Habibi et al. (Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025) demonstrated that our numerical results for spherical particle suspensions in EVP duct flows exhibit excellent agreement with existing experimental velocity profiles for finite-size particles in Carbopol duct flows.
For the sake of completeness, we present four additional validation tests. The first assesses the sensitivity of the results to domain size and grid resolution. The baseline simulation is conducted in a square duct of size
$2h \times 6h \times 2h$
, discretised using a
$160 \times 480 \times 160$
grid, corresponding to a resolution of 32 grid points per particle diameter. To evaluate the influence of streamwise periodicity on the simulation results, we extend the domain length to
$8h$
while keeping the grid resolution unchanged, resulting in a
$160 \times 640 \times 160$
mesh (denoted as ‘Long’). We also consider a high-resolution case (referred to as ‘Res40’) with the same domain size but a refined grid of
$200 \times 600 \times 200$
, corresponding to 40 Eulerian grid points per particle diameter and 5036 Lagrangian markers on the particle interface.
In all configurations, we track a particle (
$(2h)/D_p=4$
) in a Saramito fluid at
${\textit{Re}}=20$
,
${\textit{Wi}}=1$
, and
${\textit{Bi}}=1$
. We observe that the particles, initially positioned at
$Z_p, X_p = 3$
, gradually migrate towards the duct centreline. The results, shown in figure 13, demonstrate a great agreement in lateral velocity, angular velocity, migration trajectory and equilibrium position of particles across all test cases. This consistency supports the use of the reference case (‘Res32’) for the main simulations.
Figure 13. Comparison of the temporal variation of (a) vertical migration velocity (
$W_p$
), (b) angular velocity around
$x$
-axis, (c) angular velocity around
$z$
-axis and (d) vertical position of the particle (
$Z_p$
). Simulations correspond to Saramito model with
${\textit{Re}} = 20$
,
${\textit{Wi}} = 1$
,
${\textit{Bi}} = 1$
and
$\kappa = 0.25$
.
A.2. Immersed boundary method solver validation
The second validation case assesses the angular velocity of a spheroidal particle in a Couette flow, comparing simulation results with Jeffery’s analytical solution (Jeffery Reference Jeffery1922) for aspect ratios
$AR = 1/3$
and 1/2. The particles are neutrally buoyant suspended in the middle of a plane Couette flow at particle Reynolds number
${\textit{Re}}_p =$
0.1. Here,
${\textit{Re}}_p$
is defined by shear rate
$\dot {\gamma }$
and the equivalent particle diameter
$D_{eq}$
as
${\textit{Re}}_p = \rho \dot {\gamma } D_{eq}^2/\mu$
.
Simulations are conducted in a cubic domain of size
$10D_{eq} \times 10D_{eq} \times 10D_{eq}$
using resolution of 32 grid points per
$D_{eq}$
. Periodic boundary conditions are applied in the directions perpendicular to the velocity gradient, while no-slip and no-penetration conditions are enforced along the wall-normal direction. The initial particle orientation is aligned with the
$x$
-direction, with no initial angular velocity. As expected, the particles exhibit periodic tumbling motion around the spanwise axis, which is normal to the shear plane. The results, presented in figure 14(a), show excellent agreement between our DNS simulations and Jeffery’s analytical solution.
Figure 14. (a) Angular velocity of ellipsoids with aspect ratios
$AR = 1/2$
and
$1/3$
. The DNS results from the present study are compared with Jeffery’s analytical solution (Jeffery Reference Jeffery1922). The particle Reynolds number is
${\textit{Re}}_p = 0.1$
. A schematic of an oblate particle suspended in Couette flow is shown at the bottom of panel (a). (b) Angular velocity of a sphere in a VE Couette flow as a function of the Weissenberg number (
${\textit{Wi}}$
), with
${\textit{Re}}_p = 0.025$
. Blue circles (labelled ‘DNS’) denote results from the present numerical solver. The bottom of panel (b) displays a schematic of a sphere suspended in VE Couette flow.
To further validate the IBM and VE solvers, we benchmark our numerical method against a canonical problem involving the angular velocity of a neutrally buoyant sphere suspended at the centre of a 3-D Couette flow of an Oldroyd-B fluid. Experimental measurements by Snijkers et al. (Reference Snijkers, D’Avino, Maffettone, Greco, Hulsen and Vermant2011) and DNS by Goyal & Derksen (Reference Goyal and Derksen2012) characterised the particle rotation as a function of the Weissenberg number,
${\textit{Wi}} = \lambda \dot {\gamma }$
, varying from 0 to 2. The Couette cell has dimensions
$4R \times 8R \times 8R$
and is discretised with 24 grid points per particle diameter. Opposing wall velocities
$\pm U$
generate a shear rate
$\dot {\gamma } = 2U/8R$
, with periodic boundary conditions applied in the remaining directions. The particle experiences rotation solely about the
$z$
-axis, and the particle Reynolds number is fixed at 0.025 (creeping flow regime).
As shown in figure 14(b), our results (shown by ‘DNS’ label) closely match the experimental data, capturing the observed reduction in angular velocity with increasing
${\textit{Wi}}$
, consistent with the known effect of elasticity in suppressing particle rotation.
A.3. EVP solver validation
In this subsection, we assess the accuracy of the EVP solver by considering planar Couette flow of a Saramito–Giesekus fluid. The material is initially at rest (
$t = 0$
), after which a constant shear rate
$\dot {\gamma }_0$
is applied. The Weissenberg and Bingham numbers are defined as
${\textit{Wi}} = \lambda \dot {\gamma }_0$
and
${\textit{Bi}} = \tau _0 / (\mu \dot {\gamma }_0)$
, respectively. In 2-D simple shear flow, the tensorial form of the constitutive equation (see (2.3)) reduces to the following set of coupled ordinary differential equations:
\begin{align} {\textit{Wi}}\, \frac {\mathrm{d} \tau _{\textit{xx}}}{\mathrm{d} t} - 2{\textit{Wi}}\, {\tau _{\textit{xy}}} + F\left ({\tau _{\textit{xx}}} + \frac {\textit{Wi} \,\,\alpha }{(1-\beta _s)}\big({\tau _{\textit{xx}}}^2 + {\tau _{\textit{xy}}}^2\big)\right )&= 0, \end{align}
\begin{align} {\textit{Wi}}\, \frac {\mathrm{d} \tau _{\textit{yy}}}{\mathrm{d} t} + F\left ({\tau _{\textit{yy}}} + \frac {\textit{Wi} \,\,\alpha }{(1-\beta _s)}\left({\tau _{\textit{xx}}}{\tau _{\textit{xy}}} + {\tau _{\textit{xy}}}{\tau _{\textit{yy}}}\right)\right )&= 0, \end{align}
\begin{align} {\textit{Wi}}\, \frac {\mathrm{d} \tau _{\textit{xy}}}{\mathrm{d} t} - {\textit{Wi}}\, {\tau _{\textit{yy}}} + F\left ({\tau _{\textit{xy}}} + \frac {\textit{Wi} \,\,\alpha }{(1-\beta _s)}\big({\tau _{\textit{xy}}}^2 + {\tau _{\textit{yy}}}^2\big)\right )&= 1-\beta _s. \end{align}
Here,
$F = \max (0, {|\boldsymbol{\tau }_d| - Bi}/{|\boldsymbol{\tau }_d|} )$
and
$|\boldsymbol{\tau }_d|^2 = 1/2(\boldsymbol{\tau }_d : \boldsymbol{\tau }_d)$
. The initial condition for all the stress components is fixed at zero (
$\boldsymbol{\tau }(t=0) = \boldsymbol{0}$
). The above equations are solved numerically for the case of
${\textit{Re}} = 0.05$
,
${\textit{Bi}} = 0.2$
,
${\textit{Wi}} = 1$
,
$\beta _s = 1/9$
and
$\alpha = 0.2$
. Figure 15 illustrates the temporal evolution of the shear stress (
$\tau _{\textit{xy}}$
) and the first normal stress difference (
$N_1$
) obtained by our numerical code in comparison with the semi-analytical solution obtained using the framework in Saramito (Reference Saramito2007, Reference Saramito2009). The shear stress initially rises, leading to a tiny overshoot, and then it falls until it eventually reaches a steady value. It is noteworthy that the Saramito model captures the stress overshoot as it is a phenomenon that is observed in the startup shear of polymeric liquids. The reader can find more details in Schweizer & Xie (Reference Schweizer and Xie2018) and Tseng (Reference Tseng2021) who discuss the origin of the stress overshoot in the startup Couette flow of polymeric liquids from the perspective of the polymer microstructure, and how constitutive equations may capture it.
Figure 15. Evolution of (a) the first normal stress difference and (b) shear stress in time for a planar Couette flow with
${\textit{Re}} = 0.05$
,
${\textit{Bi}} = 0.2$
,
${\textit{Wi}} = 1$
and
$\alpha = 0.2$
. The solid line represents the semi-analytical solution, while the square symbols are the result of the present numerical code.
Figure 16. Inertial migration of an oblate spheroid (
$AR=1/2$
) in a Newtonian planar channel flow at
${\textit{Re}}=22$
and
$\kappa =0.3$
. (a) Particle position and orientation in the
$yz$
-plane, illustrating the transition from the initial release point to the lateral equilibrium position. Contours represent the streamwise velocity field. (b) The
$xz$
-plane showing the steady-state orientation, where the symmetry axis aligns with the vorticity direction (
$\boldsymbol{n}=[1,0,0]$
). (c) Quantitative comparison of the wall-normal trajectory
$z/2H$
as a function of streamwise displacement (
$y/(2H)$
) against the result of Nizkaya et al. (Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020).
A.4. Migration of an oblate spheroid in a Newtonian channel flow
To conclude the validation study, we examine the lateral migration of an oblate spheroid in a Newtonian planar channel flow and compare our results with existing literature. The spheroid is characterised by an aspect ratio of
$AR=1/2$
and a confinement ratio of
$\kappa =0.3$
. The channel Reynolds number is set to
${\textit{Re}}=22$
. Initially, the particle is released from rest at a wall-normal position of
$z_0/(2H)=0.2$
with an orientation axis
$\boldsymbol{n}=(0.75,0.66,0)$
. Periodic boundary conditions are employed in both the streamwise (
$y$
) and spanwise (
$x$
) directions, while a no-slip condition is imposed on the channel walls located at the
$z\hbox{-}$
axis. Figure 16(a) illustrates the initial and equilibrium configurations of the spheroid in the
$yz\hbox{-}$
plane. The background contours represent the streamwise velocity field (
$V$
). Driven by wall-induced lift forces, the particle migrates away from the boundary towards an equilibrium position located between the channel centreline and the wall. This behaviour is consistent with the classical Segré–Silberberg effect observed in the inertial migration of non-spherical particles, as previously characterised by Lashgari et al. (Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017) and Nizkaya et al. (Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020). The equilibrium state in the
$xz\hbox{-}$
plane is shown in panel (b). The spheroid eventually aligns its axis of symmetry with the vorticity direction (the
$x\hbox{-}$
axis), resulting in a steady-state orientation of
$\boldsymbol{n}=(1,0,0)$
. This alignment perfectly matches the steady-state orientations reported by Nizkaya et al. (Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020) for similar flow conditions. Finally, panel (c) provides a quantitative comparison of the particle’s trajectory, plotting the wall-normal position as a function of the downstream streamwise coordinate. Our numerical results demonstrate great agreement with the data from Nizkaya et al. (Reference Nizkaya, Gekova, Harting, Asmolov and Vinogradova2020), accurately capturing both the transient migration phase and the final lateral equilibrium position.
Appendix B. Effect of particle shape in EVP fluids
Figure 17. Particle migration in an EVP fluid with
$El = 0.05$
and
${\textit{Bi}} = 1$
for different particle shapes. (a, b) Spanwise (
$x$
) and vertical (
$z$
) trajectories of oblate, spherical and prolate particles in a Saramito–Giesekus fluid. The insets display final particle positions in the duct cross-section, overlaid on contours of square root of the trace of conformation tensor representing polymer chain stretch. (c) Cross-sectional trajectories of the same particles in a duct flow of an EVP fluid modelled by the Saramito constitutive equation.
Figure 17 presents the migration of particles with varying shapes within an EVP fluid, modelled using the Saramito–Giesekus (panels (a) and (b)) and Saramito (panel (c)) equations with an elasticity of 0.05 and
${\textit{Bi}} = 1$
. All particles have an equivalent diameter of 1. In the case of the Saramito–Giesekus fluid (a,b), shear-thinning effects drive the particles towards the duct walls. The oblate particle adopts an orientation nearly perpendicular to the duct diagonal, as detailed in § 3.4. The prolate particle migrates even closer to the wall compared with the oblate case, aligning with the flow direction to minimise its frontal area. The spherical particle is also displaced towards the upper wall, eventually making contact due to wall-directed inertial forces overcoming the wall-repulsion effects.
In Saramito EVP fluids, particles of all shapes tend to migrate towards the duct centre, ultimately entering the central plug region. This inward migration is driven by elastic effects that exceed a critical elasticity threshold (
$El = 0.05$
for the given confinement ratio). Oblate particles with aspect ratios of 1/2 and 1/3 exhibit a dynamics comparable to that shown in figure 6(a, b), and their final orientations are therefore omitted for brevity. Spherical particles become fully centred within the plug, aligning closely with its width. Prolate particles also migrate towards the core but, once inside the plug, their rotational motion ceases, and they are advected downstream at a constant velocity. In contrast to VE fluids, where particles consistently orient with their smallest surface area facing the flow to minimise drag, orientation in EVP fluids is not uniquely prescribed. Instead, it depends on both the particle initial position and orientation, which influence the angle at which it enters the plug region. For example, while a prolate particle in a VE fluid aligns with its minimal cross-sectional area facing the flow, in Saramito EVP fluids – as illustrated in panel (c) – it adopts a different, non-minimal orientation.
Appendix C. Collision model and lubrication correction
Here, we summarise the key aspects of the collision model and the lubrication correction, following the formulations presented in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015) and Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016).
C.1. Collision model
To account for solid–solid interactions, we employ a soft-sphere collision model proposed by Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015). This model is based on a linear spring–dashpot system acting in both normal and tangential directions, integrated with a Coulomb friction slider. While originally formulated for spherical particles, the model is extended to spheroids following the approach of Niazi Ardekani (Reference Niazi Ardekani2019). Near the contact point, each spheroid is locally approximated by a sphere whose radius corresponds to the local radius of curvature at the point of closest approach. This equivalent sphere is assigned the same mass as the original spheroid, allowing the use of analytical contact mechanics derived for spherical bodies. The normal contact force,
$\boldsymbol{F}_{{\textit{ij}},n}$
, acting on particle
$i$
due to contact with particle
$j$
, is directed along the line of centres
$\boldsymbol{n}_{ij} = (\boldsymbol{x}_j - \boldsymbol{x}_i) / \|\boldsymbol{x}_j - \boldsymbol{x}_i\|$
. This force is determined by the overlap distance
$\delta _{{\textit{ij}},n}$
and the relative normal velocity
$\boldsymbol{u}_{{\textit{ij}},n}$
at the contact point
\begin{equation} \boldsymbol{F}_{{\textit{ij}},n} = -k_n \delta _{{\textit{ij}},n} - \eta _n \boldsymbol{u}_{{\textit{ij}},n}. \end{equation}
The stiffness coefficient
$k_n$
and the damping coefficient
$\eta _n$
are determined by modelling the collision as a linear harmonic oscillator (Van der Hoef, van Sint Annaland & Kuipers Reference Van der Hoef, van Sint Annaland and Kuipers2004). To maintain numerical stability and allow the fluid solver to adapt to the change in particle momentum, the collision time
$T_n$
is set as a fixed multiple of the simulation time step,
$T_n = N \Delta t$
. Following the sensitivity analyses of Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015), we set
$N=8$
. The spring and dashpot coefficients are expressed as
\begin{align} k_n &= \frac {m_e(\pi ^2 + \ln ^2 e_{n,d})}{(N\Delta t)^2}, \end{align}
\begin{align} \eta _n &= -\frac {2m_e \ln e_{n,d}}{N\Delta t}, \end{align}
where
$m_e = (m_i^{-1} + m_j^{-1})^{-1}$
is the reduced mass and
$e_{n,d}$
is the dry coefficient of restitution and is set to 0.97 in this study.
The tangential contact force,
$\boldsymbol{F}_{{\textit{ij}},t}$
, and torque,
$\boldsymbol{T}_{{\textit{ij}},t}$
, is modelled using a similar spring–dashpot framework to represent the resistance to sliding and rolling. To distinguish between sticking and sliding regimes, the tangential force is regularised by a Coulomb friction law
\begin{equation} \boldsymbol{F}_{{\textit{ij}},t} = \min \left ( \| -k_t \delta _{{\textit{ij}},t} - \eta _t \boldsymbol{u}_{{\textit{ij}},t} \|, \| -\mu _c \boldsymbol{F}_{{\textit{ij}},n} \| \right ) \boldsymbol{t}_{ij}, \end{equation}
where
$\mu _c$
is the friction coefficient (fixed at 0.15),
$\delta _{{\textit{ij}},t}$
is the tangential displacement and
$\boldsymbol{t}_{ij}$
is the tangential unit vector. The tangential coefficients
$k_t$
and
$\eta _t$
are obtained by solving the harmonic oscillator equations under the requirement that the tangential and normal collision times match (
$T_t = T_n$
) (Costa et al. Reference Costa, Boersma, Westerweel and Breugem2015)
\begin{equation} k_t = \frac {m_{e,t}(\pi ^2 + \ln ^2 e_{t,d})}{(N\Delta t)^2}, \quad \eta _t = -\frac {2m_{e,t} \ln e_{t,d}}{N\Delta t}, \end{equation}
where
$e_{t,d}$
is the dry tangential coefficient of restitution, which, like the normal coefficient
$e_{n,d}$
, is a prescribed input parameter of the model and is set to 0.1. The effective tangential reduced mass,
$m_{e,t}$
, accounts for the rotational inertia of the colliding bodies and is defined as
\begin{equation} m_{e,t} = \left (1 + \frac {1}{K^2}\right )^{-1} m_e, \end{equation}
where
$K$
denotes the normalised radius of gyration for the approximating spheres (
$K^2 = 2/5$
for a homogeneous solid sphere).
Whenever the tangential force exceeds the Coulomb limit, the displacement is saturated to ensure physical consistency with the sliding regime (Costa et al. Reference Costa, Boersma, Westerweel and Breugem2015).
The total collision force
$\boldsymbol{F}_i^c$
and torque
$\boldsymbol{T}_i^c$
exerted on particle
$i$
are computed as the sum of all individual contact contributions and are added to the right-hand side of the Newton–Euler equations (see (2.5) and (2.6)) to resolve the particle motion. Particle–wall interactions are treated as the limiting case of a collision with a sphere of infinite radius. In this limit, the reduced mass simplifies to
$m_e = m_i$
and the overlap
$\delta _{iw,n}$
is defined relative to the planar wall surface. For further details regarding the collision model, the reader is referred to the work of Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015).
C.2. Lubrication model
A particle immersed in a viscous fluid experiences significant lubrication forces when approaching another body, such as a wall or another particle. These forces, which arise from the drainage of the fluid film within the gap, can be described by the analytical asymptotic solutions of Jeffrey (Reference Jeffrey1982) for spheres of different radii in the Stokes regime. While the IBM inherently resolves these hydrodynamic interactions, it underestimates the lubrication resistance when the gap width falls below the Eulerian grid resolution,
$\Delta x$
(Breugem Reference Breugem2010). To avoid the high computational cost of mesh refinement, we apply a lubrication correction to compensate for the under-resolved pressure in the narrow gap. The dominant lubrication component is the normal squeezing force, which scales with
$1/\epsilon$
(where
$\epsilon$
is the normalised gap width). In contrast, the tangential and rotational components diverge only logarithmically (
$\ln \epsilon$
) and have a negligible impact on the overall contact dynamics (Simões Costa Reference Simões Costa2017).
Similar to the collision model, spheroidal particles are locally approximated near the contact point by spheres with a radius equal to the local radius of curvature (Niazi Ardekani Reference Niazi Ardekani2019). The normal squeezing force is the dominant lubrication component, as it scales with
$1/\epsilon$
, whereas tangential and rotational terms diverge only logarithmically (
$\ln \epsilon$
) and were found to have negligible influence on contact dynamics (Costa et al. Reference Costa, Boersma, Westerweel and Breugem2015).
For gap widths below a resolution threshold,
$\epsilon \lt \epsilon _{\Delta x}$
, a correction force
$\boldsymbol{F}_{\textit{lub}}$
is incorporated into the Newton–Euler equations according to Jeffrey (Reference Jeffrey1982)
\begin{equation} \boldsymbol{F}_{\textit{lub}} = -6 \pi \mu R_i \boldsymbol{u}_{{\textit{ij}},n} \left [ \lambda (\kappa ,\varepsilon ) - \lambda (\kappa ,\varepsilon _{\Delta x}) \right ], \end{equation}
where
$R_i$
is the local radius of curvature,
$\kappa$
is the curvature ratio and
$\lambda$
is the Stokes amplification factor. To prevent the unphysical divergence of the force as
$\epsilon \to 0$
, the amplification is saturated at a roughness length scale
$\epsilon _\sigma$
, representing contact through surface asperities:
$\lambda (\kappa , \epsilon \lt \epsilon _\sigma ) = \lambda (\kappa , \epsilon _\sigma )$
(Costa et al. Reference Costa, Boersma, Westerweel and Breugem2015).
In this study, the resolution thresholds are prescribed as
$\epsilon _{\Delta x} = 0.025$
for particle–particle interactions and
$\epsilon _{\Delta x} = 0.05$
for particle–wall interactions. The value of
$\varepsilon _{\Delta x}$
is determined by simulating a sphere slowly approaching towards a planar wall or between two spheres, and identifying the gap below which the IBM no longer accurately resolves the lubrication interaction. The lubrication correction defined in (C7) is activated only when the gap falls below these specific separations. To account for surface asperities and prevent the unphysical divergence of the lubrication force at contact, a constant roughness cutoff of
$\varepsilon _\sigma = 0.001$
is applied to both particle–particle and particle–wall interactions, following the regularisation proposed by Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015).
Open access
