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