1. Introduction
Microfluidics has transformed the landscape of particle manipulation at small scales, offering precise control over suspended objects such as cells, vesicles, droplets and particles. Among the diverse strategies developed, inertial microfluidics has emerged as a robust and passive technique that relies solely on hydrodynamic forces intrinsic to the flow, rather than externally applied electromagnetic, thermal, optical or acoustic fields (Martel & Toner Reference Martel and Toner2014). This approach has attracted significant attention for its ability to control particle trajectories in confined geometries at intermediate Reynolds numbers, where inertial effects begin to influence motion in a non-negligible manner. Inertial focusing arises from a balance of lift forces generated by velocity gradients and wall interactions, which drive particles toward distinct equilibrium positions within a channel cross-section (Di Carlo et al. Reference Di Carlo, Irimia, Tompkins and Toner2007; Seo, Lean & Kole Reference Seo, Lean and Kole2007). The absence of external actuation simplifies device design and operation, while also making these systems inherently scalable for portable applications. Over the past two decades, inertial microfluidics has been successfully applied to focus, separate, enrich and trap a wide variety of particles (Ateya et al. Reference Ateya, Erickson, Howell, Hilliard, Golden and Ligler2008; Xuan, Zhu & Church Reference Xuan, Zhu and Church2010), with particular impact on biomedical diagnostics, chemical analysis and environmental monitoring (Nilsson et al. Reference Nilsson, Evander, Hammarström and Laurell2009; Pratt et al. Reference Pratt, Huang, Hawkins, Gleghorn and Kirby2011). While many studies have emphasised biological applications such as label-free cell sorting, the same principles of hydrodynamic migration remain equally relevant for the transport and organisation of passive particles in complex fluids. Most previous investigations have focused on Newtonian suspending media, where inertial lift forces are well characterised and equilibrium positions are well established. However, many real-world fluids, including blood, mucus and polymeric solutions, exhibit shear-rate-dependent viscosity, viscoelasticity, anisotropy and yield-stress behaviour, all of which can strongly influence particle–fluid interactions (Lu et al. Reference Lu, Liu, Hu and Xuan2017; Yuan et al. Reference Yuan, Zhao, Yan, Tang, Alici and Zhang2018). In this study, we focus on the shear-thinning behaviour, in which viscosity decreases with increasing shear rate. In such environments, the inertial lift forces are strongly modified by the local rheology, leading to a migration dynamics that can differ substantially from that in Newtonian systems. Understanding these effects is essential not only for advancing the fundamental physics of particle migration but also for improving the design of microfluidic systems handling non-Newtonian fluids in biomedical and industrial applications. This study addresses this gap by examining the inertial migration of passive spherical particles in shear-thinning fluids.
The lateral migration of rigid particles in Newtonian fluids under finite inertia has been extensively investigated over the past six decades. Seminal experiments by Segre & Silberberg (Reference Segre and Silberberg1961, Reference Segre and Silberberg1962a ,Reference Segre and Silberberg b ) revealed that neutrally buoyant spheres in pressure-driven pipe flows migrate toward an annular equilibrium region located around 0.6 times the pipe radius from the pipe centreline, moving further towards the wall with increasing channel Reynolds number. This phenomenon, absent in the creeping-flow limit, was attributed to inertial effects (Bretherton Reference Bretherton1962). Early theoretical analyses in unbounded shear flows, including those by Rubinow & Keller (Reference Rubinow and Keller1961) and Saffman (Reference Saffman1965), predicted lift forces arising from particle rotation and velocity slip, but typically failed to reproduce the experimentally observed annular equilibria, instead suggesting migration toward the centreline. These discrepancies were resolved by subsequent work highlighting the dominant contributions of shear gradients and wall-induced hydrodynamic interactions to the lateral lift force (Martel & Toner Reference Martel and Toner2014). Hence the slip spin force of Rubinow & Keller (Reference Rubinow and Keller1961) and the slip shear force of Saffman (Reference Saffman1965) are considered weak forces (Matas, Morris & Guazzelli Reference Matas, Morris and Guazzelli2004) in the absence of external non-hydrodynamic force/torque. The wall-directed shear-gradient-induced lift force and the centre-directed wall-induced lift-force balance each other at the equilibrium position observed in the experiment by Segre & Silberberg (Reference Segre and Silberberg1961). Analytical progress was made under weak inertial limits through Cox’s perturbation analysis (Cox & Brenner Reference Cox and Brenner1968) and the influential work by Ho & Leal (Reference Ho and Leal1974), who developed lift-force expressions for plane Poiseuille and Couette flows. These were later generalised to three-dimensional ducts by Hood, Lee & Roper (Reference Hood, Lee and Roper2015), offering accurate predictions of off-centre equilibria at moderate Reynolds numbers. Numerical simulations also played a key role in validating and extending these models. Notably, Feng, Hu & Joseph (Reference Feng, Hu and Joseph1994) demonstrated that, as Reynolds number increases, the particle equilibrium shifts toward the wall, and Matas, Morris & Guazzelli (Reference Matas, Morris and Guazzelli2009) identified additional equilibrium branches near the centreline, attributed to finite-size and confinement effects for a particle in pressure-driven flows.
A closer look into the existing literature shows that foundational work on inertial migration has predominantly focused on pressure-driven flows, while the inertial migration of particles in confined linear shear flows, which can be experimentally accessible via rheometric devices, has received comparatively less attention. This configuration, however, is crucial for developing a fundamental understanding of inertial migration mechanisms, as it allows the role of shear-induced and wall-induced effects to be examined in a controlled manner. Earlier studies have provided fundamental insights into inertial lift forces acting on particles in shear flows. A seminal contribution in this context is the work of Ho & Leal (Reference Ho and Leal1974), who employed perturbation theory to analyse the inertial migration of a neutrally buoyant sphere in plane Couette (linear shear) flow between parallel walls. Their study demonstrated that weak inertial lift arises from the interaction between shear-induced disturbances and the confining walls, thereby establishing the fundamental mechanism governing lateral migration in bounded shear flows. Using perturbation theory, Drew (Reference Drew1988) extended Saffman’s analysis to the case of a sphere translating in shear flow near a distant wall at small particle Reynolds number. Their calculations showed that the wall-induced lift force is always smaller than the corresponding unbounded lift force, and hence does not alter the migration direction. The key parameter in these studies is the particle Reynolds number, defined as
${\textit{Re}}_{\!p} = U_{\!p} a / \nu$
, where
$U_{\!p}$
is the particle velocity,
$a$
is the particle radius and
$\nu$
is the kinematic viscosity of fluid. Building on this, McLaughlin (Reference McLaughlin1991) employed matched asymptotic expansions to extend Saffman’s analysis to higher slip Reynolds numbers, defined as
${\textit{Re}}_s = U_s a / \nu$
, with
$U_s$
the slip velocity, under the regime where
${\textit{Re}}_s$
exceeds
${\textit{Re}}_{\!p}^{1/2}$
while both
${\textit{Re}}_{\!p}$
and
${\textit{Re}}_s$
remain small compared with unity. They demonstrated that, in this limit, the lift decays rapidly and approaches the quiescent flow result. Subsequently, McLaughlin (Reference McLaughlin1993) combined both wall-induced and shear-induced contributions, highlighting their delicate interplay in determining the net lift direction. In bounded shear flows with a single wall, Asmolov (Reference Asmolov1999) analytically showed that inertial lift forces act away from the wall, thereby ruling out near-wall equilibrium at small particle Reynolds number. This prediction was later corroborated by Ekanayake et al. (Reference Ekanayake, Berry, Stickland, Dunstan, Muir, Dower and Harvie2020), who provided numerical fits for lift and drag forces on wall-proximal spheres in shear flow. Complementary finite-element simulations by Feng et al. (Reference Feng, Hu and Joseph1994) for linear shear flow past a neutrally buoyant cylinder showed that the particle remains stably positioned at the centreline for a particle Reynolds number of
${\textit{Re}}_{\!p} = 0.625$
for a moderate confinement ratio,
$\kappa = 0.125$
. Here,
${\textit{Re}}_{\!p} = \kappa ^2 {\textit{Re}}_c$
, with
$\kappa = a/H$
, and
${\textit{Re}}_c = V_w H / \nu$
denotes the channel Reynolds number, where
$H$
is the channel height and
$V_w$
is the wall velocity. Recent lattice Boltzmann (LB) simulations by Fox et al. (Reference Fox, Schneider and Khair2020, Reference Fox, Schneider and Khair2021) revealed that, beyond a critical particle Reynolds number, the centreline equilibrium position becomes unstable, giving rise to symmetric off-centre equilibrium positions. The critical
${\textit{Re}}_{\!p}$
was found to depend strongly on the confinement ratio and channel geometry. Most recently, Anand & Subramanian (Reference Anand and Subramanian2023) analytically investigated the regime of low confinement ratio and low particle Reynolds number, while systematically varying the channel Reynolds number
${\textit{Re}}_c$
. They demonstrated that a supercritical pitchfork bifurcation emerges beyond a critical Reynolds number
${\textit{Re}}_{{cr}}$
, and that lateral migration persists even when
${\textit{Re}}_{\!p} \ll 1$
and
$\kappa \ll 1$
. Furthermore, they identified the channel Reynolds number
${\textit{Re}}_c$
as the appropriate control parameter governing inertial migration in such regimes, particularly for torque-free and neutrally buoyant spheres. Together, these studies establish a comprehensive understanding of inertial migration in Newtonian fluids across a wide range of geometries and flow regimes. Yet, their exclusive focus on Newtonian systems leaves open important questions regarding how non-Newtonian rheology especially shear-thinning behaviour modifies the fundamental mechanisms of lateral migration in confined shear flows at finite Reynolds number. In particular, it remains unclear whether the bifurcation behaviour, equilibrium stability and scaling laws established for Newtonian fluids persist in shear-thinning fluids, or whether new migration regimes emerge due to the nonlinear coupling between inertia and spatially varying viscosity. Addressing these questions is essential for extending inertial microfluidic concepts beyond idealised Newtonian systems toward biological and industrial polymeric fluids.
In this context, the study of particle migration in non-Newtonian fluids has gained increasing attention owing to its importance in biological systems, complex suspensions and industrial processes. Many physiological fluids, including blood, mucus and saliva, exhibit shear-thinning and viscoelastic behaviour (Merrill Reference Merrill1969; Rems et al. Reference Rems, Kawale, Lee and Boukany2016), further motivating detailed investigations of particle transport in such media. Much attention has been devoted to viscoelastic flows, where normal-stress differences and associated elastic forces drive particles toward the centreline in low-inertia regimes in pressure-driven flows (Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007; D’Avino et al. Reference D’Avino, Romeo, Greco and Maffettone2013; Kim & Kim Reference Kim and Kim2016; Raffiee et al. Reference Raffiee, Janmaleki, Saadat, Soleimani and Ebrahimi2017; Del Giudice et al. Reference Del Giudice, D’Avino, Greco, Netti and Maffettone2018; Yu Reference Yu2019).
In contrast, inelastic shear-thinning fluids, which are commonly described using power-law (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot2001) or Carreau rheology (Carreau Reference Carreau1972), have received comparatively less attention in the context of inertial migration than their Newtonian counterparts. Existing investigations of inertial migration in such fluids have been largely restricted to pressure-driven channel flows (Nie & Lin Reference Nie and Lin2015; Li & Xuan Reference Li and Xuan2019; Ouyang Reference Ouyang2019), where shear-rate-dependent viscosity modifies the balance of inertial lift forces and leads to off-centre equilibrium positions (Liu & Hu Reference Liu and Hu2017). In particular, recent numerical studies by Hu and co-workers (Hu et al. Reference Hu, Zhou and Xuan2020, Reference Hu, Zhou and Xuan2023) systematically examined the inertial migration of spherical particles in pressure-driven channel flows of shear-thinning fluids, exploring the effects of the power-law index, Reynolds number, confinement ratio and channel aspect ratio on particle trajectories and equilibrium positions. Their results demonstrated that increasing shear-thinning strength can substantially alter equilibrium locations and focusing behaviour, with particles migrating closer to the channel centreline for larger power-law indices and confinement ratios, and for lower Reynolds numbers. In these pressure-driven configurations, however, the parabolic velocity profile inherently couples shear-gradient-induced lift with wall-induced lift, such that the observed migration behaviour depends simultaneously on flow inertia, geometric confinement and fluid rheology. As a result, the direct contribution of shear-thinning rheology to inertial lift is intrinsically coupled with shear-gradient and geometric effects in such flows. Using linear shear flows in a controlled setting (Rust & Manga Reference Rust and Manga2002) it is possible to isolate and disentangle these contributions. Nevertheless, particle migration in shear-thinning linear shear flows remains comparatively underexplored. The motion of a spherical particle placed at the centre of an unbounded shear flow of shear-thinning fluid under the low Carreau number limit has been studied by Datt & Elfring (Reference Datt and Elfring2018). In the absence of any external force/torque, this analytical study showed that the sphere rotates with the local angular velocity of the fluid, and that there is no quantitative difference compared with a Newtonian fluid. Notable, this study is confined to the Stokes-flow regime and does not address inertial migration or confinement effects. In confined geometries, plane Couette flow offers a particularly attractive configuration, as it generates a uniform shear field in the absence of velocity curvature, thereby eliminating shear-gradient-induced lift and isolating the combined influences of inertia, wall interactions and shear-dependent viscosity. This simplification allows a systematic examination of how shear thinning alters the local balance between inertial and viscous stresses, and consequently the direction and stability of equilibrium positions on particle migration. Despite its experimental accessibility using sliding-plate rheometers (Rust & Manga Reference Rust and Manga2002), inertial migration in planar shear flows of shear-thinning fluids remains largely unexplored. In contrast, the corresponding problem for Newtonian fluids has received considerable attention.
For Newtonian fluids, inertial migration and the bifurcation of equilibrium positions in plane Couette flows have been systematically investigated in recent years (Fox et al. Reference Fox, Schneider and Khair2020, Reference Fox, Schneider and Khair2021; Anand & Subramanian Reference Anand and Subramanian2023), providing a foundation for inertial focusing and separation strategies for a particle in a shear flow. However, shear-thinning fluids, with their shear-rate-dependent viscosity, offer an even richer framework in which inertial migration mechanisms can be amplified or qualitatively altered. An additional motivation for considering linear shear flows is their physical relevance to the near-wall particle dynamics. Particles migrating sufficiently close to a channel or blood vessel wall experience a locally near-linear velocity profile over their own length scale, even in globally parabolic flows, making simple shear a physically relevant approximation in many practical situations. This near-wall linearisation is particularly important in physiological flows such as blood, where deformable red blood cells migrate away from the wall to form a cell-free layer, while stiffer particles such as leukocytes marginate toward the wall and are subjected predominantly to simple shear (Ahmed et al. Reference Ahmed, Baasch, Jang, Pane, Dual and Nelson2017; Alapan et al. Reference Alapan, Bozuyuk, Erkoc, Karacakol and Sitti2020). These observations again underscore the relevance of understanding particle migration in linear shear environments, particularly in complex shear-thinning fluids, where local rheology can strongly influence the migration dynamics. In this study, we address these gaps by systematically investigating the inertial migration of a neutrally buoyant spherical particle suspended in a shear-thinning Carreau fluid (Carreau Reference Carreau1972) under plane Couette flow. Using smoothed profile–LB simulations, we examine how variations in key dimensionless parameters governing shear-thinning rheology influence lateral migration, equilibrium positions and stability. Our results demonstrate that, even in the absence of shear-gradient-induced lift, a shear-thinning rheology can destabilise the centreline equilibrium and reduce the critical Reynolds number for off-centred migration, highlighting plane Couette flow as a simplified yet insightful configuration complementary to pressure-driven studies and offering guidance for the design of inertial microfluidic separation devices.
2. Problem statement and methodology
2.1. Problem description
We consider the motion of a single neutrally buoyant spherical particle of radius
$ a$
in plane Couette flow of an incompressible, inelastic shear-thinning fluid. The fluid rheology is described by the Carreau model (Carreau Reference Carreau1972) with constant density
$ \rho$
. The flow is confined between two infinite parallel plane walls located at
$ z = 0$
and
$ z = H$
. The top wall moves with a constant velocity
$ V_w \boldsymbol{e}_x$
, while the bottom wall is stationary, resulting in a linear shear profile in the absence of the particle, characterised by a shear rate
$ G = V_w/H$
. The computational domain is a rectangular parallelepiped of size
$ L \times W \times H$
, where
$L$
,
$W$
and
$H$
denote the streamwise, spanwise and wall-normal directions, respectively, as shown in figure 1. We non-dimensionalise the problem using the channel height
$ H$
as the length scale, top wall velocity
$V_{w}$
as the velocity scale and inverse shear rate
$ G^{-1}$
as the time scale. The particle is initially located at
$ Z_0$
and is free to translate and rotate under hydrodynamic forces and torques. The flow field and particle motion are governed by the channel Reynolds number, defined as
$ {\textit{Re}}_c = V_w H / \nu _0$
, where
$ \nu _0$
is the zero-shear kinematic viscosity of the fluid. The geometric confinement is quantified by the confinement ratio
$ \kappa = a/H$
. The particle Reynolds number, defined as
$ {\textit{Re}}_{\!p} = {\textit{Re}}_c \kappa ^2$
, varies with
$ {\textit{Re}}_c$
, allowing for the investigation of inertial effects beyond the creeping-flow regime. In the present study, an in-house solver based on the smoothed profile–lattice Boltzmann method (SP–LBM) is employed to simulate the coupled dynamics of the particle and the fluid (Jafari, Yamamoto & Rahnama Reference Jafari, Yamamoto and Rahnama2011). The particle’s translation and rotation are governed by Newton’s equations of motion, and the fluid–particle interactions are handled through the smoothed profile method (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005). The local shear-rate-dependent viscosity is updated dynamically according to the Carreau model to capture shear-thinning effects.
Schematic of the computational domain for particle in a confined shear flow. A neutrally buoyant spherical particle of radius
$a$
is initially placed at a distance
$z_0$
from the bottom wall between two parallel plates in a three-dimensional domain of length
$L$
(streamwise,
$x$
-direction), width
$W$
(spanwise,
$y$
-direction) and height
$H$
(wall-normal,
$z$
-direction). The bottom plate is stationary, while the top plate moves in the
$x$
-direction with velocity
$V_w$
, thereby generating a shear flow.
2.2. Carreau fluid rheology model
A Carreau fluid is a subclass of inelastic non-Newtonian fluid that captures shear-thinning behaviour (Carreau Reference Carreau1972). This constitutive model is commonly used to represent biological fluids (e.g. blood, mucus) and polymeric liquids (e.g. Xanthan gum), providing a realistic description of shear-thinning effects. The deviatoric part of the stress tensor is expressed as
$\sigma _{\alpha \beta } = 2 \mu (\dot {\gamma }) S_{\alpha \beta }$
, where
$S_{\alpha \beta }$
is the rate-of-strain tensor. The shear-rate-dependent viscosity is defined by the Carreau constitutive relation
\begin{equation} \mu (\dot {\gamma }) = \mu _{\infty } + (\mu _{0} - \mu _{\infty }) \left [ 1 + (\lambda |\dot {\gamma }|)^2 \right ]^{\frac {n-1}{2}}, \end{equation}
where
$\mu _0$
and
$\mu _\infty$
denote the dynamic viscosities at zero- and infinite-shear rates, respectively,
$\lambda$
is the time constant,
$n$
is the power-law index and
$|\dot {\gamma }|$
is the magnitude of the shear rate. The rheological behaviour of the Carreau fluid is further characterised by three key non-dimensional parameters: (i) the Carreau number
$Cu = \lambda G$
, which compares the fluid relaxation time with the characteristic flow time scale; (ii) the viscosity ratio
$\beta = \mu _{\infty }/\mu _0$
, which quantifies the degree of shear thinning; and (iii) the power-law index
$n$
, which governs the sharpness of the viscosity transition with shear rate.
2.3. Fluid and flow parameters
The parameter ranges explored in this study are motivated by both microfluidic and physiological applications. Typical channel Reynolds numbers in inertial microfluidic devices and biological flows span
${\textit{Re}}_c = \mathcal{O}(10^{-2}-10^2)$
, while shear rates encountered in blood and polymeric solutions correspond to Carreau numbers in the range
$\mathcal{O}(10^{-2}$
–
$10^1)$
. Using representative physiological values for blood density
$\rho \approx 1050\,\mathrm{kg\,m^{-3}}$
and dynamic viscosity
$\mu \approx 4.5 \times 10^{-3}\,\mathrm{Pa\,s}$
, Reynolds numbers in veins and arteries with characteristic diameters of a few millimetres and mean velocities of order
$10^{-2}$
–
$10^{-1} \,\mathrm{m\,s^{-1}}$
(Bozuyuk, Ozturk & Sitti Reference Bozuyuk, Ozturk and Sitti2023) fall in the range
${\textit{Re}}_{c} \sim \mathcal{O}(10)$
–
$\mathcal{O}(10^2)$
. In smaller vessels such as arterioles and venules, with diameters of tens of micrometres and velocities of order
$10^{-3}$
–
$10^{-2}\,\mathrm{m\,s^{-1}}$
, the Reynolds number decreases to
${\textit{Re}}_{c} \sim \mathcal{O}(10^{-2})$
–
$\mathcal{O}(1)$
, while capillary flows are essentially inertia free. Similar Reynolds number regimes have been reported in recent computational studies of microscale transport and microrobotic locomotion in blood vessels, depending on vessel size and flow conditions (Bozuyuk et al. Reference Bozuyuk, Ozturk and Sitti2023). Further, the characteristic shear-thinning time scale
$\lambda$
of biological and polymeric fluids spans a wide range depending on fluid composition and microstructure. For blood, reported values of
$\lambda$
range from fractions of a second to a few seconds, reflecting red blood cell deformation and aggregation dynamics (Cho & Kensey Reference Cho and Kensey1991; Miladi Reference Miladi2011), whereas larger relaxation times have been observed in other biological fluids such as mucus (Li et al. Reference Li2008) and in dilute aqueous polymer solutions, including Xanthan gum, which are widely used as experimental model fluids for shear-thinning rheology (Boyko & Stone Reference Boyko and Stone2021).
Representative values of shear-thinning parameters for few physiological and polymeric fluids, including zero-shear and infinite-shear viscosities, viscosity ratio, power-law index and relaxation time, are summarised in table 1. These values provide a practical guide for the parameter ranges considered in the present simulations. However, to systematically assess the influence of shear-thinning rheology, we conduct a broad parametric study spanning a wide range of Carreau numbers, viscosity ratios and power-law indices. Specifically, the Carreau number span
$Cu = 0$
–
$100$
, viscosity ratio
$\beta = 0.1$
–
$1$
and the power-law index
$n = 0.25$
–
$1$
, allowing exploration of fluids from weakly to strongly shear thinning. For example, blood exhibits high-shear to zero-shear viscosity ratios under physiological shear conditions (Cho & Kensey Reference Cho and Kensey1991), while polymer solutions and melts display comparable reductions in effective viscosity at elevated shear rates (Larson Reference Larson1999). Similarly, the Xanthan gum solution studied by Boyko & Stone (Reference Boyko and Stone2021) shows pronounced shear-thinning behaviour with
$\mu _0 = 11.9\,\mathrm{Pa\,s}$
,
$\mu _\infty = 2.31 \times 10^{-3}\,\mathrm{Pa\,s}$
,
$\beta = 1.94 \times 10^{-4}$
,
$n = 0.279$
and
$\lambda = 30.4\,\mathrm{s}$
under microchannel flow conditions. These examples motivate us to simulate systems over a wide range of Reynolds number, Carreau number, viscosity ratio and power-law index relevant to realistic inertial microfluidic and physiological flows.
Representative shear-thinning fluids and typical Carreau model parameters reported in the literature.
2.4. Lattice Boltzmann method
To solve the fluid flow, we employ the LBM coupled with the smoothed profile method (SPM) to simulate the interaction between a rigid spherical particle and a Carreau fluid under shear flow. The LBM is a lattice-based numerical approach for simulating fluid flows governed by the Navier–Stokes equations (Succi Reference Succi2001). In LBM, virtual packets of mesoscopic particles propagate and collide on a discrete lattice (Aidun & Clausen Reference Aidun and Clausen2010). The fluid domain is discretised into a series of Eulerian nodes at positions
$\boldsymbol{x} = x \boldsymbol{e}_x + y \boldsymbol{e}_y + z \boldsymbol{e}_z$
and times
$t$
. The evolution of the distribution function
$f_i$
is governed by the discretised Boltzmann equation with the Bhatnagar–Gross–Krook (BGK) collision operator (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017)
\begin{equation} f_{i}\left ( \boldsymbol{x} + \boldsymbol{c}_{i} \Delta t, t + \Delta t \right ) - f_{i} \left ( \boldsymbol{x}, t \right ) = -\frac {\Delta t}{\tau } \left [ f_{i}\left ( \boldsymbol{x}, t \right ) - f_{i}^{\textit{eq}} \left ( \boldsymbol{x}, t \right ) \right ] + F_{i}(\boldsymbol{x},t) \Delta t, \end{equation}
where
$f_i(\boldsymbol{x}, t)$
is the particle distribution function along the
$i$
th lattice velocity
$\boldsymbol{c}_i$
,
$\tau$
is the relaxation time and
$\Delta t$
is the time step. Here,
$F_i(\boldsymbol{x},t)$
is the forcing term, which incorporates external forces using Guo’s forcing scheme (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002)
\begin{equation} F_i(\boldsymbol{x},t) = w_i \left [ \frac {(\boldsymbol{c}_i - \boldsymbol{u})}{c_s^2} + \frac {(\boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{u}) \boldsymbol{c}_i}{c_s^4} \right ] \boldsymbol{\cdot }\boldsymbol{F}, \end{equation}
where
$w_i$
are the lattice weights,
$\boldsymbol{F}$
is the external force,
$\boldsymbol{u}$
is the macroscopic velocity and
$c_s = ({1}/{\sqrt {3}}) ({\Delta x}/{\Delta t})$
is the lattice speed of sound. The equilibrium distribution
$f_i^{\textit{eq}}$
is given by (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017)
\begin{equation} f_i^{\textit{eq}} = \rho \,w_i \left [ 1 + \frac {\boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{u}}{c_s^2} + \frac {(\boldsymbol{c}_i \boldsymbol{\cdot }\boldsymbol{u})^2}{2 c_s^4} - \frac {\boldsymbol{u}^2}{2 c_s^2} \right ]. \end{equation}
The evolution proceeds in two steps: streaming, where
$f_i$
propagates along
$\boldsymbol{c}_i$
to neighbouring nodes, and collision, where
$f_i$
relaxes toward
$f_i^{\textit{eq}}$
, typically given by a second-order expansion of the Maxwell–Boltzmann distribution following (2.4). The lattice is discretised using the three-dimensional, nineteen-velocity lattice model (D3Q19), with discrete velocities
$\boldsymbol{c}_i$
and weights
$w_i$
summarised in table 2. The relaxation time
$\tau$
is related to the fluid kinematic viscosity
$\nu$
via
$\nu = c_s^2 (\tau - ({1}/{2})) \Delta t$
. To model the Carreau fluid, the local viscosity
$\mu (\dot {\gamma })$
is updated based on the local shear-rate magnitude
$|\dot {\gamma }| = \sqrt {2 \sum _{\alpha ,\beta } S_{\alpha \beta } S_{\alpha \beta }}$
, which is computed from the rate-of-strain tensor
$S_{\alpha \beta } = -( {1}/({2 \tau c_s^2})) \sum _i c_{i\alpha } c_{i\beta } f_i^{\textit{neq}}$
, where
$f_i^{\textit{neq}}$
are the non-equilibrium distributions (Wang & Ho Reference Wang and Ho2011). Therefore, the relaxation time is updated as
$\tau (\dot {\gamma }) = ({\mu (\dot {\gamma })}/({\rho c_s^2 \Delta t})) + ({1}/{2})$
, ensuring self-consistent adaptation to the local shear rate. It is important thing to note that, since the shear rate
$\dot {\gamma }$
depends on the relaxation time
$\tau$
through the non-equilibrium distribution functions, the above relation for
$\tau (\dot {\gamma })$
leads to a nonlinear implicit equation for
$\tau$
in terms of
$\dot {\gamma }$
. This equation may be solved iteratively at each lattice node (Phillips & Roberts Reference Phillips and Roberts2011), but such an approach can be computationally expensive, particularly for large-scale simulations. As an efficient alternative, the relaxation time from the previous time step is used while computing
$\dot {\gamma }$
. This approximation has been shown to provide stable and accurate results for non-Newtonian LB simulations (Phillips & Roberts Reference Phillips and Roberts2011). In the present simulations, this explicit formulation is employed to reduce computational cost while maintaining numerical stability and accuracy. Macroscopic quantities are recovered as
$\rho (\boldsymbol{x},t) = \sum _i f_i(\boldsymbol{x},t)$
and
$\rho \boldsymbol{u}(\boldsymbol{x},t) = \sum _i f_i(\boldsymbol{x},t) \boldsymbol{c}_i + ( {1}/{2})\boldsymbol{F} \Delta t$
. A detailed description of the numerical implementation, including the complete SP–LBM algorithm, the viscosity update procedure and the individual computational steps, is provided in Appendix A. At this point, it is worth emphasising that two distinct relaxation times are employed in the present formulation. The reference relaxation time
$\tau _{0}$
is defined based on the zero-shear viscosity
$\mu _{0}$
, whereas the relaxation time
$\tau (\dot {\gamma })$
is evaluated locally from the shear-rate-dependent apparent viscosity
$\mu (\dot {\gamma })$
prescribed by the Carreau model. Consequently,
$\tau$
varies spatially according to the local shear rate
$\dot {\gamma }$
, reflecting the underlying shear-thinning rheology of the fluid.
Discrete velocities
$\boldsymbol{c}_i$
and corresponding weights
$w_i$
for the D3Q19 lattice.
2.5. Particle–fluid coupling using the smoothed profile method
To couple particles with the host fluid, we adopt the SPM (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Yamamoto, Molina & Nakayama Reference Yamamoto, Molina and Nakayama2021). In this approach, the fluid–particle interface is represented as a diffuse interface of finite thickness, rather than a sharp boundary, which allows for smooth coupling between the Eulerian fluid fields and the Lagrangian particle dynamics. Within this framework, the no-slip and no-penetration conditions at the particle surface are not imposed directly; instead, they are satisfied through the introduction of a body force to the momentum equation, that mediates the interaction between the particle and the surrounding fluid. Following Jafari et al. (Reference Jafari, Yamamoto and Rahnama2011), Mino et al. (Reference Mino2017) and Sobhani, Bazargan & Sadeghy (Reference Sobhani, Bazargan and Sadeghy2019), the LBM is used to solve for the flow field on fixed Eulerian grid, while particle motion is handled using the SPM. The profile function
$\phi$
is defined as
\begin{equation} \phi (\boldsymbol{x}) = \frac {1}{2} \left [ 1 + \tanh \left ( \frac {a - \|\boldsymbol{x} - \boldsymbol{X}\|}{\xi } \right ) \right ]\!, \end{equation}
where
$\boldsymbol{X} = X \boldsymbol{e}_x + Y \boldsymbol{e}_y + Z \boldsymbol{e}_z$
is the centroid position of the particle,
$\boldsymbol{x}$
is the Eulerian fluid node,
$a$
is the particle radius and
$\xi$
is the interfacial thickness. Unless otherwise stated, we set
$\xi = 2$
lattice units following Nakayama & Yamamoto (Reference Nakayama and Yamamoto2005) and Jafari et al. (Reference Jafari, Yamamoto and Rahnama2011), which ensures sufficient smoothness of the interface while maintaining accuracy. The interaction force is given by
\begin{equation} \boldsymbol{f}_{\!p} = \phi \frac {\boldsymbol{u} - \boldsymbol{u}_{\!p}}{\Delta t}, \quad \text{with } \boldsymbol{u}_{\!p} = \boldsymbol{U} + \boldsymbol{\varOmega } \times (\boldsymbol{x} - \boldsymbol{X}), \end{equation}
where
$\boldsymbol{U}$
and
$\boldsymbol{\varOmega }$
represent translational and angular velocities of the particle, respectively. The hydrodynamic force is
$\boldsymbol{f}_h = -\boldsymbol{f}_{\!p}$
, which is applied in the LB collision step using Guo’s forcing scheme (Guo et al. Reference Guo, Zheng and Shi2002). The total hydrodynamic force and torque on the particle are computed as follows:
\begin{equation} \boldsymbol{F}_{\!H} = \sum _{\boldsymbol{x} \in V_{\!p}} \rho (\boldsymbol{x}) \boldsymbol{f}_{\!p}, \quad \boldsymbol{T}_{\!H} = \sum _{\boldsymbol{x} \in V_{\!p}} (\boldsymbol{x} - \boldsymbol{X}) \times \rho (\boldsymbol{x}) \boldsymbol{f}_{\!p}, \end{equation}
where
$V_{\!p}$
denotes the particle domain. Particle velocities are updated using the improved SPM (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Mino et al. Reference Mino2017)
\begin{equation} \boldsymbol{U}(t+\Delta t) = \boldsymbol{U}(t) + \frac {\boldsymbol{F}_{\!H} + \boldsymbol{F}_{\!H,\textit{in}}}{M_{\!p}} \Delta t, \quad \boldsymbol{\varOmega }(t+\Delta t) = \boldsymbol{\varOmega }(t) + \frac {\boldsymbol{T}_{\!H} + \boldsymbol{T}_{\!H,\textit{in}}}{I_{\!p}} \Delta t, \end{equation}
where
$M_{\!p}$
and
$I_{\!p}$
are the particle mass and moment of inertia, respectively, and
$\boldsymbol{F}_{\!H,\textit{in}}$
and
$\boldsymbol{T}_{\!H,\textit{in}}$
are compensation terms accounting for internal mass effects (Mino et al. Reference Mino2017). Finally, the particle position is updated using the Crank–Nicolson scheme (Sobhani et al. Reference Sobhani, Bazargan and Sadeghy2019)
\begin{equation} \boldsymbol{X}(t+\Delta t) = \boldsymbol{X}(t) + \frac {\boldsymbol{U}(t+\Delta t) + \boldsymbol{U}(t)}{2} \Delta t. \end{equation}
2.6. Numerical set-up and robustness
The numerical parameters are chosen to ensure accuracy and stability while adequately resolving both particle motion and flow field. All simulations are performed with a fixed particle radius of
$a = 30$
lattice units. The relaxation time at zero shear rate,
$\tau _0 = 0.8$
, corresponds to the fluid viscosity at zero shear,
$\mu _0$
. For shear-thinning fluids, the local relaxation time
$\tau (\dot {\gamma })$
is updated dynamically based on the instantaneous shear rate to capture the shear-dependent viscosity, while
$\tau _0$
serves as the reference value for the Newtonian limit at vanishing shear. The domain aspect ratio is set to
$ AR = L/H = W/H = 2$
, ensuring that periodicity in the streamwise and spanwise directions does not influence the particle dynamics. The computational domain is discretised using a uniform three-dimensional lattice of size
$ N_x \times N_y \times N_z$
, where
$N_x$
,
$N_y$
and
$N_z$
denote the number of lattice nodes in the
$x$
,
$y$
and
$z$
directions, respectively. Periodic boundary conditions are applied in the flow and vorticity directions, while no-slip and no-penetration boundary conditions are imposed at the top and bottom walls using the mid-grid bounce-back scheme (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017). To achieve different channel Reynolds numbers, the wall velocity
$V_w$
is adjusted accordingly while maintaining a low lattice Mach number,
$Ma = V_w / c_s \lt 0.1$
, to minimise compressibility effects. All simulations are carried out in lattice units with
$\Delta x = \Delta t = 1$
. Diffusive scaling is enforced through the mapping to physical units, such that the physical lattice spacing and time step satisfy
$\Delta x_{\textit{ph}} \propto Ma$
and
$\Delta t_{\textit{ph}} \propto Ma^2$
, ensuring convergence to the incompressible Navier–Stokes equations (Krüger et al. Reference Krüger, Varnik and Raabe2009; Luo, Lallemand & d’Humières Reference Luo, Lallemand and d’Humières2011). The grid, time-step and domain-independence tests are conducted for an extreme parameter set considered (
${\textit{Re}}_c = 300$
,
$Cu = 100$
,
$\beta = 0.1$
,
$n = 0.25$
and
$\kappa = 0.1$
). These tests confirm that the particle migration trajectories and equilibrium positions remain unchanged with increased lattice resolution, reduced effective time step or larger domain aspect ratio. Detailed numerical verification, including Mach number constraints, grid resolution, time-step dependence and domain-size independence, is provided in Appendix B. Overall, these tests demonstrate that the simulations reported here are numerically robust and insensitive to further refinement of the lattice, time step or domain size.
2.7. Validation of the numerical implementation
To validate the accuracy of our numerical implementation, we consider the inertial migration of a neutrally buoyant spherical particle in plane Couette flow of a Newtonian fluid. Simulations are performed over a range of channel Reynolds numbers
$({\textit{Re}}_c)$
while maintaining a small confinement ratio
$\kappa = 0.05$
, such that both the particle Reynolds number
$({\textit{Re}}_{\!p})$
and the confinement ratio
$(\kappa )$
remain sufficiently small. This ensures negligible particle-scale inertia, consistent with the assumptions underlying theoretical predictions of Anand & Subramanian (Reference Anand and Subramanian2023). Figure 2 shows the non-dimensional equilibrium position of the particle as a function of
${\textit{Re}}_c$
. It can be seen that, up to a certain channel Reynolds number, the particle remains at the channel centreline at steady state, whereas beyond a certain channel Reynolds number defined as the critical Reynolds number,
${\textit{Re}}_{{cr}}$
, it migrates to off-centre equilibrium positions. The equilibrium positions obtained from our simulations are compared with the supercritical pitchfork bifurcation behaviour reported by Anand & Subramanian (Reference Anand and Subramanian2023), showing excellent agreement and demonstrating the accuracy and robustness of the implemented SP–LBM framework in capturing particle migration phenomena in confined Newtonian flows. However, at higher
${\textit{Re}}_c$
(e.g.
${\textit{Re}}_c \gt 200$
), a slight deviation from the theoretical predictions is observed. This arises because the theory assumes both
$\kappa \ll 1$
and
${\textit{Re}}_{\!p} \ll 1$
, whereas in our simulations, although
$\kappa = 0.05$
and
${\textit{Re}}_{\!p} \lt 1$
, these assumptions become only marginally valid at larger
${\textit{Re}}_c$
. To more closely capture the theoretical predictions, an even smaller confinement ratio would be required, which in turn necessitates a larger computational domain and significantly higher computational cost. Further, the inset in figure 2 illustrates the square-root scaling near the bifurcation threshold, where
$\epsilon = ({\textit{Re}}_c - {\textit{Re}}_{{cr}})/{\textit{Re}}_{{cr}}$
. The scaling behaviour from our simulations agrees well with the theoretical predictions of Anand & Subramanian (Reference Anand and Subramanian2023), further validating the numerical framework.
Validation of the simulated bifurcation behaviour of a neutrally buoyant sphere in confined shear flow against the theoretical predictions of Anand & Subramanian (Reference Anand and Subramanian2023). The plot shows the normalised equilibrium position of the particle,
$Z_{\textit{eq}}/H$
, as a function of the channel Reynolds number,
${\textit{Re}}_c$
. The results conform to a supercritical pitchfork bifurcation, with the particle migrating from the channel centreline to off-centre equilibrium positions as
${\textit{Re}}_c$
exceeds the critical threshold
${\textit{Re}}_{{cr}}$
. The inset illustrates the square-root scaling of the bifurcated branches in the vicinity of
${\textit{Re}}_{{cr}}$
, where
$\epsilon = ({\textit{Re}}_c - {\textit{Re}}_{{cr}})/{\textit{Re}}_{{cr}}$
, in agreement with the asymptotic theory of Anand & Subramanian (Reference Anand and Subramanian2023). The parameters used in the simulations are confinement ratio
$\kappa = 0.05$
and particle Reynolds number
${\textit{Re}}_{\!p} \lt 1$
. The inset validates the scaling plot (
$Z_{\textit{eq}}/H-0.5$
versus
$\epsilon$
) reported in figure 5 of Anand & Subramanian (Reference Anand and Subramanian2023).
3. Results and discussion
3.1. Effect of Carreau number, viscosity ratio and power-law index
We begin by examining the influence of shear thinning on the lateral migration of a neutrally buoyant spherical particle in plane Couette flow, with a fixed confinement ratio of
$\kappa = 0.2$
(which is maintained throughout this section unless otherwise specified). The extent of shear thinning is characterised by the Carreau number (
$Cu$
), which represents the relative importance of the shear rate to the intrinsic time scale of the fluid; larger values of
$Cu$
correspond to stronger shear-thinning effects. The viscosity ratio (
$\beta$
) quantifies the contrast between the infinite-shear and zero-shear viscosities, thereby controlling the strength of viscosity variation between regions of high and low shear rate. The power-law index (
$n$
) governs the degree of shear thinning in the intermediate shear-rate regime, with smaller values of
$n$
corresponding to more pronounced shear-thinning behaviour. We first examine the effect of the Carreau number by analysing the bifurcation behaviour of the equilibrium particle position. Figure 3(a) presents the regime diagram in the
$({\textit{Re}}_c, Cu)$
parameter space for
$\beta = 0.1$
and
$n = 0.25$
. The solid demarcating line represents the critical condition separating two distinct migration regimes. For a fixed Carreau number, this boundary corresponds to a critical Reynolds number
${\textit{Re}}_{{cr}}$
; conversely, for a fixed Reynolds number, it defines a critical Carreau number
$Cu_{{cr}}$
. Below this boundary (i.e.
${\textit{Re}}_c \lt {\textit{Re}}_{{cr}}$
for fixed
$Cu$
, or
$Cu \lt Cu_{{cr}}$
for fixed
${\textit{Re}}_c$
), the particle ultimately migrates to the channel centreline irrespective of its initial release position. In contrast, beyond the critical condition (
${\textit{Re}}_c \gt {\textit{Re}}_{{cr}}$
or
$Cu \gt Cu_{{cr}}$
), the centreline becomes unstable and the particle migrates towards an off-centre equilibrium position. The final equilibrium location depends on the initial release position. If the particle is released between the centreline and the lower half of the channel, it migrates to a stable off-centre position within that region. Similarly, if released between the centreline and the upper wall, it approaches a symmetric off-centre equilibrium in the upper half of the channel. Thus, the bifurcation is supercritical in nature: once the critical threshold (
${\textit{Re}}_{{cr}}$
for fixed
$Cu$
or
$Cu_{{cr}}$
for fixed
${\textit{Re}}_c$
) is exceeded, the centreline loses stability and symmetric off-centre equilibria emerge. Further insight into the nature of the bifurcation is obtained by examining the dependence of the equilibrium position on the Carreau number. Figure 3(b) shows the equilibrium position
$Z_{{\textit{eq}}}/H$
as a function of
$Cu$
for a fixed channel Reynolds number
${\textit{Re}}_c = 100$
, with
$\beta = 0.1$
and
$n = 0.25$
. This analysis is aimed at characterising the scaling behaviour of particle migration in a shear-thinning fluid near the onset of instability. For sufficiently small
$Cu$
, the steady equilibrium remains at the channel centreline (
$Z_{{\textit{eq}}}/H = 0.5$
), indicating stability of the symmetric state. Beyond a threshold value of
$Cu$
, the centreline loses stability and the particle migrates to an off-centre equilibrium position (
$Z_{{\textit{eq}}}/H \neq 0.5$
), with the direction determined by the initial release location. Furthermore, increasing
$Cu$
drives the particle progressively farther from the centreline and closer to the wall, reflecting the strengthening influence of shear-thinning effects. The value of
$Cu$
at which the transition from centreline to off-centre equilibrium occurs is defined as the critical Carreau number,
$Cu_{{cr}}$
. For the present parameter set, the transition is observed between
$Cu = 1$
and
$Cu = 2$
, yielding an estimate
$Cu_{{cr}} \approx 1.2$
. This estimate is obtained from simulations performed with refined parameter sampling in this interval (using increments of Carreau number
$\Delta Cu = 0.1$
) near the bifurcation point. In general,
$Cu_{{cr}}$
depends on
${\textit{Re}}_c$
,
$n$
,
$\beta$
and
$\kappa$
; for example, figure 3(a) shows that
$Cu_{{cr}}$
varies with the channel Reynolds number
${\textit{Re}}_c$
. In the present discussion, however, we restrict attention to this representative case. For
$Cu \lt Cu_{{cr}}$
, the particle migrates to the centreline irrespective of its initial position, whereas for
$Cu \gt Cu_{{cr}}$
, symmetric off-centre equilibria emerge and the final state depends on the initial condition. The inset of figure 3(b) highlights the near-critical behaviour using the reduced control parameter
$\delta = (Cu - Cu_{{cr}})/Cu_{{cr}}$
, with the vertical axis plotted as
$Z_{{\textit{eq}}}/H - 0.5$
. The observed square-root variation of
$\delta$
with
$Z_{{\textit{eq}}}/H - 0.5$
near onset is consistent with the canonical scaling of a supercritical pitchfork bifurcation. Additional insight is obtained by examining the variation of the equilibrium position with the channel Reynolds number. Figure 3(c) shows the steady-state equilibrium position
$Z_{{\textit{eq}}}/H$
as a function of
${\textit{Re}}_c$
for several values of
$Cu$
, with
$n = 0.25$
and
$\beta = 0.1$
. For each fixed
$Cu$
, the particle remains at the channel centreline (
$Z_{{\textit{eq}}}/H = 0.5$
) up to a critical Reynolds number, beyond which it migrates to an off-centre equilibrium position, thereby forming a bifurcation diagram in
${\textit{Re}}_c$
. The Reynolds number at which this transition occurs is defined as the critical Reynolds number,
${\textit{Re}}_{{cr}}$
. For example, in the Newtonian limit (
$Cu = 0$
) and for the present confinement ratio
$\kappa = 0.2$
, the critical value is approximately
${\textit{Re}}_{{cr}} \approx 153$
. For
${\textit{Re}}_c \gt {\textit{Re}}_{{cr}}$
, the particle migrates away from the centreline and settles at an off-centre equilibrium position (
$Z_{{\textit{eq}}}/H \neq 0.5$
). As
$Cu$
increases, this critical Reynolds number
${\textit{Re}}_{{cr}}$
decreases monotonically, indicating that shear-thinning effects reduce the inertial threshold required to destabilise the centreline equilibrium. Consistent with the trends observed earlier, for a fixed
${\textit{Re}}_c$
, larger values of
$Cu$
also shift the off-centre equilibrium progressively closer to the channel wall.
(a) Regime diagram in the
$({\textit{Re}}_c,\,Cu)$
parameter space for
$\kappa = 0.2$
,
$\beta = 0.1$
and
$n = 0.25$
. The solid line represents the critical Reynolds number,
${\textit{Re}}_{{cr}}$
, separating two distinct migration regimes. For
${\textit{Re}}_c$
values to the right of the boundary, the particle migrates to an off-centre equilibrium position, whereas for
${\textit{Re}}_c$
values to the left of the boundary, the particle returns to the channel centreline. (b) Normalised equilibrium position of the particle,
$Z_{{\textit{eq}}}/H$
, as a function of the Carreau number,
$Cu$
. The results exhibit a supercritical pitchfork bifurcation, in which the particle migrates from the channel centreline to off-centre equilibrium positions as
$Cu$
exceeds the critical value,
$Cu_{{cr}}$
. Note that this
$Cu_{{cr}}$
is a function of
${\textit{Re}}_c$
. Here, the bifurcation shown corresponds to
${\textit{Re}}_c = 100$
,
$\kappa = 0.2$
,
$\beta = 0.1$
and
$n = 0.25$
. Inset: variation of
$Z_{{\textit{eq}}}/H - 0.5$
with
$\delta = (Cu - Cu_{{cr}})/Cu_{{cr}}$
near the bifurcation point, showing square-root scaling close to onset. (c) Normalised equilibrium position across the channel as a function of the channel Reynolds number
${\textit{Re}}_c$
for different Carreau numbers, with
$\beta = 0.1$
,
$\kappa = 0.2$
and
$n = 0.25$
.
(a) Variation of the critical Reynolds number
${\textit{Re}}_{{cr}}$
with the Carreau number
$Cu$
for different viscosity ratios
$\beta$
, with
$\kappa = 0.2$
and
$n = 0.25$
. (b) Dependence of the critical Reynolds number
${\textit{Re}}_{{cr}}$
on the Carreau number
$Cu$
for different power-law indices
$n$
, with
$\beta = 0.5$
and
$\kappa = 0.2$
.
Having established the influence of the Carreau number on the critical Reynolds number, we now examine how the remaining rheological parameters such as viscosity ratio
$\beta$
and power-law index
$n$
affect the onset of lateral migration. The dependence of
${\textit{Re}}_{{cr}}$
on
$Cu$
for different viscosity ratios is shown in figure 4(a) for
$\beta = 0.9$
,
$0.5$
and
$0.1$
, with
$n = 0.25$
. In each case,
${\textit{Re}}_{{cr}}$
is determined from simulations performed with unit increments of
${\textit{Re}}_c$
near the transition point. For all viscosity ratios,
${\textit{Re}}_{{cr}}$
decreases monotonically with increasing
$Cu$
. Moreover, the reduction becomes more pronounced as
$\beta$
decreases, corresponding to a larger viscosity contrast between low- and high-shear regions and therefore stronger shear-thinning effects. We next consider the role of the power-law index
$n$
. Figure 4(b) shows the variation of
${\textit{Re}}_{{cr}}$
with
$Cu$
for
$n = 0.25$
,
$0.5$
and
$0.75$
, at a fixed viscosity ratio
$\beta = 0.5$
. The results indicate that
${\textit{Re}}_{{cr}}$
decreases as
$n$
decreases, with the sensitivity becoming increasingly significant at lower values of
$n$
. Since smaller
$n$
corresponds to stronger shear-thinning behaviour, this trend is consistent with the variations observed with
$Cu$
and
$\beta$
. Taken together, these results demonstrate that enhanced shear thinning, whether achieved by increasing
$Cu$
, decreasing
$\beta$
or reducing
$n$
, lowers the critical Reynolds number required to destabilise the centreline equilibrium and promotes particle migration at weaker inertial forcing.
Dimensionless lift force at steady state on a neutrally buoyant sphere in confined Couette flow. (a) Variation of the non-dimensional lift force with transverse position for different Carreau numbers, with
$\kappa = 0.2$
,
$n = 0.25$
,
$\beta = 0.1$
and
${\textit{Re}}_c = 150$
. (b) Effect of viscosity ratio
$\beta$
on the non-dimensional lift force at
${\textit{Re}}_c = 150$
,
$n = 0.25$
,
$Cu = 10$
and
$\kappa = 0.2$
. (c) Effect of the power-law index
$n$
on the non-dimensional lift force for
$\kappa = 0.2$
,
${\textit{Re}}_c = 150$
,
$\beta = 0.5$
and
$Cu = 10$
. Insets in all panels show enlarged views near the zero crossing of the normalised lift force.
To further understand the bifurcation behaviour described above, we examine the stability of both the centreline and off-centre equilibrium positions by analysing the lateral force acting on the particle. Specifically, we compute the non-dimensional lift force as a function of the particle’s lateral position across the channel height, as shown in figure 5(a–c). The lift force is evaluated by placing the particle at various
$z$
-locations while allowing free translation and rotation in all directions except
$z$
, where translation is constrained. The resulting lift force is normalised using
$F_s = \rho V_w^2 a^4 / H^2$
(Ho & Leal Reference Ho and Leal1974), allowing direct comparison across different flow conditions and rheological parameters. Figure 5 presents the lift-force profiles for variations in the Carreau number (
$Cu$
), viscosity ratio (
$\beta$
) and power-law index (
$n$
), with all other parameters held constant. Figure 5(a) illustrates the effect of
$Cu$
for
$n = 0.25$
,
$\beta = 0.1$
and
${\textit{Re}}_c = 150$
. For a Newtonian fluid (
$Cu = 0$
), the lift-force profile exhibits a single zero crossing at the channel centreline, with a negative slope, indicating that the centreline is a stable equilibrium. As
$Cu$
increases, reflecting stronger shear-thinning effects, two additional off-centre zero crossings emerge. These off-centre crossings have negative slopes, identifying them as stable equilibrium positions, while the slope at the centreline becomes positive, signalling the loss of its stability. This transition indicates that stronger shear thinning promotes off-centre migration, even when the particle initially resides at the centreline. Figure 5(b) explores the role of the viscosity ratio
$\beta$
for
$Cu = 10$
,
$n = 0.25$
and
${\textit{Re}}_c = 150$
. For higher
$\beta$
values (closer to unity), the lift-force profiles resemble those of Newtonian fluids, with a single stable equilibrium at the centreline. Reducing
$\beta$
, which increases the contrast between zero and infinite shear viscosities, destabilises the centreline and promotes the appearance of off-centre stable equilibria. This demonstrates that the distribution of shear-dependent viscosity across the channel can significantly influence particle migration, with stronger viscosity contrasts favouring lateral displacement from the centre. Figure 5(c) highlights the influence of the power-law index
$n$
for
$Cu = 10$
,
$\beta = 0.5$
and
${\textit{Re}}_c = 150$
. As
$n$
decreases, corresponding to stronger shear-thinning behaviour, the centreline equilibrium is progressively destabilised, and stable off-centre equilibria emerge. This shows that it is not only the magnitude of shear thinning (via
$Cu$
), but that also the intrinsic rheological character of the fluid (via
$\beta$
and
$n$
) critically governs the equilibrium positions of the particle. Collectively, these results reveal a consistent and unified trend: increasing the intensity of shear thinning, whether by increasing
$Cu$
, decreasing
$\beta$
or reducing
$n$
, leads to an earlier onset of lateral migration and the formation of stable off-centre equilibrium positions. By systematically varying one parameter at a time, the lift-force profiles provide a clear and quantitative picture of how shear thinning modulates particle stability, highlighting the intricate interplay between fluid rheology and inertial effects in confined flows. To place these findings in context, it is useful to compare them with Newtonian fluids. Previous studies have shown that the lift force decreases with increasing Reynolds number, a phenomenon attributed to inertial screening (Fox, Schneider & Khair Reference Fox, Schneider and Khair2021). As Reynolds number increases, the disturbance flow induced by the particle becomes increasingly confined to its vicinity, thereby reducing hydrodynamic interactions with the channel walls and leading to a weaker lift force. In shear-thinning fluids at fixed
${\textit{Re}}_c$
, we observe a qualitatively similar reduction in lift force on increasing
$Cu$
, decreasing
$\beta$
or reducing
$n$
(see figure 5
a–c), although the underlying mechanism is different. Here, the shear-thinning rheology reduces the viscosity in regions of high shear rate near the particle surface, thereby diminishing local viscous stresses. This rheological effect confines the disturbance flow more closely to the particle, effectively screening particle–wall interactions which in turn reduces the overall lift force. These insights emphasise that even subtle changes in the fluid’s shear-rate-dependent properties can markedly alter the behaviour of suspended particles, offering a predictive framework for understanding particle migration in complex fluids.
Contours of dimensionless viscosity at steady state in plane Couette flow for different Carreau numbers, with
$\beta = 0.1$
,
$\kappa = 0.2$
,
$n = 0.25$
and
${\textit{Re}}_c = 150$
. Panels show (a)
$Cu = 0.1$
, (b)
$Cu = 1.0$
, (c)
$Cu = 10$
and (d)
$Cu = 100$
.
Next, to examine the reduction in the critical Reynolds number
${\textit{Re}}_{{cr}}$
with increasing
$Cu$
and decreasing
$\beta$
and
$n$
, kinematic viscosity
$\nu (x,z)$
contours are plotted, where the viscosity is normalised by the zero-shear-rate viscosity
$\nu _{0}$
. Figure 6(a–d) presents steady-state viscosity contours for increasing Carreau numbers,
$Cu = 0.1,\, 1,\, 10$
and
$100$
, at fixed
$\beta = 0.1$
,
$n = 0.25$
and
${\textit{Re}}_c = 150$
. These contours illustrate how varying
$Cu$
affects the local viscosity field around the particle and, consequently, its equilibrium. At low Carreau number (
$Cu = 0.1$
, figure 6
a), the fluid exhibits very weak shear thinning, resulting in a nearly uniform viscosity field across the channel. In this regime, viscosity is only slightly reduced in the vicinity of the particle due to localised shear. Consequently, the particle remains near the channel centreline, as the weak and nearly symmetric viscosity gradients generate negligible lift asymmetry. As
$Cu$
increases (figure 6
b–d), the shear-thinning response becomes progressively stronger, leading to a pronounced reduction in viscosity near the particle surface, where local velocity gradients are largest. This reduction in viscosity locally around the particle modifies the local stress distribution which eventually destabilises the centreline equilibrium, causing the particle to migrate towards an off-centre equilibrium position, the exact location of which depends on the initial particle position. The destabilisation of the centreline can be rationalised through the Carreau constitutive relation: larger local velocity gradients induce a stronger reduction in viscosity near the particle surface, which decreases the local viscous force. This reduction effectively makes local inertial force stronger than the viscous force, thereby accelerating the onset of inertial migration away from the centreline. At low
$Cu$
, where the fluid behaviour is nearly Newtonian, the centreline remains a stable equilibrium. However, as
$Cu$
increases, the enhanced shear thinning generates progressively stronger viscosity gradients, which in turn promote off-centre migration. This mechanism is reflected quantitatively in the critical Reynolds number,
${\textit{Re}}_{{cr}}$
, for lateral migration. For instance, at
$\beta = 0.9$
(figure 4
a), the Newtonian limit (
$Cu \to 0$
) yields
${\textit{Re}}_{{cr}} \approx 153$
, whereas at
$Cu = 10$
,
${\textit{Re}}_{{cr}}$
decreases markedly to
${\textit{Re}}_{{cr}}\approx 96$
. This reduction demonstrates that local shear thinning lowers the effective viscous resistance, facilitating an earlier onset of off-centre equilibrium positions as
$Cu$
increases. Overall, the viscosity contours provide a clear visualisation of how shear thinning modulates the local stress field around the particle, thereby governing its migration behaviour in confined flows.
Contours of the local Reynolds number
${\textit{Re}}_{\ell }(\boldsymbol{x}) = V_{w}H/\nu (x,z)$
on the mid-plane for
$Cu = 0.1$
,
$1$
,
$10$
and
$100$
at
${\textit{Re}}_c = 150$
,
$\beta = 0.1$
,
$n = 0.25$
and
$\kappa = 0.2$
(panels a–d). As
$Cu$
increases, shear thinning enhances
${\textit{Re}}_{\ell }$
in the vicinity of the particle through a reduction of viscous stresses near the particle surface.
To rationalise this migration using the local Reynolds number distribution, we quantify the spatial variation of
${\textit{Re}}_{\ell }(\boldsymbol{x})$
, following Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2012) and Patel, Rothstein & Modarres-Sadeghi (Reference Patel, Rothstein and Modarres-Sadeghi2022), in the mid-plane (i.e.
$xz$
-plane) passing through the particle centre. The local Reynolds number is defined as
${\textit{Re}}_{\ell }(\boldsymbol{x}) = V_w H / \nu (x,z)$
, which provides a local measure of the balance between inertial and viscous stresses and can therefore be interpreted as an effective local Reynolds number. Contours of
${\textit{Re}}_{\ell }(\boldsymbol{x})$
at steady state for
$Cu = 0.1$
,
$1$
,
$10$
and
$100$
are shown in figure 7(a–d) for
${\textit{Re}}_c = 150$
,
$\beta = 0.1$
and
$n = 0.25$
. As the Carreau number increases from
$Cu = 0.1$
to
$100$
,
${\textit{Re}}_{\ell }$
increases both around the particle surface and in the bulk region far away from the particle. Away from the particle, the distribution of
${\textit{Re}}_{\ell }$
is nearly uniform. This uniform far-field value corresponds to the Reynolds number based on the effective base viscosity associated with the imposed background shear rate according to the Carreau relation. The increase in the far-field
${\textit{Re}}_{\ell }$
with increasing
$Cu$
is due to the reduction in the base viscosity for a fixed background shear rate in the absence of the particle, as predicted by the Carreau rheology. This behaviour is also evident in the viscosity contour plots shown in figure 6(a–d). Closer to the particle surface, however, the distribution of
${\textit{Re}}_{\ell }(\boldsymbol{x})$
shows a significant local increase compared with its far-field value for each Carreau number. This local amplification arises from the enhanced shear rates generated around the particle, which further reduce the viscosity in shear-thinning fluids. Consequently, the local viscous stresses decrease, leading to a locally more inertia-dominated region in the vicinity of the particle. Although the global channel Reynolds number
${\textit{Re}}_c$
is fixed, the particle experiences a locally higher effective Reynolds number over a finite region surrounding it. In comparison with the Newtonian case at the same
${\textit{Re}}_c$
, the shear-thinning-induced reduction in viscous stress causes the particle to effectively experience a higher Reynolds number flow in its immediate neighbourhood. As a result, the equilibrium position shifts away from the centreline. For example, configurations that remain centred in the Newtonian limit shift to an off-centre equilibrium position at finite
$Cu$
, as shown in figure 4(a). Specifically, while the Newtonian case (
$Cu \to 0$
) yields
${\textit{Re}}_{{cr}} \approx 153$
, increasing the Carreau number to
$Cu = 10$
reduces the critical Reynolds number to
${\textit{Re}}_{{cr}} \approx 96$
, consistent with the enhanced local inertial effects discussed above. Overall these results indicates that even in the absence of flow curvature, shear thinning alone can generate localised inertia-dominated regions around the particle, destabilising the centreline equilibrium and causing off-centre migration.
Classical studies by Ho & Leal (Reference Ho and Leal1974) and Schonberg & Hinch (Reference Schonberg and Hinch1989) examined inertial migration of particles in shear and Poiseuille flows at small but finite Reynolds numbers using matched asymptotic expansions. In their framework, the flow is divided into distinct regions, one in the immediate vicinity of the particle (called the Stokes region) where viscous and pressure stresses dominate and the Stokes equations are valid, and the other far from the particle (called the Oseen region) where inertial stresses dominate. Consistency of the asymptotic analysis requires an intermediate region in which viscous and inertial stresses are comparable, which enables matching of the inner and outer solutions. In the absence of such a region, the matched asymptotic expansion approach breaks down. Further, Hood et al. (Reference Hood, Lee and Roper2015) revisited this problem numerically and found that, even at moderate channel Reynolds numbers (
${\textit{Re}}_c=10$
), viscous and pressure stresses exceed inertial contributions at all radial distances from the particle, leaving no region of co-dominance and contradicting the predictions of Ho & Leal (Reference Ho and Leal1974) and Schonberg & Hinch (Reference Schonberg and Hinch1989). At higher channel Reynolds numbers (
${\textit{Re}}_c=50$
and
$80$
), inertial stresses increase but scale linearly with
${\textit{Re}}_c$
, acting as a passive correction rather than establishing a new dominant balance. Further, in Newtonian shear flows, particle migration at finite Reynolds numbers has been interpreted in terms of inertial screening (Ho & Leal Reference Ho and Leal1974; Fox et al. Reference Fox, Schneider and Khair2021; Anand & Subramanian Reference Anand and Subramanian2023). Fox et al. (Reference Fox, Schneider and Khair2021) also showed that increasing particle Reynolds number confines the particle-induced velocity disturbance to a thin region near the particle surface, thereby weakening long-range hydrodynamic interactions and particle–wall coupling. Existing studies (Ho & Leal Reference Ho and Leal1974; Anand & Subramanian Reference Anand and Subramanian2023) show that, for ambient shear flow, the characteristic inertial screening length scales as
$\ell _s \sim a {\textit{Re}}_{\!p}^{-1/2}$
, or equivalently in channel scaling,
$\ell _s \sim H {\textit{Re}}_c^{-1/2}$
. When
${\textit{Re}}_c \ll 1$
,
$\ell _s \gg H$
and the walls lie within the inner Stokes region, with inertia acting as a regular perturbation. When
${\textit{Re}}_c = O(1)$
or larger,
$\ell _s \lesssim H$
and the walls lie in the Oseen region. In the present case, one further insight may be drawn. Since
${\textit{Re}}_c \gt O(1)$
, it follows that
$\ell _s \lesssim H$
, and the confining walls lie in the Oseen region. Consequently, inertial effects influence the global disturbance field rather than acting solely as a regular perturbation. Furthermore, owing to shear thinning, the viscosity decreases with increasing
$Cu$
, leading to a reduction in the effective base viscosity at the imposed shear rate as per Carreau fluid rheology. As a result, even far from the particle, the effective local Reynolds number increases in correspondence with the channel Reynolds number
${\textit{Re}}_c$
. In addition, the shear rate is significantly elevated in the immediate vicinity of the particle, producing a further local reduction in viscosity and a pronounced enhancement of the local Reynolds number, as demonstrated by the spatial distributions of
${\textit{Re}}_{\ell }$
shown in figure 7. Notably, this near-particle amplification of
${\textit{Re}}_{\ell }$
exceeds its far-field value. Such a localised enhancement of inertia is consistent with a stronger spatial localisation of the particle-induced disturbance. Consequently, relative to a Newtonian fluid at the same nominal channel Reynolds number, the effective inertial screening length,
$\ell _s$
, is further reduced in the shear-thinning case. This reduction in the screening length, together with the associated attenuation of the disturbance field, is reflected in the variation of the non-dimensional lift force across the channel height, as shown in figure 5(a), obtained at a fixed channel Reynolds number
${\textit{Re}}_c = 150$
with
$\kappa$
,
$n$
and
$\beta$
held constant. As
$Cu$
increases from
$0$
(Newtonian) to
$100$
, the magnitude of the lift force decreases at all fixed positions
$z/H$
, leading to progressively flatter lift-force profiles at higher
$Cu$
. This reduction in lift is consistent with enhanced inertial screening of the particle-induced velocity disturbance. At larger
$Cu$
, the disturbance becomes increasingly confined to the immediate vicinity of the particle surface, thereby weakening long-range hydrodynamic interactions between the particle and the confining wall, particularly wall-mediated interactions responsible for lift are weakened and reducing the net lift force. A similar trend is observed upon decreasing
$\beta$
(figure 5
b) or decreasing the power-law index
$n$
(figure 5
c), both of which intensify shear thinning. In each case, stronger shear thinning leads to a further reduction in the lift force, consistent with a greater spatial localisation of the particle-induced disturbance. Overall, our results show that, even in the absence of flow curvature, shear thinning alone is sufficient to destabilise the centreline equilibrium and trigger off-centre inertial migration. This transition can occur at relatively low channel Reynolds numbers, as shear thinning reduces the critical Reynolds number
${\textit{Re}}_{{cr}}$
compared with a Newtonian fluid. The mechanism arises from the formation of localised inertia-dominated regions around the particle, where reduced viscous stresses and increased inertial effects produce a significant local rise in
${\textit{Re}}_{\ell }$
. The resulting migration behaviour differs fundamentally from that observed in Newtonian channel flows. Overall these findings highlight the important role of the shear-thinning rheology in modifying the local stress balance and in turn affecting the nature of inertial migration.
Time-dependent trajectories of a particle released from different initial positions. (a) Effect of the Carreau number at
$n = 0.25$
,
$\kappa = 0.2$
,
$\beta = 0.1$
and
${\textit{Re}}_c = 150$
. (b) Effect of the viscosity ratio at
$Cu = 10$
,
$n = 0.25$
,
$\kappa = 0.2$
and
${\textit{Re}}_c = 150$
. (c) Effect of the power-law index at
$Cu = 10$
,
$\beta = 0.5$
,
$\kappa = 0.2$
and
${\textit{Re}}_c = 150$
. (d) Influence of the Carreau number,
$Cu$
, on the particle migration time. The evolution of the normalised particle position,
$Z/H$
, with
$Gt$
illustrates the time taken to attain the equilibrium position. Results correspond to
${\textit{Re}}_c = 50$
and
$\beta = 0.5$
,
$n = 0.25$
for
$\kappa = 0.2$
.
Figure 8(a–c) presents the lateral migration trajectories of particles released from two different initial positions,
$Z_0/H=0.4$
(below the centreline) and
$Z_0/H=0.6$
(above the centreline), within the channel. The particles are free to translate and rotate, and their trajectories are followed until they reach steady-state equilibrium positions in the
$z$
-direction. In figure 8(a), the Carreau number is varied while fixing the other parameters as
$\beta =0.1$
,
$n=0.25$
and
${\textit{Re}}_c=150$
. For
$Cu=1$
, the particle migrates to the centreline (
$Z/H=0.5$
), independent of the initial position. As
$Cu$
increases, however, two stable equilibria appear symmetrically about the channel centreline; a particle released from
$Z_0/H=0.4$
migrates to an off-centre position between the lower wall and the centreline, whereas a particle released from
$Z_0/H=0.6$
stabilises symmetrically between the centreline and the upper wall. This demonstrates that sufficiently strong non-Newtonian effects break the uniqueness of the centreline equilibrium, and the final state becomes sensitive to the initial release height. Figure 8(b) illustrates the role of the viscosity ratio
$\beta$
for
$Cu=10$
,
$n=0.25$
and
${\textit{Re}}_c=150$
. At larger values of
$\beta$
, the centreline remains the sole equilibrium position, and trajectories are independent of initial lateral positions. Reducing
$\beta$
, however, produces two off-centre equilibria symmetrically about the channel centreline, such that a particle released near the lower wall remains trapped below the centreline while one released near the upper wall stabilises above it. A similar trend is observed when varying the power-law index
$n$
(figure 8
c) at fixed
$Cu=10$
,
$\beta =0.5$
and
${\textit{Re}}_c=150$
: for
$n=1$
(Newtonian case), particles invariably migrate to the centreline, but as
$n$
decreases, corresponding to stronger shear-thinning behaviour, the centreline equilibrium becomes unstable and stable off-centre equilibria emerge. To further understand the migration dynamics in comparison with the Newtonian case, we consider a parameter set for which the particle migrates to the centreline, namely
${\textit{Re}}_c = 50$
,
$\beta = 0.5$
and
$n = 0.25$
, as shown in figure 8(d). It can be observed that as the Carreau number increases from
$Cu = 0.1$
to
$100$
, the time required for the particle to reach the centreline decreases. In addition to that, we have observed that the streamwise distance travelled by the particle before reaching its equilibrium position is also reduced. This behaviour indicates that in a shear-thinning fluid the particle experiences a lower effective hydrodynamic force than its Newtonian counterpart. The reduction in viscosity with increasing shear rate decreases the overall viscous resistance acting on the particle, allowing it to migrate more rapidly and over a shorter streamwise distance toward the centreline. Taken together, these trajectories provide direct numerical confirmation of the multiple stable equilibrium positions predicted by the lift-force analysis shown in figure 5(a–c). The viscosity contours (figure 6) and the local Reynolds number distributions (figure 7) further clarify the underlying mechanism, demonstrating that shear thinning enhances local inertial effects through viscosity reduction in high-shear regions around the particle. As a consequence, the migration dynamics and equilibrium locations differ fundamentally from the Newtonian case, and particle migration is governed primarily by the coupled influence of fluid rheology and inertia. All the above results correspond to a confinement ratio
$\kappa = 0.2$
. It is therefore important to examine how these equilibrium states and migration characteristics are modified as the degree of confinement is varied. This aspect is addressed in the following section.
3.2. Effect of confinement ratio
To elucidate the role of geometric confinement on inertial migration, we examine the effect of the confinement ratio
$\kappa = a/H$
, which quantifies the particle size relative to the channel height. In the present study,
$\kappa$
is varied between
$0.05$
and
$0.2$
, corresponding to weak confinement (
$\kappa = 0.05, 0.1$
) and strong confinement (
$\kappa = 0.2$
). Figure 9(a) shows the normalised stable equilibrium position as a function of the Reynolds number
${\textit{Re}}_c$
for different confinement ratios at fixed rheological parameters
$Cu = 10$
,
$\beta = 0.5$
and
$n = 0.25$
. It is evident that the critical Reynolds number
${\textit{Re}}_{{cr}}$
increases with increasing confinement. In particular, for strong confinement (
$\kappa = 0.2$
), the onset of off-centre migration occurs at a significantly higher
${\textit{Re}}_c$
compared with the weakly confined cases (
$\kappa = 0.05, 0.1$
), indicating a stabilising influence of the channel walls. The influence of shear-thinning at a fixed confinement is illustrated in figure 9(b), which compares equilibrium positions for two Carreau numbers at
$\kappa = 0.1$
, with
$\beta = 0.5$
and
$n = 0.25$
. At low
${\textit{Re}}_c$
, the particle remains centred; however, beyond a critical Reynolds number
${\textit{Re}}_{{cr}}$
, a symmetry-breaking bifurcation occurs and stable off-centre equilibrium positions emerge. Importantly, increasing the Carreau number leads to a systematic reduction in
${\textit{Re}}_{{cr}}$
, consistent with the trends observed under stronger confinement. This demonstrates that shear thinning accelerates the onset of inertial migration irrespective of the confinement level. Figure 9(c) summarises these trends by showing the variation of the critical Reynolds number
${\textit{Re}}_{{cr}}$
with the Carreau number
$Cu$
for confinement ratios
$\kappa = 0.05, 0.1$
and
$0.2$
, at fixed rheological parameters
$n = 0.25$
and
$\beta = 0.5$
. In the limit of very small
$Cu$
, corresponding to a flow approaching the Newtonian regime,
${\textit{Re}}_{{cr}}$
decreases with decreasing confinement ratio and approaches closer to
${\textit{Re}}_{{cr}} \approx 148$
specifically when
$\kappa = 0.05$
. This value is in good agreement with the analytically obtained critical Reynolds number reported by Anand & Subramanian (Reference Anand and Subramanian2023) for a neutrally buoyant particle in a Newtonian fluid. As the Carreau number increases,
${\textit{Re}}_{{cr}}$
decreases markedly, indicating an enhanced role of shear thinning in promoting the instability. Furthermore, for a fixed value of
$Cu$
, reducing the confinement ratio, equivalently increasing the channel height relative to the particle size, leads to a further reduction in
${\textit{Re}}_{{cr}}$
. Overall, the dependence of
${\textit{Re}}_{{cr}}$
on
$\kappa$
reflects a balance between wall-induced hydrodynamic interactions and viscosity-driven asymmetries. Under strong confinement, enhanced particle–wall interactions generate stabilising wall-induced lift forces that delay the onset of bifurcation. In contrast, weak confinement diminishes wall effects, allowing non-Newtonian modifications of the flow field to dominate the lateral force balance and trigger earlier symmetry breaking. It is further observed from figure 9(c) that, at lower Carreau numbers, corresponding to the near-Newtonian limit, the variation of the critical Reynolds number
${\textit{Re}}_{{cr}}$
with confinement ratio is relatively weak. However, as
$Cu$
increases, a pronounced dependence of
${\textit{Re}}_{{cr}}$
on confinement emerges, indicating that confinement effects become significantly amplified in the shear-thinning regime.
(a) Normalised equilibrium position with channel Reynolds number
${\textit{Re}}_c$
for three confinement ratios,
$\kappa = 0.05$
,
$0.1$
and
$0.2$
, at
$\beta = 0.5$
,
$Cu = 10$
and
$n = 0.25$
. (b) Influence of Carreau number
$Cu$
on equilibrium position for
$\kappa = 0.1$
,
$\beta = 0.5$
,
$n = 0.25$
, showing the effect of shear-thinning intensity at fixed confinement. (c) Dependence of the critical Reynolds number
${\textit{Re}}_{{cr}}$
on the Carreau number
$Cu$
for different confinement ratios at
$\kappa$
for
$n = 0.25$
and
$\beta = 0.5$
.
(a) Variation of the normalised particle translational velocity in the streamwise direction,
$U/V_w$
, with the Carreau number,
$Cu$
, at
${\textit{Re}}_c = 150$
,
$\beta = 0.1$
and
$n = 0.25$
, for different confinement ratios
$\kappa$
. The particle is initially released from the lower half of the channel. (b) Corresponding variation of the normalised particle angular velocity in
$z$
direction,
$\varOmega /G$
, under the same conditions.
Next, we look at the translational and angular velocities of the particle at steady state. Figure 10(a, b) shows the magnitude of particle translational velocity
$(U=|\boldsymbol{U}|)$
normalised by the wall velocity
$V_w$
, and the magnitude of angular velocity
$(\varOmega = |\boldsymbol{\varOmega }|)$
normalised by the imposed shear rate
$G$
, for varying Carreau numbers as a function of confinement ratio
$\kappa$
. The results correspond to
${\textit{Re}}_c = 150$
,
$\beta = 0.1$
and
$n = 0.25$
. It is observed that both the translational and angular velocities decrease with increasing
$Cu$
. Furthermore, as the confinement ratio
$\kappa$
decreases (i.e. for wider channels), the translational velocity decreases, whereas the rotational velocity increases. This behaviour is linked to the shift in the equilibrium position of the particle. For the present parameter set, all cases except
$\kappa = 0.2$
at
$Cu = 0.1$
migrate toward off-centre equilibrium positions between the bottom wall and the centreline, since the particle is initially released in the lower half of the channel. As
$Cu$
increases and
$\kappa$
decreases, the equilibrium position moves progressively closer to the wall. The reduced local streamwise fluid velocity near the wall leads to a decrease in translational speed, while the enhanced velocity gradients in the wall vicinity promote increased particle rotation.
4. Conclusion
We have numerically investigated the inertial migration of a neutrally buoyant spherical particle in plane Couette flow of inelastic shear-thinning fluid governed by the Carreau model. The influence of key physical parameters, including the channel Reynolds number (
${\textit{Re}}_c$
), Carreau number (
$Cu$
), viscosity ratio (
$\beta$
), power-law index (
$n$
) and confinement ratio (
$\kappa$
), has been systematically examined to assess their impact on lateral migration and equilibrium positioning over a wide range of flow conditions. In Newtonian fluids, particles undergoing inertial migration remain stably centred at low Reynolds numbers, up to a critical channel Reynolds number,
${\textit{Re}}_{{cr}}$
, beyond which two symmetric off-centre equilibrium positions emerge through a supercritical pitchfork bifurcation (Anand & Subramanian Reference Anand and Subramanian2023). In contrast, for shear-thinning fluids, this bifurcation occurs at significantly lower values of the channel Reynolds number (i.e.
${\textit{Re}}_{{cr}}$
reduces), leading to an earlier transition from centreline to off-centre equilibria that remain symmetric about the channel centreline. This shift in the bifurcation threshold is strongly influenced by the degree of shear thinning, characterised by larger
$Cu$
and smaller values of
$\beta$
and
$n$
. Enhanced shear thinning reduces the local viscosity in the vicinity of the particle, resulting in a reduction of the viscous stress which leads to increase in the local Reynolds number around the particle. These rheology-induced modifications alter the distribution of stresses around the particle and reshape the lateral inertial lift-force profile, thereby destabilising the centreline equilibrium. It is also worth noting that, at fixed
${\textit{Re}}_c$
, increasing the Carreau number drives the particle progressively closer to the channel wall, further highlighting the role of spatially varying viscosity in governing the migration dynamics. Overall, the migration behaviour exhibits a pronounced dependence on
$Cu$
,
$\beta$
and
$n$
, with the destabilisation of the centreline equilibrium and the emergence of stable off-centre equilibria occurring under appropriate combinations of rheological parameters. A detailed force-balance analysis reveals that the stability of an equilibrium position is determined by the slope of the lateral force near its zero crossing; a negative slope corresponds to a stable equilibrium, whereas a positive slope indicates unstable equilibrium. This provides a clear physical interpretation of how the shear-thinning rheology modifies the inertial lift-force landscape through its interplay with inertia. Furthermore, the effect of geometric confinement, quantified by the confinement ratio
$\kappa$
, is also examined and identified as a critical factor in inertial migration, as it strongly influences the balance between wall-induced hydrodynamic interactions and inertia-driven lateral forces acting on the particle. Weak confinement (smaller
$\kappa$
), corresponding to larger channel heights relative to the particle diameter, reduces wall-induced interactions and promotes an earlier onset of equilibrium bifurcation by allowing viscosity-driven asymmetries to dominate the lateral force balance. In contrast, strong confinement enhances wall-induced lift forces, stabilising the centreline trajectory and delaying or suppressing the emergence of off-centre equilibria. Another notable observation is that the reduction in the critical Reynolds number remains weak with decreasing confinement ratio at low Carreau numbers, i.e. in the near-Newtonian limit. However, as the Carreau number increases, this reduction becomes significantly more pronounced, indicating that the combined effects of shear thinning and geometric confinement strongly influence the stability threshold.
It is worthwhile commenting on the role of fluid elasticity, which has not been considered in the present study. The analysis here is restricted to inelastic shear-thinning fluids, whereas for many biological and polymeric liquids viscoelasticity plays a central role and can qualitatively modify particle migration. Extending the present framework to incorporate both a shear-dependent viscosity and fluid elasticity thus represents a natural and important direction for future work. The destabilisation mechanism identified here differs fundamentally from that expected in viscoelastic fluids. In shear-thinning media, the reduction of viscosity in high-shear regions surrounding the particle enhances local inertial stresses, leading to an earlier loss of stability of the centreline equilibrium. Migration is therefore promoted indirectly through inertia amplified by spatial variations in viscosity. By contrast, viscoelastic fluids generate normal-stress differences that can drive cross-stream migration even in the absence of inertia. In plane Couette flow, numerical and experimental studies have demonstrated that elastic normal stresses render the centre plane unstable and drive particles toward the nearest wall (D’Avino et al. Reference D’Avino, Maffettone, Greco and Hulsen2010; Caserta et al. Reference Caserta, D’Avino, Greco, Guido and Maffettone2011). At small Deborah numbers, the migration velocity scales linearly with the transverse displacement, reflecting a normal-stress imbalance across the particle surface. Elasticity alone can therefore produce wall-directed migration through a purely stress-driven mechanism. At finite Reynolds numbers, however, inertia and shear-thinning effects both act to destabilise the centreline equilibrium, while elasticity independently promotes wall-ward motion in the case of shear flow. The combined action of these mechanisms may therefore reinforce one another, potentially displacing particles further away from the centreline and toward the walls, depending on fluid and flow parameters. To the best of our knowledge, this coupled elasto-inertial regime in plane Couette flow at moderate to high Reynolds numbers has not yet been systematically investigated, and remains an open and compelling problem.
Lastly, it should be noted that the present numerical framework does not disentangle the individual contributions of wall-induced lift and shear-thinning-induced lift. In confined shear flows, wall proximity, background shear and rheological variations are inherently coupled, and the observed migration behaviour emerges from their combined action. As a result, it is not possible to independently quantify the respective roles of geometric confinement and viscosity stratification in generating the net transverse force. A systematic decomposition of these effects would require the development of an analytical framework capable of isolating distinct physical mechanisms asymptotically. Classical theoretical treatments of inertial lift in wall-bounded shear flows (Ho & Leal Reference Ho and Leal1974; Schonberg & Hinch Reference Schonberg and Hinch1989; Hood et al. Reference Hood, Lee and Roper2015; Anand & Subramanian Reference Anand and Subramanian2023) explicitly incorporate the presence of boundaries, typically through reflection procedures or matched asymptotic expansions that satisfy no-slip conditions at the walls. These studies therefore quantify inertial lift in the presence of wall effects and shear-profile curvature. In contrast, foundational analyses such as those of Saffman (Reference Saffman1965) and McLaughlin (Reference McLaughlin1991) consider particles in unbounded linear shear flows and derive inertial lift forces in the absence of confining boundaries. The methodology employed in these works provides a framework for isolating shear-induced inertial lift independent of wall interactions. A comparable analytical treatment adapted to shear-thinning channel flows would allow a clearer separation of wall-induced and rheology-induced contributions, and represents a promising direction for future investigation. Overall, the present results provide new insight into rheologically modulated inertial migration in confined shear flows and highlight the coupled influence of shear-thinning behaviour and geometric confinement on particle trajectories. These findings may offer useful guidance for the design and optimisation of microfluidic devices intended for particle focusing, separation and transport in complex fluids.
Acknowledgements
The authors gratefully acknowledge the PARAM Pravega high-performance computing system, located at the Supercomputer Education and Research Centre of the Indian Institute of Science.
Funding
The authors acknowledge the generous support of the Saroj Poddar Trust and the Anusandhan National Research Foundation (Grant No. SRG/2021/000074).
Declaration of interest
The authors declare no conflict of interest.
Appendix A. Algorithmic details of the SP–LBM implementation
The LBM coupled with the SPM (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005) is employed to simulate the interaction between a rigid spherical particle and a shear-thinning fluid. The fluid flow is resolved on a fixed Eulerian lattice using the LBM, while the particle dynamics is described in a Lagrangian framework. The coupling between the fluid and particle phases is achieved through a body-force formulation based on the SPM (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Jafari et al. Reference Jafari, Yamamoto and Rahnama2011). The inertial (added-mass) effects are incorporated by following Mino et al. (Reference Mino2017). We have implemented the algorithm on CUDA-capable GPU. For clarity and reproducibility, the algorithm is summarised as follows.
Step 1: initialisation
At time
$ t=0$
, both the fluid and the particle are assumed to be at rest. The following quantities are initialised:
\begin{align} &\text{Particle position:} && \boldsymbol{X}(0) = \boldsymbol{X}_{0}, \nonumber \\ &\text{Particle translational velocity:} && \boldsymbol{U}(0) = \boldsymbol{0}, \nonumber \\ &\text{Particle angular velocity:} && \boldsymbol{\varOmega }(0) = \boldsymbol{0}, \nonumber \\ &\text{Fluid velocity:} && \boldsymbol{u}(\boldsymbol{x},0) = \boldsymbol{0}, \nonumber \\ &\text{Fluid density:} && \rho (\boldsymbol{x},0) = 1, \nonumber \\ &\text{Hydrodynamic force:} && \boldsymbol{F}_{\!H}(0) = \boldsymbol{0}, \nonumber \\ &\text{Hydrodynamic torque:} && \boldsymbol{T}_{\!H}(0) = \boldsymbol{0}, \nonumber \\ &\text{Distribution functions:} && f_i(\boldsymbol{x},0) = f_i^{\textit{eq}}. \end{align}
Step 2: particle velocity and position update
The translational and angular velocities of the particle are updated using the hydrodynamic forces and torques, including the inertial compensation terms
\begin{align} \boldsymbol{U}(t+\Delta t) &= \boldsymbol{U}(t) + \frac {1}{M_{\!p}} \left [\boldsymbol{F}_{\!H}(t)+\boldsymbol{F}_{\!H,\textit{in}}(t) \right ]\,\Delta t, \end{align}
\begin{align} \boldsymbol{\varOmega }(t+\Delta t) &= \boldsymbol{\varOmega }(t) + \frac {1}{I_{\!p}} \left [\boldsymbol{T}_{\!H}(t)+\boldsymbol{T}_{H,\textit{in}}(t)\right ]\Delta t. \end{align}
The particle position is updated using a Crank–Nicolson scheme
\begin{equation} \boldsymbol{X}(t+\Delta t) = \boldsymbol{X}(t) + \frac {1}{2} \left [\boldsymbol{U}(t+\Delta t)+\boldsymbol{U}(t) \right ]\,\Delta t. \end{equation}
Step 3: smoothed profile function
Following the particle update, the diffuse particle–fluid interface is described using the smoothed profile function
\begin{equation} \phi (\boldsymbol{x},t+\Delta t) = \frac {1}{2} \left [ 1+ \tanh \!\left ( \frac {a-\left |\boldsymbol{x}-\boldsymbol{X}(t+\Delta t)\right |}{\xi } \right ) \right ]\!, \end{equation}
where
$a = 30 \Delta x$
is the particle radius and
$\xi = 2 \Delta x$
denotes the interfacial thickness.
Step 4: particle velocity field
The rigid-body velocity field inside the particle domain is evaluated as
\begin{equation} \phi (\boldsymbol{x},t+\Delta t)\boldsymbol{u}_{\!p}(\boldsymbol{x},t+\Delta t) = \phi (\boldsymbol{x},t+\Delta t) \!\left [ \boldsymbol{U}(t+\Delta t) + \boldsymbol{\varOmega }(t+\Delta t) \times \!\bigl (\boldsymbol{x}-\boldsymbol{X}(t+\Delta t)\bigr ) \right ]\!. \end{equation}
Step 5: equilibrium distribution and forcing term
The equilibrium distribution function is given by
\begin{equation} f_i^{\textit{eq}} = \rho \,w_i \left [ 1+ \frac {\boldsymbol{c}_i\boldsymbol{\cdot }\boldsymbol{u}}{c_s^2} + \frac {(\boldsymbol{c}_i\boldsymbol{\cdot }\boldsymbol{u})^2}{2c_s^4} - \frac {\boldsymbol{u}^2}{2c_s^2} \right ]\!, \end{equation}
where
$ w_i$
are the lattice weights and
$ c_s$
is the lattice speed of sound. External and particle-induced forces are incorporated using Guo’s forcing scheme (Guo et al. Reference Guo, Zheng and Shi2002), in which the discrete forcing term is expressed as
\begin{equation} F_i = w_i \left [ \frac {\boldsymbol{c}_i-\boldsymbol{u}}{c_s^2} + \frac {(\boldsymbol{c}_i\boldsymbol{\cdot }\boldsymbol{u})\boldsymbol{c}_i}{c_s^4} \right ] \boldsymbol{\cdot }\boldsymbol{F}. \end{equation}
For shear-flow configurations considered in the present study, the total body-force density is expressed as
$ \boldsymbol{F}(\boldsymbol{x},t) = \phi (\boldsymbol{x},t+\Delta t)\,\boldsymbol{f}_{\!p},$
where
$ \phi (\boldsymbol{x},t+\Delta t)$
is the smoothed profile function distinguishing the particle and fluid regions, and
$ \boldsymbol{f}_{\!p}$
represents the force density enforcing particle rigidity and momentum exchange between the particle and the surrounding fluid.
Step 6: collision and streaming
The LB evolution is decomposed into two successive steps: collision and streaming. The collision process is local and instantaneous, and therefore does not advance the time level. During this step, the distribution function relaxes toward its equilibrium value under the action of the BGK collision operator and the forcing term. The post-collision distribution function, denoted by
$ f_i^{\prime}$
, is given by
\begin{equation} f_i^{\prime}(\boldsymbol{x},t) = f_i(\boldsymbol{x},t) - \frac {\Delta t}{\tau } \left [ f_i(\boldsymbol{x},t)-f_i^{\textit{eq}}(\boldsymbol{x},t) \right ] + F_i(\boldsymbol{x},t)\,\Delta t, \end{equation}
where
$ f_i^{\prime}(\boldsymbol{x},t)$
represents the post-collision (pre-streaming) distribution function. In the subsequent streaming step, the post-collision distributions are propagated along the corresponding lattice velocity directions
$ \boldsymbol{c}_i$
to neighbouring lattice nodes, which advances the solution in time
\begin{equation} f_i(\boldsymbol{x}+\boldsymbol{c}_i\Delta t,t+\Delta t) = f_i^{\prime}(\boldsymbol{x},t). \end{equation}
Following the streaming step, a mid-grid (half-way) bounce-back boundary condition (Ladd Reference Ladd1994) is imposed at the top and bottom walls. In this formulation, the wall is located midway between fluid and solid lattice nodes, such that particle populations effectively travel only half a lattice link before being reflected back along the opposite direction. The velocity reversal is thus realised within the streaming process, thereby enforcing the no-slip condition at the solid boundaries (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017).
Step 7: macroscopic fields and shear-thinning viscosity
The macroscopic density and the intermediate velocity field are computed as
\begin{align} \rho &= \sum _i f_i, \end{align}
\begin{align} \boldsymbol{u}^* &= \frac {1}{\rho } \left ( \sum _i f_i \boldsymbol{c}_i + \frac {1}{2}\boldsymbol{F}\Delta t \right )\!, \end{align}
where
$\boldsymbol{F}$
denotes the total body-force density acting on the fluid, including both external and particle-induced contributions. The rate of strain tensor is obtained locally from the non-equilibrium distribution functions, consistent with the inherent locality of the LBM and preserving second-order accuracy (Krüger et al. Reference Krüger, Varnik and Raabe2009) as
\begin{equation} S_{\alpha \beta } = -\frac {1}{2\tau c_s^2} \sum _i c_{i\alpha } c_{i\beta } f_i^{\textit{neq}}, \end{equation}
with
$f_i^{\textit{neq}} = f_i - f_i^{{\textit{eq}}}$
. This local formulation avoids finite-difference approximations of velocity gradients, thereby reducing numerical diffusion while retaining second-order accuracy (Krüger et al. Reference Krüger, Varnik and Raabe2009). The magnitude of local shear rate is defined as
\begin{equation} |\dot {\gamma }| = \sqrt {2 S_{\alpha \beta } S_{\alpha \beta }}. \end{equation}
For shear-thinning fluids, the dynamic viscosity
$\mu (\dot {\gamma })$
is described using the Carreau model
\begin{equation} \mu (\dot {\gamma }) = \mu _\infty + \left (\mu _0 - \mu _\infty \right ) \left [ 1 + \left (\lambda \dot {|\gamma |}\right )^2 \right ]^{\frac {n-1}{2}}, \end{equation}
where
$\mu _0$
and
$\mu _\infty$
are the zero- and infinite-shear viscosities, respectively,
$\lambda$
is the characteristic time constant and
$n\lt 1$
is the power-law index. The relaxation time is related to the local dynamic viscosity through
\begin{equation} \tau (\dot {\gamma }) = \dfrac {\mu (\dot {\gamma })}{\rho c_s^2 \Delta t} + \dfrac {1}{2} .\end{equation}
It is important to note that (A13) contains the relaxation time
$\tau$
, which itself depends on the local shear rate through (A16). The shear-rate magnitude, computed from the strain-rate tensor (see (A14)), therefore depends implicitly on
$\tau$
. Consequently, at a given time step, the evaluation of
$\tau$
via (A16) requires knowledge of the shear rate, which in turn depends on
$\tau$
through (A13). This renders the formulation formally implicit at each lattice node. One possible approach is to resolve this nonlinear coupling iteratively at every time step and at every node until convergence is achieved. However, such a fully implicit treatment substantially increases the computational cost. In the present study, we adopt a computationally efficient explicit strategy, whereby the strain-rate tensor at the current time step is evaluated by taking the right-hand side quantities in (A13) from the previous time step. The resulting shear rate is then used to update the viscosity and the relaxation time through (A16). This approximation decouples the nonlinear update procedure while retaining numerical stability and accuracy, and has been widely employed in non-Newtonian LB implementations (Phillips & Roberts Reference Phillips and Roberts2011; Wang & Ho Reference Wang and Ho2011).
Step 8: particle–fluid interaction force
The interaction force enforcing the no-slip condition at the particle interface is computed as
\begin{equation} \phi (\boldsymbol{x},t+\Delta t)\,\boldsymbol{f}_{\!p} = \phi (\boldsymbol{x},t+\Delta t)\, \frac {\rho }{\Delta t} \left [ \boldsymbol{u}_{\!p}(\boldsymbol{x},t+\Delta t)-\boldsymbol{u}^* \right ]. \end{equation}
Step 9: velocity correction
The macroscopic fluid velocity is corrected according to
\begin{equation} \boldsymbol{u} = \boldsymbol{u}^* + \frac {1}{2\rho } \phi (\boldsymbol{x},t+\Delta t)\,\boldsymbol{f}_{\!p}\,\Delta t. \end{equation}
Step 10: hydrodynamic force and torque
The hydrodynamic force and torque acting on the particle are obtained as
\begin{align} \boldsymbol{F}_{\!H}(t+\Delta t) &= -\frac {1}{\Delta t} \sum _{\boldsymbol{x}\in V_{\!p}} \phi (\boldsymbol{x},t+\Delta t)\,\boldsymbol{f}_{\!p}, \end{align}
\begin{align} \boldsymbol{T}_{\!H}(t+\Delta t) &= -\frac {1}{\Delta t} \sum _{\boldsymbol{x}\in V_{\!p}} \phi (\boldsymbol{x},t+\Delta t) \left ( \boldsymbol{x}-\boldsymbol{X}(t+\Delta t) \right ) \times \boldsymbol{f}_{\!p}, \end{align}
where
$ V_{\!p}$
denotes the particle volume.
Step 11: time advancement
The simulation advances to the next time step, and steps 2–10 are repeated until a steady state is reached.
Appendix B. Grid, time-step and domain-size verification
This appendix documents the numerical verification procedures adopted to ensure that the reported numerical simulation results are independent of grid resolution, time step and computational domain size. All verification tests are carried out for an extreme parameter set considered in the study, thereby providing a stringent assessment of numerical robustness.
B.1. Numerical considerations for LB simulations
The selection of numerical parameters in LB simulations is guided by considerations of accuracy, stability and computational efficiency. In the LBM, the Mach number
$Ma$
plays a role analogous to the Courant–Friedrichs–Lewy number (Wang et al. Reference Wang, Wang, Lallemand and Luo2013). The Mach number
$Ma = V_w/c_s$
, where
$V_w$
is a characteristic velocity and
$c_s$
is the lattice sound speed, must therefore remain sufficiently small (typically
$Ma \lesssim 0.1$
) to suppress compressibility effects (Krüger et al. Reference Krüger, Varnik and Raabe2009). A uniform Cartesian lattice is employed in the standard LBM, with identical spacing in all directions (
$\Delta x = \Delta y = \Delta z = 1$
). The lattice nodes are located at the centres of the grid cells and are separated by one lattice unit, so that the total number of nodes directly corresponds to the physical domain size in each direction. For a prescribed channel Reynolds number
${\textit{Re}}_c = V_w H/\nu _0$
, the Mach number, the relaxation time
$\tau _0$
(corresponding to the zero-shear-rate viscosity
$\mu _0$
) and the particle radius
$a$
(or equivalently the channel height
$H$
) determine the physical lattice spacing
$\Delta x_{\textit{ph}}$
and the physical time step
$\Delta t_{\textit{ph}}$
. Under diffusive scaling,
$\Delta x_{\textit{ph}} \propto Ma$
and
$\Delta t_{\textit{ph}} \propto Ma^2$
(Krüger et al. Reference Krüger, Varnik and Raabe2009; Luo et al. Reference Luo, Lallemand and d’Humières2011), ensuring convergence to the incompressible Navier–Stokes equations as
$\Delta x_{\textit{ph}} \rightarrow 0$
. Accordingly, maintaining
$Ma \lesssim 0.1$
,
$\tau _{0} \in [0.7,1.0]$
and a sufficiently large computational domain provides a practical balance between numerical stability, accuracy and computational cost. In the present study, grid-resolution, time-step and domain-aspect-ratio independence tests are performed for an extreme case considered for the study, namely
${\textit{Re}}_c = 300$
,
$Cu = 100$
,
$\beta = 0.1$
,
$n = 0.25$
and
$\kappa = 0.1$
, with the particle initially located at
$Z_0/H = 0.4$
. These tests confirm that the results are insensitive to further refinements of the lattice resolution, time step and domain size, ensuring that the chosen numerical parameters reliably resolve both particle and flow dynamics for all other cases reported in this study. Independence of the grid, aspect ratio and time step is summarised as follows.
Numerical independence studies for the extreme parameter set. (a) Grid-resolution independence for different particle radii
$a$
. (b) Time-step independence obtained by varying the wall velocity
$V_w$
while keeping
$a = 30$
fixed. (c) Aspect-ratio independence for
$AR = 1, 2$
and
$3$
with the particle radius fixed at
$a = 30$
. All cases correspond to
${\textit{Re}}_c = 300$
,
$Cu = 100$
,
$\beta = 0.1$
,
$n = 0.25$
and
$\kappa = 0.1$
.
B.2. Grid-resolution independence
Grid independence is assessed by varying the particle radius
$a$
while keeping the characteristic wall velocity fixed at
$V_w = 0.1$
and the domain aspect ratio at
$AR = 2$
, as shown in figure 11(a). The BGK relaxation time is maintained within the stable range
$0.7 \leqslant \tau _{0} \leqslant 0.9$
, ensuring numerical stability and second-order accuracy (Succi Reference Succi2001; Krüger et al. Reference Krüger, Varnik and Raabe2009; Luo et al. Reference Luo, Lallemand and d’Humières2011). The particle migration trajectories and final equilibrium positions exhibit negligible changes with increasing lattice resolution, with variations in the equilibrium position remaining below
$0.05\,\%$
relative to the equilibrium position obtained for larger radius
$a = 35$
shown in figure 11(a). These results indicate that a particle radius of
$a = 30$
lattice units is sufficient to accurately resolve both the flow field and particle dynamics, and this value is therefore adopted for all subsequent simulations.
B.3. Time-step independence
Time-step independence is examined by fixing the particle radius at
$a = 30$
, which also fixes the domain size, and varying the wall velocity
$V_w$
shown in figure 11(b). This procedure effectively varies the Mach number and the relaxation time corresponds to the zero-shear-rate viscosity
$\tau _{0}$
while keeping all other parameters unchanged. The transient migration dynamics and final equilibrium positions are found to be nearly identical across all tested cases. Based on these results,
$V_w = 0.1$
and
$\tau _{0} = 0.8$
are selected for all simulations. In addition to providing a favourable balance between accuracy and efficiency, choosing
$\tau _{0} \approx 0.8$
is known to improve the accuracy of the bounce-back boundary condition in enforcing the no-slip condition at wall (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017) and also provides balanced efficiency, stability and accuracy at a given lattice size of numerical simulations with the LBM (Zhao Reference Zhao2013).
B.4. Aspect-ratio independence
The influence of the computational domain aspect ratio is assessed by performing simulations with
$AR = L/H = 1$
,
$2$
and
$3$
, while keeping all physical and numerical parameters fixed at
${\textit{Re}}_c = 300$
,
$Cu = 100$
,
$\beta = 0.1$
,
$n = 0.25$
,
$\kappa = 0.1$
,
$a = 30$
,
$\tau _{0} = 0.8$
and
$V_w = 0.1$
, as shown in figure 11(c). In all cases, the span-to-height ratio is maintained at
$W/H = 2$
. The streamwise domain length
$L$
is varied to achieve different
$AR$
, while the channel height is fixed according
$H=a/\kappa$
. Consequently, the number of lattice nodes in the streamwise direction (
$N_x$
) varies with
$L$
, whereas the number of nodes in the wall-normal (
$N_z$
) and spanwise (
$N_y$
) directions remains fixed, consistent with the fixed values of
$H$
and
$W$
, respectively. Particle trajectories, migration dynamics and final equilibrium positions are found to be indistinguishable for
$AR = 2$
and
$3$
, while a slight deviation in the equilibrium position is observed for
$AR = 1$
, as highlighted in the inset of figure 11(c). In particular, the equilibrium position for
$AR = 1$
deviates by approximately
$6.5\,\%$
relative to the reference case, indicating the presence of finite-domain effects at this aspect ratio. This suggests that the domain corresponding to
$AR = 1$
is marginally insufficient to fully eliminate finite-domain effects, whereas
$AR = 2$
provides an adequate representation of an effectively unbounded flow in the streamwise direction. It is also worth noting that, although the channel Reynolds number is
${\textit{Re}}_c = 300$
, the corresponding particle Reynolds number remains moderate, with
${\textit{Re}}_{\!p} \approx 3$
for
$\kappa = 0.1$
,
${\textit{Re}}_{\!p} \approx 12$
for
$\kappa = 0.2$
and
${\textit{Re}}_{\!p} \approx 0.75$
for
$\kappa = 0.05$
. As a result, no pronounced turbulent wake forms around the particle even at the highest Reynolds number considered, and an aspect ratio of
$AR = 2$
is sufficient to accurately capture the particle-induced disturbance flow. This choice is consistent with values employed in previous studies of inertial migration of particles in Newtonian shear flows (Fox, Schneider & Khair Reference Fox, Schneider and Khair2020). At substantially higher channel Reynolds numbers, a larger streamwise domain extent may be required. For the Reynolds numbers considered in the present study, an aspect ratio of
$AR = 2$
is sufficient to accurately capture the particle-induced flow.
B.5. Summary of numerical robustness
The combined grid-resolution, time-step, Mach-number and aspect-ratio independence tests demonstrate that the numerical results reported in this work are insensitive to further refinement of the spatial discretisation, temporal resolution or domain size. Since all other simulations correspond to lower Reynolds numbers or weaker rheological effects than the reference case
${\textit{Re}}_c = 300$
, they necessarily operate at equal or smaller Mach numbers and are therefore expected to be at least as accurate. All simulations presented in this study are performed with fixed values
$a = 30$
and
$\tau _{0} = 0.8$
, and
$\tau (\dot {\gamma })$
will vary accordingly, thereby ensuring consistency and numerical robustness across the parameter space.
a
z0
L
x
W
y
H
z
x
Vw
ci
wi
Zeq/H
Rec
Rec
Recr
Recr
ϵ=(Rec−Recr)/Recr
κ=0.05
Rep<1
Zeq/H−0.5
ϵ
(Rec,Cu)
κ=0.2
β=0.1
n=0.25
Recr
Rec
Rec
Zeq/H
Cu
Cu
Cucr
Cucr
Rec
Rec=100
κ=0.2
β=0.1
n=0.25
Zeq/H−0.5
δ=(Cu−Cucr)/Cucr
Rec
β=0.1
κ=0.2
n=0.25
Recr
Cu
β
κ=0.2
n=0.25
Recr
Cu
n
β=0.5
κ=0.2
κ=0.2
n=0.25
β=0.1
Rec=150
β
Rec=150
n=0.25
Cu=10
κ=0.2
n
κ=0.2
Rec=150
β=0.5
Cu=10
β=0.1
κ=0.2
n=0.25
Rec=150
Cu=0.1
Cu=1.0
Cu=10
Cu=100
Reℓ(x)=VwH/ν(x,z)
Cu=0.1
1
10
100
Rec=150
β=0.1
n=0.25
κ=0.2
Cu
Reℓ
n=0.25
κ=0.2
β=0.1
Rec=150
Cu=10
n=0.25
κ=0.2
Rec=150
Cu=10
β=0.5
κ=0.2
Rec=150
Cu
Z/H
Gt
Rec=50
β=0.5
n=0.25
κ=0.2
Rec
κ=0.05
0.1
0.2
β=0.5
Cu=10
n=0.25
Cu
κ=0.1
β=0.5
n=0.25
Recr
Cu
κ
n=0.25
β=0.5
U/Vw
Cu
Rec=150
β=0.1
n=0.25
κ
z
Ω/G
a
Vw
a=30
AR=1,2
3
a=30
Rec=300
Cu=100
β=0.1
n=0.25
κ=0.1
Open access
