1. Introduction
Precise manipulation of particle migration in microfluidic systems is fundamental to biomedical diagnostics, chemical analysis and lab-on-a-chip technologies (Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2018; Xia et al. Reference Xia, Wu, Zheng, Zhang and Wang2021; Battat, Weitz & Whitesides Reference Battat, Weitz and Whitesides2022; Xue et al. Reference Xue, Yin, Xu, Tian, Su and Hu2025). In Newtonian fluids, inertial effects make particles migrate to equilibrium positions in Poiseuille flow. Ho & Leal’s (Reference Ho and Leal1974) work on hydrodynamic interactions laid the groundwork for modern microfluidics, and pioneering work by Di Carlo et al. (Reference Di Carlo, Irimia, Tompkins and Toner2007, Reference Di Carlo, Edd, Humphry, Stone and Toner2009) initiated the inertial manipulation of particles in microfluidics. Classical theories of particle migration in Newtonian fluids, rooted in inertial lift forces, have been widely discussed and established (
$F_{i}=C_{i}\rho a^4\dot \gamma ^2$
, where
$C_{i}$
,
$\rho$
,
$a$
and
$\dot \gamma$
denote the inertial lift coefficient, fluid density, particle radius and shear rate, respectively) (Asmolov Reference Asmolov1999; Matas, Morris & Guazzelli Reference Matas, Morris and Guazzelli2004; Martel & Toner Reference Martel and Toner2014; Liu et al. Reference Liu, Xue, Sun and Hu2016; Asmolov et al. Reference Asmolov, Dubov, Nizkaya, Harting and Vinogradova2018). On the other hand, non-Newtonian fluids characterised by elastic stress responses introduce different migration mechanisms (Boyer, Pouliquen & Guazzelli Reference Boyer, Pouliquen and Guazzelli2011; Yang et al. Reference Yang, Kim, Lee, Lee and Kim2011; Lu et al. Reference Lu, Liu, Hu and Xuan2017; Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2018; Yuan et al. Reference Yuan, Zhao, Yan, Tang, Alici, Zhang and Li2018). These fluids enable particle migration driven by elastic forces, with theoretical frameworks predicting elastic lifts proportional to the first normal stress difference (Leshansky et al. Reference Leshansky, Bransky, Korin and Dinnar2007) competing with inertial lifts (
$\boldsymbol{F}_{e}=C_{e}a^3\boldsymbol{\nabla }N_1=C_{e}a^3\eta \lambda \boldsymbol{\nabla }\dot {\gamma }^2$
, where
$C_{e}$
,
$\eta$
,
$\lambda$
and
$N_1$
denote the elastic lift coefficient, solution viscosity, relaxation time and the first normal stress, respectively), which has been confirmed by computational works (Raffiee, Ardekani & Dabiri Reference Raffiee, Ardekani and Dabiri2019; Yu et al. Reference Yu, Wang, Lin and Hu2019; Yokoyama et al. Reference Yokoyama, Yamashita, Higashi, Miki, Itano and Sugihara-Seki2021). By taking advantage of this competition, it has been demonstrated that microparticles of different sizes can be separated (Liu et al. Reference Liu, Xue, Chen, Shan, Tian and Hu2015; Holzner, Stavrakis & DeMello Reference Holzner, Stavrakis and DeMello2017; Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2018).
A key phenomenon that has attracted attention in microfluidic systems recently is slip-induced migration, where particles exhibit cross-streamline motion due to the interplay between hydrodynamic forces the particle’s additive slip velocity (Li & Xuan Reference Li and Xuan2018, Reference Li and Xuan2023; Choudhary et al. Reference Choudhary, Li, Renganathan, Xuan and Pushpavanam2020; Lochab & Prakash Reference Lochab and Prakash2021; Cao et al. Reference Cao, Dvoriashyna, Liu, Lauga and Wu2022). Typically, the slip velocity originates from inherent properties of the particles, such as self-generated chemical forces (e.g. due to biological activity) or externally applied driving forces (e.g. electric fields or magnetic fields). These forces, acting on the particle, impart a velocity component that is parallel to the direction of the local flow. The non-fluid forces acting on the particles are balanced by the flow resistance of the particles. Consequently, the particles do not strictly follow the flow field; they may lead (as shown in figure 1
a) or lag relative to the local fluid velocity. The resulting velocity difference is defined as the slip velocity,
$u_s$
. For example, Cao et al. (Reference Cao, Dvoriashyna, Liu, Lauga and Wu2022) observed the phenomenon of Escherichia coli migrating to the tube wall in DNA solutions when swimming against flow direction. Li & Xuan (Reference Li and Xuan2018, Reference Li and Xuan2023) observed the migration of electrophoretic particles in polyethyleneoxide (PEO) solution, and they showed that the leading and lagging particles move towards the channel centre and the channel wall in viscoelastic flow, respectively, whereas the opposite behaviour is observed in a Newtonian fluid. Recently, this method has been demonstrated to manipulate nanoscale particles (Qi, Ma & Hu Reference Qi, Ma and Hu2025). Differences also exist in the aggregation and separation effects of particles in PEO solutions of different chain lengths (Choudhary et al. Reference Choudhary, Li, Renganathan, Xuan and Pushpavanam2020; Li & Xuan Reference Li and Xuan2023). Recent theoretical and computational studies demonstrate that the slip velocity modifies these forces, leading to migration reversal. Einarsson & Mehlig (Reference Einarsson and Mehlig2017) and Khair & Kabarowski (Reference Khair and Kabarowski2020) showed that the slipping particles driven by an external force are subjected to opposite slip-induced inertial and elastic lift in unbounded weakly inertial and weakly viscoelastic shear flows. They considered the slip particles in the infinite space shear flow and derived the analytical expression of the inertial and elastic lift induced by the slip. Zhang et al. (Reference Zhang, Murch, Einarsson and Shaqfeh2020) calculated the lift of slip particles in viscoelastic solutions using the immersed boundary method (IBM) and found that the imbalance of polymer stresses on both sides of the particles was responsible for the lift. Despite progress, critical gaps remain in our understanding of the mechanistic interplay between polymer conformation, the slip velocity and lifts. Besides, existing theoretical and computational fluid dynamics (CFD) models often neglect the microscale polymer dynamics and its impact on macroscopic flow of the particle, which makes it difficult to design the required polymer solution.
(a) Schematic representation of the slip velocity-induced lateral inertial and elastic lift. Taking the particle with a leading velocity (
$u_s\gt 0$
) as an example, the particle is subjected to an additional inertial lift directed towards the wall and an additional elastic lift directed towards the centre of the channel. (b) Schematic representation of the computational model. (c) The ensemble averages of the first normal stress differences,
$\langle \tau _{xx}-\tau _{yy}\rangle$
, during the flow in the microchannel at different
$\Delta x$
. (d–e) Comparison of the results between this work and conventional CFD for the lateral lift in Newtonian (d) and viscoelastic fluids (e).
The aim of this study is to develop a mesoscale model for investigating how the slip velocity alters polymer stress asymmetry and drives migration reversal in viscoelastic flows, while quantifying the competing roles of inertial and elastic forces under different slip velocity, particle size and polymer molecular weight conditions. Theoretically, we develop a theory of slip-induced elastic lift in dilute polymer solutions and semi-quantitatively predict the slip-induced elastic lift of particles. Methodologically, we employ a hybrid lattice Boltzmann method (LBM) and coarse-grained molecular dynamics (CGMD) framework to address the multiscale interactions between polymer chains and particles. Unlike conventional CFD methods, this approach helps to tune the molecular characteristics of the polymer to capture the stretching and stress relaxation of the polymer from a mesoscale perspective, which is crucial for predicting and designing migration in weak viscoelastic environments.
2. Problem formulation, computational model and experiments
2.1. Problem formulation
The problem of fluid in a microchannel is governed by the Navier–Stokes equations
\begin{equation} \rho \left(\frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{u}\right)=\boldsymbol{\nabla }\boldsymbol{\cdot }(-p\boldsymbol{I}+\boldsymbol{\tau }_v+\boldsymbol{\tau }_e)+\boldsymbol{R},\quad \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}
where
$p$
,
$\boldsymbol{\tau }_v$
,
$\boldsymbol{\tau }_e$
and
$\boldsymbol{u}$
denote the pressure, viscous stress, elastic stress and velocity field of the flow, and
$\boldsymbol{R}$
is the external force term, respectively. In the microchannel, no-slip and non-penetration boundary conditions exist at the stationary walls of the channel. The undisturbed velocity field at the small Reynolds number is denoted as
$\boldsymbol{u}_{\!f}$
. Consider the flow of a uniformly rigid medium particle with neutral buoyancy in a microchannel. The particle moves with the velocity
$\boldsymbol{u}_0$
at the centre-of-mass (COM),
$\boldsymbol{x}_{\!p}$
, and rotates with the angular velocity
$\boldsymbol{\omega }$
in the absence of an external force field. With an additional external force field exists, the equilibrium of the external force and the hydrodynamic force causes the steady velocity of the particle to change, resulting in an additional slip velocity
$\boldsymbol{u}_s$
. The disturbed velocity field is defined as
$\boldsymbol{u}_d=\boldsymbol{u}-\boldsymbol{u}_{\!f}$
. Furthermore, Morrison Jr (Reference Morrison1970) demonstrated that the disturbed velocity field
$\boldsymbol{u}_d$
surrounding a forced particle decays as
$1/\|\boldsymbol{x}-\boldsymbol{x}_{\!p}\|$
. Cox & Brenner (Reference Cox and Brenner1968) and Ho & Leal (Reference Ho and Leal1976) established that, for small particles within weakly inertial and viscoelastic fluid regimes, the disturbed velocity
$\boldsymbol{u}_d$
exhibits a vanishing behaviour in Poiseuille flow at large distances from the particle, i.e.
$\boldsymbol{u}_d \rightarrow 0$
as
$\|\boldsymbol{x}-\boldsymbol{x}_{\!p}\|\rightarrow \infty$
. Therefore, neglecting inertioelastic coupling, the velocity field satisfies
\begin{equation} \boldsymbol{u}=\boldsymbol{u}_0+\boldsymbol{\omega }\times (\boldsymbol{x}-\boldsymbol{x}_{\!p})+\boldsymbol{u}_s \quad \mbox{as }\quad \|\boldsymbol{x}-\boldsymbol{x}_{\!p}\|\leqslant a, \end{equation}
\begin{equation} \boldsymbol{u} \rightarrow \boldsymbol{u}_{\!f} \quad \mbox{as }\quad \|\boldsymbol{x}-\boldsymbol{x}_{\!p}\|\rightarrow \infty. \end{equation}
In accordance with Einarsson & Mehlig (Reference Einarsson and Mehlig2017), the effects of inertioelastic coupling can be neglected when the particle Reynolds number is much smaller than the Deborah number (
${\textit{Re}}_p=\rho u_ca/\eta \ll {\textit{De}}=\lambda u_c/a$
, where
$u_c$
is the characteristic velocity of the particle), which is equivalent to requiring that
$\rho a^2/(\eta \lambda )\ll 1$
.
For small particles under the Stokes flow, we assume that the local shear rate around the particles is approximately equal to the unperturbed flow field and that the weakly viscoelastic flow satisfies the condition in which the second normal stress coefficient is much smaller than the first. In a microchannel, the transverse forces acting on a particle consist of inertial and elastic lifts, as well as additional slip-induced contributions
\begin{equation} \boldsymbol{F}_{t}=\boldsymbol{F}_{\! i}+\boldsymbol{F}_{e}+\boldsymbol{F}_{\! i,s}+\boldsymbol{F}_{\! e,s}. \end{equation}
According to Einarsson & Mehlig (Reference Einarsson and Mehlig2017) and Khair & Kabarowski (Reference Khair and Kabarowski2020), the slip-induced inertial and elastic lifts are perpendicular to
$\boldsymbol{u}_s$
, and also perpendicular to the streamlines, and can be expressed as
\begin{equation} F_{i,s}=C_{i,s}a\eta u_s\sqrt {{\textit{Re}}_l}, \end{equation}
\begin{equation} \boldsymbol{F}_{\!e,s}=C_{e,s}a\eta {{Wi}_l} \hat {\varOmega }\times \boldsymbol{u}_s + O({Wi}_l^2), \end{equation}
where
$C$
, and
$\hat {\varOmega }$
are the coefficient and the unit vector along the vorticity direction, respectively. Here,
${\textit{Re}}_l=\rho \dot {\gamma }a^2/\eta$
and
${Wi}_l=\lambda \dot {\gamma }$
are the local Reynolds and Weissenberg numbers without slip velocity disturbance. Einarsson & Mehlig (Reference Einarsson and Mehlig2017) demonstrated that the second term
$O({Wi}_l^2)$
is perpendicular to the first term for
$\boldsymbol{u}_s$
in any three-dimensional direction. In particular, the second term
$O({Wi}_l^2)=0$
vanishes when
$\boldsymbol{u}_s$
is parallel to the mainstream direction according to their calculation. Consequently, truncating (2.6) at
$O({Wi}_l)$
yields
\begin{equation} \boldsymbol{F}_{\!e,s} =C_{e,s}a\eta {Wi}_l \hat {\varOmega }\times \boldsymbol{u}_s. \end{equation}
For the dilute solution of the ideal chains, the Zimm model shows the following relationships (Zimm Reference Zimm1956; Briels & Padding Reference Briels and Padding2007; Hiemenz & Lodge Reference Hiemenz and Lodge2007):
\begin{equation} \lambda =\frac {\eta _sb^3M^{3/2}}{\sqrt {3\pi }k_{{B}}T}, \end{equation}
\begin{equation} \eta (\dot {\gamma })=\eta _s+\frac {cRT\lambda }{Mm}\sum _{q=1}^M \frac {q^{-3/2}}{1+\dot {\gamma }^2 \lambda ^2q^{-3}}, \end{equation}
where
$\eta _s$
,
$b$
,
$M$
,
$k_{\textrm {B}}$
,
$T$
,
$c$
,
$R$
,
$m$
and
$q$
are solvent viscosity, Kuhn length of each monomer, number of monomer links of each polymer, Boltzmann’s constant, temperature, concentration, universal gas constant and molar mass of each monomer and the index of summation, respectively.
With the substitution of (2.8) and (2.9) into (2.7), one can obtain
\begin{align} \begin{split} \boldsymbol{F}_{\!e,s}=C_{e,s}\frac {\eta _s^2b^3M^{3/2}\dot \gamma }{\sqrt {3\pi }k_{ { B}}T}\left(1+\frac {cRb^3M^{1/2}\sqrt {3\pi }k_{{B}}T^2}{m}\sum _{q=1}^{M}\frac {q^{-3/2}}{3\pi k^2_{ { B}}T^2+\dot \gamma ^2\eta _s^2b^6M^3q^{-3}}\right)a\hat {\varOmega }\times \boldsymbol{u}_s. \end{split} \end{align}
2.2. Computational model
On the mesoscopic scale, viscoelastic fluids are obtained by dissolving elastic polymer chains in Newtonian fluids. We apply an IBM that combines the LBM and CGMD to model viscoelastic fluids and particles (Mackay & Denniston Reference Mackay and Denniston2013; Mackay, Ollila & Denniston Reference Mackay, Ollila and Denniston2013; Denniston et al. Reference Denniston, Afrasiabian, Cole-André, Mackay, Ollila and Whitehead2022), as shown in figure 1(b).
Essentially, our LBM-CGMD-IBM approach solves the governing equations of viscoelastic fluids at the continuum scale (2.1) by explicitly resolving the polymer dynamics and incorporating its contributions into the mesoscopic LBM framework. The LBM provides a numerical solution to the Navier–Stokes equations in the hydrodynamic limit (Chen & Doolen Reference Chen and Doolen1998; Adhikari et al. Reference Adhikari, Stratford, Cates and Wagner2005) by computing solvent density, velocity, pressure and viscous stress (
$\boldsymbol{\tau }_v$
) while inherently enforcing continuity (
$\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}=0$
). The elastic stress tensor (
$\boldsymbol{\tau }_e$
) originates directly from the polymer chains modelled via CGMD.
The IBM couples the CGMD and LBM solvers by interpolating velocities and distributing forces, thereby imposing polymer-induced stresses on fluid motion. External forces acting on the particle beads are incorporated into the LBM as part of the momentum source term,
$\boldsymbol{R}$
, completing the two-way coupling between the solute and the solvent. This strategy enables the combined solver to capture the response of the Newtonian solvent and the non-Newtonian elastic stresses of the polymer solution within the same framework. Similar hybrid methods have been successfully applied to polymer simulations (Ahlrichs & Dünweg Reference Ahlrichs and Dünweg1999; Ding et al. Reference Ding, Duan, Lu and Shi2015; Shi et al. Reference Shi, Wu, Chen, Xu, Yang and Ding2024) and blood flow (Ye, Shen & Li Reference Ye, Shen and Li2018, Reference Ye, Shen and Li2019a
, Reference Ye, Shen, Wei and Li2021).
Specifically, the LBM computes the fluid density, velocity and pressure, as well as the solvent viscous stress term (
$\tau _v$
). It also enforces the continuity equation (
$\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}=0$
) through its structural properties. The elastic stress tensor (
$\tau _e$
) in (2.1) arises directly from the polymer chains modelled in the CGMD framework. In this framework, the chains are constructed from beads connected by harmonic bonds. Rigid particles are also represented as bead assemblies within this framework. The IBM couples the dynamics of polymers and particles to the solvent by interpolating velocities and distributing forces to ensure two-way coupling. This mechanism enables the macroscopic effects of polymer-induced elastic stresses to be applied to fluid motion simulations using the LBM. Furthermore, the LBM includes external forces acting on particle beads in CGMD as contributions to the momentum source term,
$\boldsymbol{R}$
, thereby closing the loop of multi-scale coupling.
2.2.1. Lattice Boltzmann method for fluid flow
The LBM simulates incompressible, Newtonian flow on a rectangular mesh by discretising terms in the linearised Boltzmann equation
\begin{equation} (\partial _t+e_{kl}\partial _l)f_k=-\frac {1}{\tau }(f_k-f_k^{eq})-F_{l}\partial _{p_l}f_k, \end{equation}
where
$f_k$
,
$f_k^{eq}$
,
$p_l$
and
$\tau$
represent the distribution function corresponding to a discrete velocity, the lattice equilibria, a component of the momentum density and the non-dimensional relaxation time, respectively. The external force density acting on the fluids is represented by
$F_{l}=\partial p_l/\partial t$
.
By introducing auxiliary distribution functions
\begin{equation} \bar {f}_k(\boldsymbol{x},t)=f_k(\boldsymbol{x},t)-\frac {\Delta t}{2}R_k(\boldsymbol{x},t), \end{equation}
a finite difference scheme is derived
\begin{equation} \bar {f}_k(\boldsymbol{x}+\boldsymbol{e_k}\Delta t,t+\Delta t/2)=\bar {f}_k(\boldsymbol{x},t)+\frac {\Delta t}{\tau +\Delta t/2}\big[-(\bar {f}_k-f_k^{eq})-\tau F_{l}\partial _{p_l}f_k^{eq}\big], \end{equation}
where
$\Delta t$
is the time step, and
$R_k(\boldsymbol{x},t)$
is the right-hand side of (2.11). The local equilibrium distribution function is given by
\begin{equation} f_k^{eq} = \rho w_{c_k}\! \left [ 1 + \frac {\boldsymbol{c}_k \boldsymbol{\cdot }\boldsymbol{u}}{c_s^2} + \frac {(\boldsymbol{c}_k \boldsymbol{\cdot }\boldsymbol{u})^2}{2c_s^4} - \frac {u^2}{2c_s^2} \right ]\!, \end{equation}
with
$ c_s = (\Delta x/\Delta t)/\sqrt {3}$
representing the sound speed, and
$\Delta x$
denoting the lattice spacing. The weight factors
$w_{c_k}$
depend on the magnitude of
$\boldsymbol{c}_k$
. In our simulations, we use the three-dimensional, 19-velocity model (D3Q19), where the set of
$\boldsymbol{c}_k$
includes the rest velocity (
$k=0$
), the 6 nearest neighbours (
$k=1{-}6$
) and the 12 next-nearest neighbours (
$k=7{-}18$
) on a simple cubic lattice. The corresponding weight factors for the D3Q19 model are
$ w_0 = 1/3$
,
$w_k = 1/18$
for
$(k=1{-}6)$
, and
$ w_k = 1/36$
for
$(k=7{-}18)$
.
The external forcing term is discretised by (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002)
\begin{equation} F_{l,k} = \left ( 1 - \frac {1}{2\tau } \right ) w_k\! \left [ \frac {\boldsymbol{e}_k - \boldsymbol{u}}{c_s^2} + \frac {(\boldsymbol{e}_k \boldsymbol{\cdot }\boldsymbol{u})}{c_s^4} \boldsymbol{e}_k \right ] \boldsymbol{\cdot }\boldsymbol{F}_l. \end{equation}
Once the particle density distributions are known throughout the fluid domain, fluid properties such as solvent density, shear viscosity and fluid velocity can be calculated as
\begin{equation} \rho _s=\sum _k \bar {f}_k, \end{equation}
\begin{equation} \eta _s=\rho _sc_s^2\! \left(\tau -\frac {\Delta t}{2}\right)\!, \end{equation}
\begin{equation} \boldsymbol{u}=\frac {1}{\rho _s}\! \left(\sum _k \bar {f}_k\boldsymbol{e}_k+\frac {\Delta t}{2}\boldsymbol{F}\right)\!. \end{equation}
The half-step fluid velocity is defined as
\begin{equation} \boldsymbol{\bar {u}}=\frac {1}{\rho _s}\sum _k \bar {f}_k\boldsymbol{e}_k, \end{equation}
which is crucial for constructing the particle–fluid coupling force.
2.2.2. Coarse-grained molecular dynamics for polymers and particle
The polymer solution are obtained by immersing coarse-grained (CG) chains and particles in to the Newtonian fluid. Every two adjacent beads in the polymer chain are connected with a harmonic bond
\begin{equation} V_{b}=k(r-r_0)^2, \end{equation}
where
$k$
is the spring strength and
$r_0={0.4}\,{\unicode{x03BC} }\textrm {m}$
is the equilibrium length of CG beads. Each chain consists of
$M+1$
beads, corresponding to the polymer length,
$L_p=Mr_0$
. And a weak repulsion between beads maintains the uniformity of the system, represented by
\begin{equation} \boldsymbol{F}^a_{\textit{ij}}=a_{\textit{ij}}\boldsymbol{\hat {r}}_{\textit{ij}}, r_{\textit{ij}}\lt r_c, \end{equation}
where
$a_{\textit{ij}}$
,
$r_c$
,
$r_{\textit{ij}}$
,
$\boldsymbol{\hat {r}}_{\textit{ij}}$
are the conservative force coefficient, cutoff radius, the distance and the unit vector between beads. Here,
$a_{\textit{ij}}=75k_{ { B}}T/([L]^3\rho _N)=2.5k_{{B}}T$
(Groot & Warren Reference Groot and Warren1997; Ye et al. Reference Ye, Pan, Huang and Liu2019b
), where
$[L]={1}\,\rm {\unicode{x03BC} m}$
is the length unit of the simulation system in this work. The rigid particle is modelled by
$N$
beads and has a mass of
$m_p$
and a COM of
$\boldsymbol{X}_{\!p}=(X_{\!p},Y_{\!p},Z_{\!p})$
.
2.2.3. Immersed boundary method for fluid–structure interaction
The CG beads and the LBM fluid are coupled through the IBM. This method ensures the no-slip condition between the beads and the fluid by interpolating the velocity and force distributions between the Lagrangian and Eulerian mesh points (Mittal & Iaccarino Reference Mittal and Iaccarino2005; Verzicco Reference Verzicco2023). The CG bead
$X$
, with Lagrangian coordinate
$\boldsymbol{n}$
, moves at the same velocity as the surrounding fluid. This velocity can be interpolated from the fluid velocity using a smoothed Dirac delta function
\begin{equation} \frac {\partial \boldsymbol{X}(\boldsymbol{n},t)}{\partial t}=\boldsymbol{U}(\boldsymbol{X}(\boldsymbol{n},t))=\int _\varPi \boldsymbol{u}(\boldsymbol{x},t)\delta (\boldsymbol{x}-\boldsymbol{X}(\boldsymbol{n},t)){\rm d}\varPi, \end{equation}
where
$\varPi$
represents the simulation domain. This condition governs the movement and deformation of the CG polymer and the particle beads. The force density on the CG beads is distributed to the surrounding fluid mesh points, satisfying
\begin{equation} \boldsymbol{f}(\boldsymbol{x},t)=\int _\varPi \boldsymbol{F}(\boldsymbol{n},t)\delta (\boldsymbol{x}-\boldsymbol{X}(\boldsymbol{n},t)){\rm d}\varPi. \end{equation}
The coupling between the CG beads and the fluid is achieved using a constraint-based algorithm that ensures the velocity of each bead matches the interpolated velocity of the fluid at the bead’s position at the end of each time step. Specifically, the fluid velocity is interpolated to the off-lattice position of a bead
$i$
using a mass-weighted average to ensure momentum conservation. The interpolated half-step fluid velocity
$\bar {\boldsymbol {u}}_i$
is calculated as
\begin{equation} \bar {\boldsymbol{u}}_i = \frac {\sum _k \rho (\boldsymbol{x}_k) \bar {\boldsymbol{u}}(\boldsymbol{x}_k) \xi _{\textit{ik}}}{\sum _k \rho (\boldsymbol{x}_k) \xi _{\textit{ik}}}, \end{equation}
where
$\xi _{\textit{ik}}=\phi (\Delta x_{\textit{ik}})\phi (\Delta y_{\textit{ik}})\phi (\Delta z_{\textit{ik}})$
is the interpolation weight from fluid lattice site
$k$
to bead
$i$
. In this study, the 4-point stencil is used for the interpolation function
\begin{equation} \phi (x) = \left \{ \begin{array}{ll} \dfrac {3}{2}|x|^3-\dfrac {5}{2}|x|^2+1, & |x|\lt 1, \\ [5pt] -\dfrac {1}{2}|x|^3+\dfrac {5}{2}|x|^2-4|x|+2, & 1\lt |x|\lt 2.\end{array} \right . \end{equation}
The forces on each CG beads are derived from
\begin{equation} \boldsymbol{F}=-\boldsymbol{\nabla }V-\boldsymbol{F}^H=\boldsymbol{F}^M-\boldsymbol{F}^H=-\boldsymbol{\nabla }V_b+\boldsymbol{F}^a-\boldsymbol{F}^H, \end{equation}
where
$\boldsymbol{F}^H$
is a constraining force applied to the CG bead to enforce velocity matching, and
$\boldsymbol{F}^M$
is the molecular dynamics force. The constraining force is derived from the constraint condition and is given by
\begin{equation} \boldsymbol{F}^H = \varGamma \left ( \bar {\boldsymbol{u}} - \boldsymbol{U}(t+\Delta t/2) \right ), \end{equation}
where
$\varGamma$
is the force coupling constant, calculated according to Denniston et al. (Reference Denniston, Afrasiabian, Cole-André, Mackay, Ollila and Whitehead2022)
\begin{equation} \boldsymbol{\varGamma }=\frac {2m_i m_n^i}{(m_i+m_n^i)\Delta t}, \end{equation}
where
$m_i$
is the mass of the CG bead, and
$m_n^i$
is a representative fluid mass at the bead location, given by
\begin{equation} m_n^i = \frac { \big ( \sum _j \zeta _{\textit{ik}} \big )^2\sum _k \rho _k \xi _{\textit{ik}} \Delta x^3}{ \sum _j \xi _{\textit{ik}} \zeta _{\textit{ik}} }, \end{equation}
where
$\zeta _{\textit{ik}}$
is the spreading weight
\begin{equation} \zeta _{\textit{ik}} = \frac {\xi _{\textit{ik}} |\xi _{\textit{ik}}|}{\sum _s |\xi _{sk}|}. \end{equation}
The motion of the CG beads is solved by the velocity-Verlet algorithm (Verlet Reference Verlet1967; Swope et al. Reference Swope, Andersen, Berens and Wilson1982) with a time step
$\Delta t={0.02}\,{\unicode{x03BC} }\textrm {s}$
. Specifically, the bead velocity at the end of the integration step is computed from the half-step velocity
\begin{equation} \boldsymbol{U}_{\!i}(t+\Delta t)=\boldsymbol{U}_{\!i}(t+\Delta t/2)+\frac {\Delta t}{2m_i}\big(\boldsymbol{F}^M_i-\boldsymbol{F}^H_i\big), \end{equation}
and the velocity of fluid site
$k$
is updated by
\begin{equation} \boldsymbol{u}_k= \bar {\boldsymbol{u}}_k+\frac {\Delta t}{2\rho _{s,k}\Delta x^3}\sum _i\left(\boldsymbol{F}_{\!i}^H+\frac {m_n^i}{m_i+m_n^i}\boldsymbol{F}^M_i\right)\frac {\zeta _{\textit{ik}}}{\sum _k\zeta _{\textit{ik}}}. \end{equation}
2.2.4. Model set-up and validation
The LBM and CGMD are both solved in LAMMPS (Thompson et al. Reference Thompson2022). Unless otherwise stated, we set
$T={300}\,\textrm {K}$
,
$\rho _s={1000}\,{\rm {kg\,m}^{-3}}$
,
$\eta _s={1.0}\,{\textrm {mPa}\boldsymbol\, \textrm {s}}$
,
$k={0.6}\,\textrm {nJ}$
,
$m={0.0003}\,\textrm {pg}$
,
$M=19$
, block ratio
$\kappa =2a/W=0.1$
, the number density
$\rho _N={30}\,{\unicode{x03BC} }{\textrm {m}^{-3}}$
, the density of viscoelastic fluid
$\rho =m\rho _N+\rho _s={1009}\,\rm {kg\,m}^{-3}$
and the concentration
$c=m\rho _N/(m\rho _N+\rho _s)\approx 0.89\,\%$
, which is consistent with common PEO dilute solutions.
Based on previous simulations (Schlijper, Hoogerbrugge & Manke Reference Schlijper, Hoogerbrugge and Manke1995; Kong et al. Reference Kong, Manke, Madden and Schlijper1997; Salerno et al. Reference Salerno, Agrawal, Perahia and Grest2016; Asano, Watanabe & Noguchi Reference Asano, Watanabe and Noguchi2018), we set
$M=O(10)$
, in which each CG bead captures the collective behaviour of multiple monomers. This approach simplifies the polymer representation while retaining essential physical characteristics. Our work focuses on the collective polymer behaviour and its influence on particle lift forces rather than on resolving the fine-scale dynamics of individual monomers. Employing this CG method allows us to balance computational efficiency and physical fidelity. We use an Eulerian fluid mesh resolution
$\Delta x=1.25r_0={0.5}\,{\unicode{x03BC} }\textrm {m}$
, in which case the corresponding elastic number is
${El}={Wi}/{\textit{Re}}=0.78$
. In our simulations,
$\textit{Re}=\rho u_m W/\eta$
and
${Wi}=2\lambda u_m/W$
are the Reynolds and Weissenberg numbers of the flow, where the apparent viscosity
$\eta$
is determined by a Poiseuille flow test and
$u_m$
is the maximum value of the flow velocity. The calculated relation time
$\lambda =(\tau _{xx}-\tau _{yy})/(2\dot {\gamma }\tau _{xy})\approx {0.46}\,\textrm {ms}$
from shear flow tests agrees well with the measured value of
${0.36}\,\textrm {ms}$
for 1000 kDa PEO in experiments (Li & Xuan Reference Li and Xuan2023). For a typical PEO monomer with a length of
${{\sim} 0.58}\,\textrm {nm}$
and a molecular weight of
${{\sim} 44}\,\textrm {Da}$
(Hammouda & Ho Reference Hammouda and Ho2007), the contour length of a single
${600}\,\textrm {kDa}$
-chain is approximately
${7.9}\,{\unicode{x03BC} }\textrm {m}$
. In a later section of this study, values of
$M$
ranging from
$9$
to
$24$
correspond to polymer chain lengths from
$3.6$
to
${9.6}\,{\unicode{x03BC} }\textrm {m}$
, which is consistent with the calculated contour length. The CG chains are fully flexible because no bond-angle constraints are imposed. The average radius of gyration in our simulations,
$\langle R_g \rangle = \sqrt {\sum m_i r_i^2/\sum m_i} =0.31 {-}{0.59}\,{\unicode{x03BC} }\textrm {m}$
, agrees well with the reported experimental values of
$ 0.4{-}{0.8}\,{\unicode{x03BC} }\textrm {m}$
for PEO clusters (Devanand & Selser Reference Devanand and Selser1991).
The size of the channel is set to
$L=H=W={50}\,{\unicode{x03BC} }\textrm {m}$
,
$D=\sqrt {H^2+W^2}$
with periodic boundary conditions applied at both the inlet and outlet of the computational domain, allowing us to ignore the periodic effects on flow. The slip velocity is introduced by applying an additional force along the flow direction of the particles. In this work, it is important to clarify that, while particles naturally experience slip due to their finite size and interactions with the fluid, the slip velocity we discuss is primarily imposed by external forces. To determine the lift force acting on a rigid particle at a given lateral position (
$2y/H$
,
$2z/H$
), the lateral velocities are set to zero. Under these conditions, the particle, initially placed at the specified coordinates, begins to translate and rotate due to the forces exerted by the surrounding fluid. The system reaches equilibrium when the translational velocity in the x-direction and the angular velocities in the y- and z-directions stabilise. This occurs after approximately
$150\,000$
time steps. At equilibrium, the hydrodynamic force at various lateral positions can be calculated by evaluating the lateral forces exerted by the fluid. The lateral lifts acting on the rigid particle in Newtonian (inertial lift force,
$\boldsymbol{F}_{\!i}$
) and viscoelastic fluids (total lift force,
$\boldsymbol{F}_{t}$
) are obtained by ensemble averaging of the corresponding forces for over
$300\,000$
time steps, during which the particle moves a distance of approximately
$600$
to
${800}\,{\unicode{x03BC} }\textrm {m}$
. We assume that the inertial lift force (
$\boldsymbol{F}_{\!i}$
) experienced by the particles in the viscoelastic fluid is consistent with that in a Newtonian fluid, provided that
$F_t\gt \gt F_i$
. In this work, the typical values are
$\rho a^2/(\eta \lambda )\approx 0.046{-}0.183$
. Consequently, inertioelastic coupling is neglected, and the lift forces are obtained by linear superposition, as shown in (2.4). The elastic lift force is defined as
$\boldsymbol{F}_{e}=\boldsymbol{F}_{t}-\boldsymbol{F}_{\!i}$
. Additionally, we can get a new total lift and inertial lift when the slip velocity
$u_s\neq 0$
. Specifically, the total lift is given by
$\boldsymbol{F}_{\!t,s}=\boldsymbol{F}_{t}(u_s\neq 0)-\boldsymbol{F}_{t}(u_s=0)$
, and the inertial lift is given by
$\boldsymbol{F}_{\!i,s}=\boldsymbol{F}_{\!i}(u_s\neq 0)-\boldsymbol{F}_{\!i}(u_s=0)$
. Thus, the elastic lift induced by slip can be expressed as
$\boldsymbol{F}_{\!e,s}=\boldsymbol{F}_{\!t,s}-\boldsymbol{F}_{\!i,s}$
. This is similar to the method of calculating the lift force of particles using the conventional CFD method (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Chen et al. Reference Chen, Xue, Zhang, Hu, Jiang and Sun2014; Liu et al. Reference Liu, Xue, Sun and Hu2016; Raffiee et al. Reference Raffiee, Ardekani and Dabiri2019). It takes about 552 CPU core hours using AMD EPYC 9,654 2.4 GHz 192-core processors for a typical simulation case.
Grid independence studies of the flow were conducted by comparing the ensemble averages of the first normal stress differences,
$\langle \tau _{xx}-\tau _{yy}\rangle$
, during the flow in the microchannel at
$\Delta x = 0.25, 0.5, 0.75, {1}\,{\unicode{x03BC} }\textrm {m}$
. The results shown in figure 1(c) indicate that
$\langle \tau _{xx}-\tau _{yy}\rangle$
differs slightly at different resolutions but gradually converges as
$\Delta x$
decreases. The difference between
$\Delta x=0.25$
and
${0.5}\,{\unicode{x03BC} }\textrm {m}$
is negligible. Considering computational efficiency, we adopted
$\Delta x=1.25r_0={0.5}\,{\unicode{x03BC} }\textrm {m}$
in this study, which is also consistent with previously reported suitable resolutions where
$r_0=(0.6{-}0.8)\Delta x$
(MacMeccan et al. Reference MacMeccan, Clausen, Neitzel and Aidun2009; Ye et al. Reference Ye, Shen and Li2019a
). In addition, the present settings effectively capture the inertia and elastic lift of the particles. The inertial lift coefficient
$C_i=F_iW^2/(\rho u_m^2a^4)$
and the normalised total lift in viscoelastic fluids calculated in this work are consistent with the results of conventional CFD (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Raffiee et al. Reference Raffiee, Ardekani and Dabiri2019), as shown in figures 1(d) and 1(e).
Previous studies have shown that, when
$Wi\gt 0$
, elastic lift dominates and the particles reach equilibrium along the diagonal since
$\boldsymbol{F}_{e}\sim \boldsymbol{\nabla }N_1\sim \boldsymbol{\nabla }\dot \gamma ^2$
(Yang et al. Reference Yang, Kim, Lee, Lee and Kim2011; Stoecklein & Di Carlo Reference Stoecklein and Di Carlo2018; Yokoyama et al. Reference Yokoyama, Yamashita, Higashi, Miki, Itano and Sugihara-Seki2021). Figure 2(a) shows the normalised shear rate, indicating that the shear rate gradient is most pronounced along the diagonal. Hence, we calculate the lift acting on the particle with
$\kappa =0.1$
across the diagonal (
$\alpha$
direction) under
${\textit{Re}}=2.09$
and
${\textit{Re}}=13.4$
, and
$\kappa =0.3$
under
${\textit{Re}}=13.4$
, as shown in figure 2(b–d). In general, the positive
$F_i$
near the centre of the channel pushes the particle away from the centre, while it becomes negative near the wall. The elastic lift
$F_e$
exhibits a trend opposite to the inertial lift
$F_i$
. For small particles, at a low
${\textit{Re}}$
, the effect of
$F_i$
is relatively small, and
$F_t$
is similar to that of
$F_e$
. As a result, particles migrate towards the centre and the corners. As
${\textit{Re}}$
increases, the effect of the wall lift increases and the
$F_t$
value of the particles is negative throughout the computational domain, causing the particles to migrate towards the centre. For large particles,
$F_i$
and
$F_e$
are of the same order of magnitude and the competition causes
$F_t$
to be positive near the centre and negative near the channel wall. This is because that
$F_i$
is proportional to the fourth power of
$a$
, and as
$\kappa$
increases it will rapidly exceed
$F_e$
, which is proportional to the third power of the particle size. The stable position of the large particles will be at a certain position between the centre of the tube and the wall of the channel. These results are supported by previous computational (Li, McKinley & Ardekani Reference Li, McKinley and Ardekani2015; Raffiee et al. Reference Raffiee, Ardekani and Dabiri2019; Yu et al. Reference Yu, Wang, Lin and Hu2019; Yokoyama et al. Reference Yokoyama, Yamashita, Higashi, Miki, Itano and Sugihara-Seki2021) and experimental phenomena (Yang et al. Reference Yang, Kim, Lee, Lee and Kim2011; Liu et al. Reference Liu, Xue, Chen, Shan, Tian and Hu2015).
(a) The dimensionless shear rate distribution of the channel cross-section. Darker colours represent areas of larger magnitude. (b–d) The dimensionless lateral total, inertial and elastic lift under different conditions.
2.3. Experiments
In our experiments, we use a
${2}\,\textrm {cm}$
-long rectangular polydimethylsiloxane microchannels with a cross-section of
$50\times {50}\,{\unicode{x03BC} }\textrm {m}^{2}$
. Suspensions contain
${5}\,{\unicode{x03BC} }\textrm {m}$
polystyrene particles (
${\lt} 0.1\,\%$
vol) with
$0.1\,\%$
Tween 20 for dispersion in the Newtonian (
${0.01}\,\textrm {mM}$
phosphate buffer) and viscoelastic fluids (buffer with
$0.2\,\%$
PEO, 600 kDa, Sigma–Aldrich). The fluids are infused at
${18}\,{\unicode{x03BC} }\rm {l\,h^{-1}}$
for Newtonian and
${90}\,{\unicode{x03BC} }\rm {l\,h^{-1}}$
for viscoelastic fluids via a syringe pump. A
${200}\,\rm {V\,cm^{-1}}$
DC electric field is applied using a function generator (AFG31000, Tektronix) connected to a high-voltage amplifier (ATA-7050, Aigtek). Experimental images are observed using a microscope (IX73, Olympus) and captured with a high-speed camera (VEO1010L, Phantom). The images are processed using the ImageJ software package, employing z-projection with the `standard deviation’ option.
3. Results
In pressure-driven flows, particles of finite size do not precisely follow the velocity of the undisturbed fluid. Even in the Stokes regime, Faxén corrections predict curvature-induced slip, and direct simulations confirm that neutrally buoyant spheres in Poiseuille flow ‘lag the fluid slightly’ (Feng, Hu & Joseph Reference Feng, Hu and Joseph1994). In contrast, the slip velocity considered in this study is imposed externally by forces that act along the flow direction and are not part of the fluid. This is achieved in simulations by applying a body force to the particles and in experiments by subjecting them to a DC electric field that drives electrophoretic migration in the opposite direction to the field. This imposed slip velocity is significantly larger than natural finite-size slip and thus induces more pronounced effects. This allows particles to lead (
$u_s\gt 0$
) or lag (
$u_s\lt 0$
) the local flow. The resulting velocity difference is defined as the slip velocity,
$u_s$
.
Figure 3(a) shows a top view of particles at the inlet and outlet of a microchannel for Newtonian and viscoelastic fluid flows driven by combined pressure and electric fields. The bright lines are composite images created by superimposing multiple instantaneous frames in order to visualise the behaviour of the particles being focused. Under the applied electric field, particles migrate electrophoretically in the opposite direction to the field. When the slip velocity is zero (
$u_s^*=u_s/u_m=0$
), no significant particle focusing occurs at the inlet or outlet for either fluid type. This indicates that inertial and elastic lift forces alone are insufficient to achieve focusing within a
${2}\,\textrm {cm}$
-long channel. With a positive slip velocity (
$u_s^*\gt 0$
), particles migrate toward the channel walls in Newtonian fluids but focus at the centreline in viscoelastic fluids. Conversely, when the slip velocity is negative (
$u_s^*\lt 0$
),the focusing behaviour is reversed. These observations confirm that slip-induced lift is the dominant factor governing particle migration.
(a) Composite images of slip-induced migration for
${5}\,{\unicode{x03BC} }\textrm {m}$
-diameter particles in Newtonian and viscoelastic fluids, where the slip of the particles is achieved by the application of a direct current electric field. The bright lines represent composite images of the particles formed from multiple instantaneous frames. When the slip velocity
$u_s=0$
, no significant particle focusing is observed at the inlet and outlet of the
${2}\,\textrm {cm}$
-long channel for both Newtonian and viscoelastic fluids, indicating that the small elastic and inertial lift forces are insufficient for particle focusing. However, with a leading slip velocity (
$u_s\gt 0$
), particles migrate towards the channel walls in Newtonian fluids, while they focus at the centre in viscoelastic fluids. Conversely, with a slip velocity (
$u_s\lt 0$
), the focusing behaviour is reversed between the two types of fluid. The dimensionless total lift (b), inertial lift (c), elastic lift (d) and slip-induced lift (e) of particles with
$\kappa =2a/W=0.1$
,
${\textit{Re}}=2.09$
,
${Wi}=1.64$
.
To investigate these functional dependencies, we compute the lateral dimensionless lift force,
$F^*=F/(W\eta _s u_m)$
, along the diagonal of the channel cross-section as a function of the dimensionless parameter
$2\alpha /D$
. This normalisation allows us for systematic comparison of how
$F_{e,s}^*$
varies with different parameters. Given that
$|u_s^*|=0.04273$
and
${Re=2.09}$
, the total, inertial, elastic and slip-induced lifts of the particle with
$\kappa =0.1$
in the viscoelastic solution are shown in figure 3(b–e). The value of
$F_t^*$
of the particles varies dramatically at small
$u_s^*$
. Specifically,
$F_t^*\lt 0$
over the entire computational domain for
$u_s^*\gt 0$
, resulting in the particles migrating towards the centre of channel. If
$u_s^*\lt 0$
, this causes
$F_t^*\gt 0$
, causing the particles to migrate to the corners. By separating the inertial and elastic parts of the lift, we find that
$F_i^*$
caused by
$u_s^*\gt 0$
shifts upwards and becomes positive, while
$u_s^*\lt 0$
causes
$F_i^*$
to shift downwards, leading to the migration of particles towards the wall and centre, respectively, in the Newtonian fluid. However, the effect of this tendency on
$F_e^*$
is reversed. We compare
$F_{i,s}^*$
and
$F_{e,s}^*$
induced by the same
$u_s^*$
, and, as shown in the figure 3(e),
$F_{e,s}^*$
and
$F_{i,s}^*$
are in a competitive relationship, with
$F_{e,s}^*$
being an order of magnitude higher. This leads to the dominant role of
$F_{e,s}^*$
in the viscoelastic solutions, resulting in the experimental phenomenon shown in figure 3(a).
The relative velocity and shear rate (
$\gamma ^*=\partial u^*/\partial (2\alpha /D)$
) along the diagonal of the channel are plotted in figure 4(a,b). The presence of the particle perturbs the normal velocity distribution of the fluids, causing an imbalance in the shear rate on either side of the particle, which is well known to affect inertial lift. The case of
$u_s^*\gt 0$
causes the particle to have a velocity that leads the fluid, inducing an increase in the gap in shear rate between the two sides of the particle, which leads to an increase in
$F_i$
. On the other hand,
$u_s^*\lt 0$
causes the particle to lag behind the fluid, inducing the shear rate difference between the two sides of the particle to reverse, which leads
$F_i$
to change from positive to negative values. We investigated the deformation of the polymer (
$r^*=r/r_0$
) around the particle and the projection of dimensionless surface force density on the particle along the direction of
$\alpha$
(
$\boldsymbol{\tau }_{\alpha }^*=\boldsymbol{F}^*\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\hat \alpha }/(\pi \kappa ^2)$
, where
$\boldsymbol{n}$
and
$\hat \alpha$
are the unit normal vectors along the spherical normal and
$\alpha$
direction), as shown in figure 4(c).The upper side of the particle corresponds to the wall-facing region, and the lower side corresponds to the channel centre. Upstream of the particle, the polymer near the wall is compressed (
$r^*\lt 0$
), while the polymer near the channel centre is stretched (
$r^*\gt 0$
); downstream, the deformation is reversed. Compressed polymer produces a pushing force toward the COM (
$\boldsymbol{\tau }_{\alpha }^*\lt 0$
), whereas a stretched polymer produces a pulling force away from it (
$\boldsymbol{\tau }_{\alpha }^*\gt 0$
).
The dimensionless velocity (a) and shear rate (b) along the diagonal of the channel. (c) The projection of dimensionless surface force density on the particle along the direction of
$\alpha$
(
$\boldsymbol{\tau }_{\alpha }^*=\boldsymbol{F}^*\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\hat \alpha }/(\pi \kappa ^2)$
) and polymer deformation (
$r^*=r/r_0$
) under (i)
$u_s^*=u_s/u_m=0$
, (ii)
$u_s^*\gt 0$
and (iii)
$u_s^*\lt 0$
. The upper side of the particle is near the wall, and the lower side is near the centre of the channel. In (c),
$r^*\lt 0$
indicates that the polymer is compressed and it causes an elastic force on the sphere in the direction of the particle COM (
$\boldsymbol{\tau }_{\alpha }^*\lt 0$
), and the blue part
$r^*\gt 0$
indicates that the polymer is stretched and exerts a pulling force on the particles (
$\boldsymbol{\tau }_{\alpha }^*\gt 0$
).
When
$u_s^*\gt 0$
, compression of the polymer near the wall becomes more pronounced, and the affected region expands, exerting a stronger downward thrust on the particle. Meanwhile, compression near the centre-facing side is significantly reduced, which diminishes the upward thrust there. Although there are stretched regions on both sides, they have a relatively minor effect compared with the larger compressed area. The combined outcome is an additional elastic lift force that drives the particle toward the channel centre. By contrast, when
$u_s^*\lt 0$
, the deformation pattern is reversed and the resulting elastic lift pushes the particle toward the wall. Although the reported polymer deformation is only approximately
$0.001$
, which might appear small, the slip velocity generates a significant velocity gradient within an interfacial layer surrounding the particle. This gradient induces highly inhomogeneous polymer conformations and stress distributions across the surface. Specifically, the polymer near the wall and the centre of the microchannel undergoes different levels of stretching and orientation. Consequently, surface traction varies significantly across the particle’s surface, resulting in notable changes in force integration at all points.
To study the dependence of
$F_{e,s}$
on slip velocity
$u_s$
, figure 5 compares
$F_{e,s}^*$
under different
$u_s^*$
with
$\kappa =2a/W=0.1$
,
${Re}=2.09$
and
${Wi}=1.64$
. Consistent with the aforementioned results, the direction of the slip-induced elastic lift is determined by the direction of
$u_s^*$
, while it grows with
$u_s^*$
. Taking
$u_s^*\gt 0$
as an example,
$F_{e,s}^*$
increases monotonically, with a fitted exponent of
${\sim} 0.98$
, demonstrating a linear relationship with
$u_s^*$
, as illustrated in figure 5(b, c). Figure 6 demonstrates that
$F_{e,s}^*$
also increases monotonically with
$\kappa$
at fixed
$|u_s^*|=0.04273$
,
${Re}=2.09$
and
${Wi}=1.64$
. The fitted exponent of
${\sim} 0.91$
align closely with theoretical predictions, where we expect
$F_{e,s}\sim u_s$
and
$F_{e,s}\sim \kappa$
.
Dimensionless slip-induced elastic lift (
$F_{e,s}^*=F_{e,s}/(W\eta _su_m)$
) under different slip velocities (a, b) and its exponential dependence on slip velocity (c) for parameters
$\kappa =2a/W=0.1$
,
${Re}=2.09$
and
${Wi}=1.64$
.
Dimensionless slip-induced elastic lift (
$F_{e,s}^*=F_{e,s}/(W\eta _su_m)$
) under different block ratios (a, b) and its exponential dependence on block ratio (c) for parameters
$|u_s^*|=|u_s/u_m|=0.04273$
,
${Re}=2.09$
, and
${Wi}=1.64$
.
To probe the mesoscopic effect of polymer size, we varied chain length
$M$
with
$\kappa =0.1$
,
$|u_s^*|=0.04273$
and
${Re}=2.09$
at the same concentrations. Figure 7(a) shows that
$|F_{e,s}^*|$
increases monotonically with
$M$
on a power-law basis with a fitted exponent of
${\sim} 1.66$
. Taking the lift at
$M = 9$
as
$1$
, figure 7(b) normalises
$|F_{e,s}^*|$
by its value at
$M = 9$
, enabling direct comparison with the theoretical model ((2.10)). While the theory slightly underestimates the growth of
$|F_{e,s}^*|$
with
$M$
, the simulation results follow the fitted scaling closely. The discrepancy likely arises from the ideal-chain assumption in the theory, which neglects the excluded-volume effects that are captured in the CGMD simulations.
(a) Dimensionless slip-induced elastic lift (
$F_{e,s}^*=F_{e,s}/(W\eta _su_m)$
) under different polymer lengths for parameters
$\kappa =2a/W=0.1$
,
$|u_s^*|=|u_s/u_m|=0.04273$
and
${Re}=2.09$
. (b) Computational vs. theoretical results by (2.10) of dimensionless slip-induced elastic lift under different polymer lengths.
4. Discussion and conclusions
The present study elucidates the intricate interplay between slip-induced migration mechanisms and viscoelastic flow characteristics in microchannels. Our computational framework reveals that the slip velocity fundamentally alters particle migration patterns through competing inertial and elastic lift components. In Newtonian fluids, slip-induced inertial lift dominates, driving particles with forward slip toward the channel walls. Conversely, the slip-induced elastic lift becomes predominant due to polymer stress asymmetry in viscoelastic fluids, reversing the migration direction. This finding aligns seamlessly with the experimental observations of Li & Xuan (Reference Li and Xuan2018, Reference Li and Xuan2023) and our own. This reversal mechanism stems from differential polymer deformation patterns around the particles, where forward slip enhances compression-induced stresses near channel walls while reverse slip amplifies tensile stresses toward the centre.
The observed power-law dependence of elastic lift on polymer chain length suggests non-ideal-chain behaviour beyond Zimm model predictions. While our theoretical framework qualitatively captures the trend, the quantitative underestimation likely arises from neglected volume exclusion effects in the ideal-chain assumption. This discrepancy highlights the need for improved mesoscopic models incorporating excluded-volume interactions, particularly for longer polymer chains where crowding effects become significant. In addition, the Zimm model is only applicable to dilute solutions. Entanglement will occur between the polymer chains for concentrated solutions, and theories such as the tube model need to be further considered. In addition, changes in viscoelasticity can alter the velocity distribution of the flow field and may also cause deviations from the theoretical model.
Notably, the linear dependence of slip-induced elastic lift on both the slip velocity and the blocking ratio supports the theoretical scaling relationships (Einarsson & Mehlig Reference Einarsson and Mehlig2017; Khair & Kabarowski Reference Khair and Kabarowski2020). Although
$Re= 2.09$
exceeds the strict low-Reynolds-number limit, our results confirm that the slip-induced elastic lift preserves the predicted direction and scaling. These scaling laws persist across different flow regimes, suggesting broad applicability for predicting particle manipulation in microfluidic devices. The identified migration reversal mechanism provides new opportunities for viscoelasticity-enhanced particle sorting, where subtle adjustments in polymer concentration or chain length could enable precise spatial control of particles with different properties. Such a slip velocity can be induced by a variety of external force fields, including electric fields, magnetic fields, etc.
Moreover, the fundamental framework established in this work can be extended to study more complex polymer solutions, such as polyelectrolytes. Unlike commonly used electro-neutral PEO solutions, electrophoretic particles in polyelectrolyte solutions exhibit contrasting migration patterns, a phenomenon attributed to the polyelectrolyte effect (Serhatlioglu et al. Reference Serhatlioglu, Isiksacan, Özkan, Tuncel and Elbuken2020). This extension would further enrich our understanding of particle migration in diverse viscoelastic media.
Funding
This work was supported by the National Key Research and Development Program of China (Grant No. 2022YFA1203200) and the National Natural Science Foundation of China (Grant No. 12272345, 12472322).
Declaration of interests
The authors report no conflict of interest.
Author contributions
S. M. and X. Q. have contributed equally to the work.
Data availability statement
The data that support the findings of this study are available under reasonable request.
Open access
