1. Introduction
Microparticles suspended in fluids are ubiquitous in natural and engineered systems, underpinning processes such as blood flow and cell separation (Di Carlo Reference Di Carlo2009), sediment transport in geophysical flows (Meiburg & Kneller Reference Meiburg and Kneller2010) and particle-laden flows in engineered systems (Brandt & Coletti Reference Brandt and Coletti2022). The efficacy of many such processes depends critically on the transport and spatial distribution of suspended particles, governed by both their streamwise kinematics and their cross-stream migration. Among the various flow configurations used to investigate particle motion, plane Couette flow provides a simple yet representative model system. It captures the essential shear dynamics relevant to rheological devices, such as lubricated narrow gaps, and shear-based microfluidic systems, while retaining a sufficiently simple geometry to allow analytical treatment. We therefore adopt plane Couette flow as a canonical model configuration. In homogeneous shear flows, most studies of particle dynamics have focused on the weak-to-moderate inertial regime, where cross-stream migration arises from finite fluid inertia. A seminal theoretical result is that of Saffman (Reference Saffman1965), who obtained the lift on a slipping sphere in unbounded shear flow. For plane Couette flow confined between parallel walls, Ho & Leal (Reference Ho and Leal1974) showed within a point-particle framework that the centreline is the only stable equilibrium in the weak-inertia limit, a result confirmed by the more accurate calculations of Vasseur & Cox (Reference Vasseur and Cox1976). Subsequent numerical studies corroborated and extended these predictions. Direct simulations by Feng et al. (Reference Feng, Hu and Joseph1994) confirmed centreline focusing for a cylinder at moderate particle inertia, while lattice Boltzmann simulations of cylinders (Fox et al. Reference Fox, Schneider and Khair2020) and spheres (Fox et al. Reference Fox, Schneider and Khair2021) showed that increasing fluid inertia destabilises the centreline, producing two symmetric off-centre stable equilibria through a supercritical pitchfork bifurcation. Using a point-particle formulation, Anand & Subramanian (Reference Anand and Subramanian2023) showed that the onset of this instability is controlled by the channel rather than the particle Reynolds number, with critical values
${\approx} 148$
for a sphere and
${\approx} 110$
for a cylinder.
In realistic environments, however, the suspending medium is rarely homogeneous. Spatial variations of temperature or solute concentration commonly produce viscosity gradients, for example in the microvasculature and glycocalyx (Pries & Secomb Reference Pries and Secomb2005), mucus layers (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009) and biofilms (Wilking et al. Reference Wilking, Angelini, Seminara, Brenner and Weitz2011), as well as in co-flow systems (Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021) and temperature-controlled microdevices (Mao et al. Reference Mao, Yang and Cremer2002). Such gradients break the symmetries of constant-viscosity flows and introduce hydrodynamic couplings absent in homogeneous fluids (Sharma et al. Reference Sharma, Thakur, Jose and Govindarajan2026): forces couple to particle rotation and torques to translation, yielding novel dynamics for spheres (Datt & Elfring Reference Datt and Elfring2019), spheroids (Anand & Narsimhan Reference Anand and Narsimhan2024) and slender bodies (Kamal & Lauga Reference Kamal and Lauga2023). Recently, Sharma et al. (Reference Sharma, Bosler, Govindarajan and Koch2025) obtained the velocity slip of freely suspended spheres and spheroids in an unbounded shear flow with linearly stratified viscosity. Wall-driven flows such as plane Couette flow are fundamentally different: because the shear stress is uniform across the channel, a linearly stratified viscosity necessarily produces a curved undisturbed velocity profile, whose consequences for rigid-particle motion remain unexplored. Moreover, although viscosity gradients induce cross-stream migration of deformable capsules in Stokesian Couette flow (Laumann & Zimmermann Reference Laumann and Zimmermann2019), the behaviour of rigid particles, which lack deformation-induced symmetry breaking, is not yet understood. The present work addresses this gap by analysing a rigid sphere in viscosity-stratified plane Couette flow, both in the Stokes limit and at non-zero Reynolds numbers. Hereafter, dimensional (dimensionless) variables are denoted with (without) a tilde.
2. Problem formulation
A rigid, neutrally buoyant spherical particle of radius
$\tilde {R}$
is suspended at position
$\boldsymbol{\tilde {X}} = (\tilde {X},\tilde {Y},\tilde {Z})$
in an incompressible Newtonian fluid confined between two infinite parallel plates separated by a distance
$\tilde {H}$
(figure 1
a). The confinement ratio is
$\kappa = \tilde {R}/\tilde {H}$
. The upper plate translates with velocity
$\tilde {U}_w\boldsymbol{e}_x$
and the lower plate is stationary; no-slip boundary conditions hold at both walls, and the streamwise (
$\tilde {x}$
) and spanwise (
$\tilde {z}$
) directions are periodic. Driven by the imposed flow, the particle translates with velocity
$\boldsymbol{\tilde {U}} = (\tilde {U}_x, \tilde {U}_y, \tilde {U}_z)$
, and rotates with angular velocity
$\boldsymbol{\tilde {\varOmega }} = (\tilde {\varOmega }_x, \tilde {\varOmega }_y, \tilde {\varOmega }_z)$
. The fluid satisfies the variable-viscosity Navier–Stokes equations
\begin{equation} \boldsymbol{\tilde{\nabla }}\boldsymbol{\cdot }{\boldsymbol{\tilde {u}}}=0,\quad \tilde \rho \big (\partial _{\tilde t}\,\boldsymbol{\tilde {u}} +\boldsymbol{\tilde {u}} \boldsymbol{\cdot }\boldsymbol{\tilde{\nabla }} \boldsymbol{\tilde {u}}\big) =-\boldsymbol{\tilde {\nabla }}\tilde p +\boldsymbol{\tilde{\nabla }}\boldsymbol{\cdot } \big [\tilde \eta ( \boldsymbol{\tilde {x}}) \big(\boldsymbol{\tilde{\nabla }}\boldsymbol{\tilde {u}} +\boldsymbol{\tilde {\nabla }}\boldsymbol{\tilde {u}}^{\mathrm{T}}\big)\big ], \end{equation}
where
$\tilde {\rho }$
is the fluid density,
$\tilde {\eta } (\boldsymbol{\tilde {x}} )$
is the spatially varying fluid viscosity and
$(\tilde {p},\boldsymbol{\tilde {u}})$
denote the pressure and velocity fields, respectively. On the particle surface
$\tilde {\mathcal{S}}_p$
, the no-slip condition requires
$\boldsymbol{\tilde {u}}=\boldsymbol{\tilde {U}}+\boldsymbol{\tilde {\varOmega }}\times \boldsymbol{(\tilde {x} - \tilde {X})}$
. We solve for the particle kinematics,
$\boldsymbol{\tilde {U}}$
and
$\boldsymbol{\tilde {\varOmega }}$
, both in the Stokes limit and at non-zero Reynolds numbers.
(a) Schematic of a spherical particle in plane Couette flow with a prescribed viscosity field
$\tilde {\eta }(\tilde {\boldsymbol{x}})$
, represented by the background colour gradient. (b) Multi-direct-forcing immersed-boundary lattice Boltzmann methodology, showing the computational domain with Eulerian fluid nodes, no-slip wall nodes and Lagrangian particle markers.
Figure 1. Long description
Panel A: A schematic of a spherical particle in plane Couette flow. The particle is shown within a fluid represented by a background color gradient indicating the viscosity field. The flow is confined between two no-slip walls, one moving and one fixed, with arrows indicating the direction of fluid and particle movement. Panel B: A diagram of the multi-direct-forcing immersed-boundary lattice Boltzmann methodology. The computational domain includes Eulerian fluid nodes, no-slip wall nodes, and Lagrangian particle markers. The periodic boundaries and the distribution of nodes and markers are shown, with arrows indicating the direction of fluid flow.
The fluid viscosity is prescribed to vary linearly across the channel as
$\tilde {\eta }(\boldsymbol{\tilde {x}})= \tilde {\eta }_0(1+\epsilon \,\tilde {y}/\tilde {H})$
, where
$\epsilon$
is the dimensionless viscosity gradient and
$\tilde \eta _0\equiv \tilde \eta (0)$
the reference viscosity. Such profiles arise frequently in microfluidic gradient generators (Jeon et al. Reference Jeon, Dertinger, Chiu, Choi, Stroock and Whitesides2000), thermally stratified flows (Mao et al. Reference Mao, Yang and Cremer2002) and the microvasculature (Fedosov et al. Reference Fedosov, Caswell, Popel and Karniadakis2010). The flow is characterised by the particle Reynolds number
$\textit{Re}_p := \tilde {\rho }\tilde {G}\tilde {R}^2/\tilde \eta _0$
and the channel Reynolds number
$\textit{Re}_c := \textit{Re}_p/\kappa ^{2}$
, where
$\tilde {G}=\tilde {U}_w/\tilde {H}$
is the nominal shear rate. Lengths, velocities and times are scaled by
$\tilde H$
,
$\tilde U_w$
and
$\tilde {H}/\tilde {U}_w$
. We consider
$\epsilon \in [-0.5,1]$
. A representative value
$\epsilon =0.5$
gives a wall-to-wall viscosity ratio
$\tilde \eta (\tilde H)/\tilde \eta (0)=1.5$
, i.e. a dimensional gradient of
$5\ \mathrm{Pa\,s\,m^{-1}}$
in a
$100$
μm channel filled with an aqueous solution (
$\tilde \eta _0=1$
mPa s). Such modest stratifications are experimentally achievable through, for example, a
${\approx} 15\,\%$
transverse difference in glycerol mass fraction (Cheng Reference Cheng2008), a
${\approx} 20$
K temperature difference in water (Kestin et al. Reference Kestin, Sokolov and Wakeham1978) or co-flowing fluids in Y-shaped channels (Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021). We take the viscosity profile to be frozen during focusing. For a solute-driven gradient (e.g. glycerol in water, diffusivity
$\tilde D\approx 10^{-9}\ \mathrm{m^2\,s^{-1}}$
), the transverse diffusion time
$\tilde t_D\sim \tilde H^2/\tilde D\approx 10$
s exceeds the focusing time over
$\textit{Re}_c\in [10,250]$
(
${\approx} 0.8$
s to
$12$
ms) by at least an order of magnitude, justifying a stationary viscosity field. Thermal diffusivity, by contrast, is about two orders of magnitude larger, so a thermal gradient relaxes on time scales comparable to or shorter than focusing and must be sustained by the wall temperatures. We also assume that the particle does not disturb the viscosity field; the effect of relaxing this assumption is discussed in subsequent sections.
3. Solution methodology
3.1. Reciprocal theorem in the Stokes limit
In the creeping-flow limit, the system reduces to the variable-viscosity Stokes equations. To obtain the particle velocities, we apply the Lorentz reciprocal theorem (Masoud & Stone Reference Masoud and Stone2019) using an auxiliary problem of a sphere translating (
$\boldsymbol{U}_0$
) and rotating (
$\boldsymbol{\varOmega }_0$
) in an unbounded fluid of constant viscosity
$\tilde {\eta }_0$
, yielding
\begin{equation} \boldsymbol{F}_0\boldsymbol{\cdot }\boldsymbol{U} +\boldsymbol{T}_0\boldsymbol{\cdot }\boldsymbol{\varOmega } - \int _{\mathcal{S}_{p}}(\boldsymbol{\sigma }_0\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{u}_\infty \,\mathrm{d} \mathcal{S} = 2\int _{\mathcal{V}}\left (\frac {\tilde {\eta }}{\tilde {\eta }_0} - 1 \right )\,\boldsymbol{E}' : \boldsymbol{E}_0\,\mathrm{d} \mathcal{V}. \end{equation}
This exact identity holds for a force- and torque-free particle; the subscript ‘
$0$
’ denotes the auxiliary problem, with stress
$\boldsymbol \sigma _0$
, force
$\boldsymbol F_0$
and torque
$\boldsymbol T_0$
,
$\boldsymbol n$
is the outward normal,
$\boldsymbol u_\infty$
is the undisturbed Couette profile in the stratified fluid,
$\mathcal S_p$
is the particle surface and
$\mathcal V$
is the fluid volume. In (3.1),
$\boldsymbol{ E}'$
is the rate of strain of the disturbance field
$\boldsymbol{ u}-\boldsymbol{ u}_\infty$
, which for a force- and torque-free sphere decays as a stresslet and so renders the volume integral convergent, whereas the total rate of strain would diverge. We treat the domain as unbounded (
$ \kappa \ll 1$
), so that the walls enter only through the undisturbed flow. Because the volume integrand of (3.1), weighted by the volume element, decays only as
$r^{-2}$
(
$(\tilde {\eta }/\tilde {\eta }_0-1)\sim r$
,
$\boldsymbol{E}'\sim r^{-3}$
,
$\boldsymbol{E}_0\sim r^{-2}$
and
$\mathrm{d}\mathcal{V}\sim r^2\,\mathrm{d}r$
), distances
$\sim H$
contribute an
$\mathcal{O}(\kappa )$
fraction. The analysis below is therefore the leading term of an expansion in
$\kappa$
, with
$\mathcal{O}(\kappa )$
wall-reflection corrections. We further assume weak gradients (
$|\epsilon |\ll 1$
), for which we expand the dynamic variables as:
$\{\boldsymbol{{U}}, \boldsymbol{{\varOmega }}, \boldsymbol{{u}}_\infty , \boldsymbol{{E}}^\prime \}=\{\boldsymbol{{U}}, \boldsymbol{{\varOmega }}, \boldsymbol{{u}}_\infty , \boldsymbol{{E}}^\prime \}^{(0)}+\epsilon \{\boldsymbol{{U}}, \boldsymbol{{\varOmega }}, \boldsymbol{{u}}_\infty , \boldsymbol{{E}}^\prime \}^{(1)}+\mathcal{O}(\epsilon ^{2})$
. Equation (3.1) then splits into a leading-order balance that recovers the classical Faxén relations and an
$\mathcal{O}(\epsilon )$
balance driven entirely by the leading-order disturbance, each solvable in closed form.
3.2. Immersed-boundary lattice Boltzmann method at non-zero Reynolds number
To relax the assumptions of weak confinement, weak gradients and vanishing inertia, we develop an in-house three-dimensional immersed-boundary lattice Boltzmann method (IB-LBM) solver, following Krüger et al. (Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017) for the lattice Boltzmann algorithm, Suzuki & Inamuro (Reference Suzuki and Inamuro2011) for the multi-direct forcing and the treatment of the fluid enclosed by the particle and Inamuro et al. (Reference Inamuro, Yoshino and Suzuki2021) for the rigid-body coupling. In lattice units (
$\Delta \hat x=\Delta \hat t=1$
, hats denoting lattice quantities), the distribution functions
$\hat f_i$
associated with the D3Q19 velocity set
$\hat {\boldsymbol c}_i$
obey the single-relaxation-time (Bhatnagar–Gross–Krook) equation:
\begin{equation} \hat f_i(\hat {\boldsymbol x}+\hat {\boldsymbol c}_i,\hat t+1)-\hat f_i(\hat {\boldsymbol x},\hat t)=-\frac {\hat f_i-\hat f_i^{\textit{eq}}}{\hat \tau (\hat y)}+\hat S_i,\qquad \hat \tau (\hat y)=\frac 12+\frac {\hat \nu _0}{\hat c_s^{2}}\left (1+\epsilon \frac {\hat y}{\hat H}\right )\!, \end{equation}
where
$\hat f_i^{\textit{eq}}=w_i\hat \rho \,[1+\hat {\boldsymbol c}_i\boldsymbol{\cdot }\hat {\boldsymbol u}/\hat c_s^2+(\hat {\boldsymbol c}_i\boldsymbol{\cdot }\hat {\boldsymbol u})^2/(2\hat c_s^4)-\hat {\boldsymbol u}\boldsymbol{\cdot }\hat {\boldsymbol u}/(2\hat c_s^2)]$
is the second-order equilibrium distribution with lattice weights
$w_i$
,
$\hat S_i$
is the forcing term of Guo et al. (Reference Guo, Zheng and Shi2002) carrying the immersed-boundary force density
$\hat {\boldsymbol F}$
,
$\hat c_s=1/\sqrt 3$
and
$\hat \nu _0=\hat G\hat R^{2}/\textit{Re}_p$
; the density and velocity follow as
$\hat \rho =\sum _i\hat f_i$
and
$\hat \rho \hat {\boldsymbol u}=\sum _i\hat {\boldsymbol c}_i\hat f_i+\hat {\boldsymbol F}/2$
. A Chapman–Enskog expansion shows that (3.2) recovers (2.1) with
$\hat \eta =\hat \rho \hat c_s^{2}(\hat \tau -1/2)$
, i.e. exactly the prescribed linear viscosity, up to
$\mathcal{O}({Ma}^2)$
errors, with
${Ma}=\hat U_w/\hat c_s$
(Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017; Laumann & Zimmermann Reference Laumann and Zimmermann2019); because
$\hat \tau$
enters the collision term algebraically, its spatial variation introduces no additional terms. The streamwise and spanwise directions are periodic, and no slip is enforced at the walls by halfway bounce-back,
$\hat f_{\bar i}=\hat f_i^{\ast }-2w_i\hat \rho \,\hat {\boldsymbol c}_i\boldsymbol{\cdot }\hat {\boldsymbol u}_w/\hat c_s^{2}$
, where
$\hat f_i^{\ast }$
is the post-collision distribution,
$\bar i$
the direction opposite to
$i$
and the last term, the momentum correction of Ladd (Reference Ladd1994), imparts the wall velocity
$\hat {\boldsymbol u}_w$
(see figure 1
b). The sphere is represented by Lagrangian markers
$\boldsymbol X_l$
on a surface Fibonacci lattice with spacing
${\approx} \Delta \hat x$
, supplemented by interior markers that enforce rigid motion of the enclosed fluid and thereby remove its spurious inertia (Suzuki & Inamuro Reference Suzuki and Inamuro2011). At each time step, the fluid velocity
$\boldsymbol U_l$
is interpolated to every marker with Peskin’s four-point regularised delta function, the marker force is set to
$\hat {\boldsymbol g}_l=(\boldsymbol U_l^{d}-\boldsymbol U_l)/\Delta \hat t$
, with
$\boldsymbol U_l^{\mathrm d}=\boldsymbol U+\boldsymbol \varOmega \times (\boldsymbol X_l-\boldsymbol X)$
the rigid-body velocity, and spread back to the lattice with the same kernel; repeating this multi-direct-forcing correction
$15$
times per step enforces the no-slip condition at every marker. The hydrodynamic force and torque, obtained from the marker forces together with the rate of change of momentum of the enclosed fluid, advance the neutrally buoyant particle through Newton’s equations of rigid-body motion. The range
$\textit{Re}_c\in [10,250]$
is spanned by adjusting
$\hat U_w$
and
$\hat \tau$
together, keeping
${Ma}\lesssim 0.1$
and
$\hat \tau _{min }\geqslant 0.53$
, with
$150\times 75\times 150$
nodes (
$\hat H=75$
,
$\hat R=15$
) for
$\textit{Re}_c\leqslant 100$
and
$300\times 150\times 300$
nodes (
$\hat H=150$
,
$\hat R=30$
) for
$\textit{Re}_c\gt 100$
at
$\kappa =0.2$
. Stable equilibria are obtained from free-particle trajectories. The lift profiles (
$F_y$
) are computed by holding the particle at fixed
$Y$
while leaving it free to translate streamwise and rotate, and unstable equilibria are identified as zero crossings of
$F_y$
with
$\partial _Y F_y\gt 0$
. Without the particle, the solver reproduces the Couette profile with the shear stress uniform to within
$1\,\%$
for
$\epsilon \in [-0.5,1]$
, a stringent test of the variable relaxation time. With the particle, it reproduces the homogeneous (
$\epsilon =0$
) trajectories of Fox et al. (Reference Fox, Schneider and Khair2021) over
$\textit{Re}_p\in [1,30]$
(figure 2
a), and the trajectories are insensitive to resolution, domain size and Mach number (figure 2
b).
(a) Validation of inertial-lift trajectories for a freely suspended, neutrally buoyant sphere in constant-viscosity plane Couette flow at
$\kappa =0.2$
. Symbols: Fox et al. (Reference Fox, Schneider and Khair2021); lines: present IB-LBM. (b) Convergence of the transverse trajectory under grid refinement (GIT), periodic-length (LIT) and shear-rate (SIT) variation at
$\textit{Re}_c=250$
,
$\epsilon =1$
and
$\kappa =0.2$
.
4. Results and discussion
4.1. Stokes regime hydrodynamics
Solving (3.1) as outlined above gives
\begin{align} \boldsymbol{U} &= \Biggl [ \,\underbrace {\left (\boldsymbol{u}_\infty ^{(0)} + \frac {R^2}{6} {\nabla} ^2 \boldsymbol{u}_\infty ^{(0)} \right )}_{\text{Leading order}} +\, \epsilon \, \underbrace {\left (\boldsymbol{u}_\infty ^{(1)} - \frac {R^2}{3} {\nabla} ^2 \boldsymbol{u}_\infty ^{(1)} \right )}_{\mathcal{O}( \epsilon )}\, \Biggr ]_{\!\boldsymbol{X}} + \mathcal{O}\big ( \epsilon ^2\big ), \\[-12pt]\nonumber \end{align}
\begin{align} \boldsymbol{\varOmega } &= \Biggl [ \, \underbrace {\frac {1}{2} \big ( \boldsymbol{\nabla } \times \boldsymbol{u}_\infty ^{(0)} \big )}_{\text{Leading order}} +\,\epsilon \, \underbrace {\frac {1}{2} \big ( \boldsymbol{\nabla } \times \boldsymbol{u}_\infty ^{(1)} \big )}_{\mathcal{O}( \epsilon )} \Biggr ]_{\!\boldsymbol{X}} + \mathcal{O}\big ( \epsilon ^2\big ), \end{align}
where
$\boldsymbol{u}_\infty = u_\infty \boldsymbol{e}_x = [\ln (1+\epsilon y)/\ln (1+\epsilon )]\boldsymbol{e}_x \approx [y- ({1}/{2})\epsilon (y^2-y)]\boldsymbol{e}_x = \boldsymbol{u}^{(0)}_\infty + \epsilon \boldsymbol{u}^{(1)}_\infty$
is the Couette velocity profile modified by the viscosity gradient. The exact profile is logarithmic; consistent with
$|\epsilon |\ll 1$
, we retain it to
$\mathcal{O}(\epsilon )$
. Since no pressure gradient is imposed, the shear stress
$\eta \,\mathrm{d} u_\infty /\mathrm{d} y$
is uniform across the channel. Both signs of
$\epsilon$
are considered, where
$\epsilon \gt 0$
(
$\epsilon \lt 0$
) corresponds to the viscosity increasing towards the moving (stationary) wall. The leading-order terms in (4.1) recover the classical Faxén relations, whereas the
$\mathcal{O}(\epsilon )$
corrections arise solely from the viscosity gradient. We note that the compact notation of (4.1) is specific to this configuration and is not a generalised Faxén law. To elucidate the physical origins of the particle motion, we separate (4.1a
) into two parts, corresponding to the surface and volume integrals of the reciprocal theorem (3.1), respectively, scaled by the auxiliary Stokes drag:
\begin{equation} U_x - u_{\infty }(\boldsymbol{X}) = \frac {R^{2}}{3}\,\epsilon = \;\; \underbrace {-\,\frac {R^{2}}{6}\,\epsilon }_{\text{Faxen kinematic lag}} + \underbrace {\frac {R^{2}}{2}\,\epsilon }_{\text{Gradient}{-}\text{disturbance coupling}} \!. \end{equation}
Equation (4.2) demonstrates that the net velocity slip (the difference between the particle velocity and the undisturbed fluid velocity at its centre) results from two competing mechanisms, reflecting the gradient-induced modification of the undisturbed flow and of the disturbance flow, respectively. (i) Faxén lag: Uniform shear stress across the channel dictates that the local shear rate scales as
$1/\eta$
. This bends the undisturbed velocity profile, concentrating shear in the less-viscous region and imparting a constant curvature (
${\nabla} ^2 u_{\infty } = -\epsilon$
). By Faxén’s first law, the spherical particle translates with the surface average of the ambient velocity, which is smaller (larger) than its centre value for
$\epsilon \gt 0$
(
$\epsilon \lt 0$
). This gives a kinematic slip of
$-\epsilon R^{2}/6$
whose sign is opposite to that of the viscosity gradient. (ii) Gradient–disturbance coupling: Competing against this kinematic lag is the coupling between the viscosity gradient and the particle’s own disturbance flow. In a constant-viscosity fluid, the disturbance stresses generated by the particle in shear flow exactly cancel. However, when viscosity increases towards the moving wall, the stress exerted by the overtaking upper fluid layer, which pushes the particle downstream, exceeds the drag from the lower layer opposing its motion. This asymmetry generates a net streamwise driving force, yielding a co-directional slip (
$+\epsilon R^2/2$
) that takes the sign of the viscosity gradient. This coupling term coincides with the unbounded linear-flow result of Sharma et al. (Reference Sharma, Bosler, Govindarajan and Koch2025), where the undisturbed flow has no curvature and hence no Faxén lag. In wall-driven stratified flows, however, the curvature is unavoidable, and because the two mechanisms oppose each other, neglecting the lag would overestimate the slip by
$50\,\%$
.
The velocity slip is independent of particle position, grows quadratically with particle size (vanishing for a point particle) and takes the sign of
$\epsilon$
: the particle leads the local flow when the viscosity increases towards the moving wall and lags when it decreases. It is entirely viscosity-gradient-induced, since a neutrally buoyant sphere in constant-viscosity Couette flow translates with the local fluid. Similarly, (4.1b
) gives
$\varOmega _z=- (1/2)[\partial _y u_\infty ]_{\boldsymbol X} =- 1/2+ (1/4)\epsilon (2Y-1)$
: the particle rotates at half the local vorticity of the gradient-modified flow, without any gradient-induced relative rotation. Figure 3 compares these predictions with IB-LBM results for both signs of the gradient. At
$\kappa =0.1$
, the variation of
$U_x$
in figure 3(a) is dominated by the change of the undisturbed velocity at the particle centre,
$\epsilon u_\infty ^{(1)}(Y)$
, while the slip in (4.2), arising from both the undisturbed-flow curvature and the disturbance flow, contributes only
$\mathcal{O}(\epsilon \kappa ^2)$
. The small deviation at the largest gradient (
${\approx} 1.2\times 10^{-3}$
at
$\epsilon =0.1$
) is consistent in magnitude with the
$\mathcal{O}(\epsilon ^2)$
truncation error of the undisturbed flow (
$|u_\infty - u_\infty ^{(0)}-\epsilon u_\infty ^{(1)}|\approx 0.7\times 10^{-3}$
at
$|\epsilon |=0.1$
,
$Y=0.55$
). In figure 3(b), the IB-LBM results reproduce the predicted slope
$\mathrm{d}\varOmega _z/\mathrm{d}\epsilon =(2Y-1)/4$
but are offset by an
$\epsilon$
-independent
${\approx} 0.7\,\%$
. This offset is a confinement effect absent from the unbounded theory: the leading far-field wall correction reduces the rotation of a torque-free sphere in shear by
${\approx} (5/16)(R/h)^3$
per wall at distance
$h$
(Goldman et al. Reference Goldman, Cox and Brenner1967), i.e. by
$0.53\,\%$
for walls at
$Y=0.55$
and
$1-Y$
, which accounts for most of the offset. Since this correction is
$\mathcal{O}(\kappa ^3)$
, with gradient-induced modifications only at
$\mathcal{O}(\epsilon \kappa ^3)$
, it does not affect the
$\mathcal{O}(\epsilon )$
predictions. Moreover, (4.1a
) yields no cross-stream velocity, and this absence of migration,
$U_y\equiv 0$
, is exact rather than a consequence of the
$\mathcal{O}(\epsilon )$
approximation. By linearity of the Stokes equations, reversing the wall speed reverses the transverse velocity,
$U_y(-U_w)=-U_y(U_w)$
. Yet the reflection
$x\mapsto -x$
leaves the sphere, the walls and
$\eta (y)$
invariant while mapping the problem with wall speed
$U_w$
onto that with
$-U_w$
, without altering
$y$
; hence
$U_y(-U_w)=U_y(U_w)$
. Both conditions hold only if
$U_y\equiv 0$
, for any confinement ratio, gradient magnitude and transverse viscosity profile.
Comparison of the
$\mathcal{O}(\epsilon )$
analytical predictions (lines) with IB-LBM results (symbols) for (a)
$U_x$
and (b)
$\varOmega _z$
at
$\textit{Re}_p=0.01$
and
$\kappa =0.1$
(
$\textit{Re}_c=1$
), for a particle at
$Y=0.55$
. Inset of (a): schematic of the gradient-modified shear profile.
Three remarks follow. First, the fore–aft symmetry of the particle shape is essential: deformable particles, whose shapes do not reverse with the flow, break this symmetry and migrate across streamlines (Laumann & Zimmermann Reference Laumann and Zimmermann2019). Second, the translation–rotation couplings of Datt & Elfring (Reference Datt and Elfring2019), a torque
$\propto \boldsymbol U_{\textit{rel}}\times \boldsymbol {\nabla }\eta$
and a force
$\propto \boldsymbol \varOmega _{\textit{rel}}\times \boldsymbol {\nabla }\eta$
, where
$\boldsymbol U_{\textit{rel}}$
and
$\boldsymbol \varOmega _{\textit{rel}}$
are the velocity slip and relative rotation, enter only at
$\mathcal{O}(\epsilon ^2)$
, since the velocity slip is
$\mathcal{O}(\epsilon )$
and the relative rotation vanishes at this order. Third, we have assumed that the particle leaves the viscosity-controlling scalar (solute concentration or temperature) unperturbed, which holds exactly for matched solute diffusivities or thermal conductivities. Otherwise, when advective transport of scalar is negligible (i.e. at low particle Péclet number,
$Pe_p=\tilde G\tilde R^2/\tilde D \ll 1$
, with
$\tilde D$
the diffusivity of the scalar), the scalar, and hence the viscosity if it depends linearly on the scalar, is harmonic, so that
$\eta (\boldsymbol x)=1+\epsilon y+\eta '({\boldsymbol x})$
with
${\nabla} ^2\eta '=0$
. One can show that the particle generates a dipolar viscosity disturbance,
$\eta '=\epsilon \chi \,(y-Y)R^3/r^3$
, where
$\chi =1/2$
for a solute-impermeable particle and
$\chi =(\tilde k-\tilde k_p)/(2\tilde k+\tilde k_p)$
for fluid and particle thermal conductivities
$\tilde k$
and
$\tilde k_p$
(Batchelor Reference Batchelor1974). The additional stress
$2\eta '({\boldsymbol E}_\infty +{\boldsymbol E}')$
, with
${\boldsymbol E}_\infty$
the undisturbed rate of strain, adds a velocity slip
$\epsilon \chi R^2/4$
to (4.2), raising it to
$11\epsilon R^2/24$
for a solute-impermeable particle, whereas
$U_y$
and the relative rotation still vanish, by the reflection symmetry above and because no pseudovector is linear in both
$\boldsymbol {\nabla }\eta$
and
$\boldsymbol E_\infty$
.
4.2. Effect of non-zero Reynolds number
At non-zero Reynolds numbers, cross-stream migration at constant viscosity arises solely through wall-induced inertial lift, as the linear Couette profile fundamentally lacks shear-gradient lift (Ho & Leal Reference Ho and Leal1974; Vasseur & Cox Reference Vasseur and Cox1976). This wall-induced force focuses the particle at the centreline at low-to-moderate
$\textit{Re}_c$
, but reverses near the centreline at higher
$\textit{Re}_c$
, driving the particle towards off-centre equilibria (Fox et al. Reference Fox, Schneider and Khair2020, Reference Fox, Schneider and Khair2021; Anand & Subramanian Reference Anand and Subramanian2023). Figure 4 presents computed trajectories for
$\epsilon = 0,\ \pm 0.5$
and
$1$
at relatively low (
$\textit{Re}_c=10$
) and moderately high (
$\textit{Re}_c=250$
) inertia. For a constant-viscosity fluid (
$\epsilon =0$
), the system features a single stable centreline equilibrium at
$\textit{Re}_c=10$
(figure 4
a). By
$\textit{Re}_c=250$
, this centreline equilibrium has destabilised, yielding two symmetric off-centre equilibria (figure 4
b). Introducing a viscosity gradient systematically displaces these equilibria towards the less-viscous region. Crucially, at
$\textit{Re}_c=250$
, strong gradients (
$\epsilon =0.5, 1$
) force the particle to stabilise exclusively in the bottom half of the channel, irrespective of its initial position, whereas for
$\epsilon =-0.5$
, two distinct stable off-centre equilibria persist. We also observe that inertia modifies the particle’s streamwise translational velocity and rotation: at stable equilibria,
$U_x$
differs from the local undisturbed velocity only by a small slip, as in the Stokes limit, so its variation with
$\textit{Re}_c$
and
$\epsilon$
primarily follows that of
$Y^{\mathrm{eq}}$
, while
$|\varOmega _z|$
increasingly falls below half the local vorticity with increasing
$\textit{Re}_c$
(Mikulencak & Morris Reference Mikulencak and Morris2004), by a fraction that is nearly independent of
$\epsilon$
. For brevity, these results are not shown here. The primary effect of inertia, however, is to restructure the transverse equilibria.
Particle trajectories from different initial positions in a constant-viscosity fluid (
$\epsilon =0$
) and under viscosity gradients (
$\epsilon =\pm 0.5,\,1$
) for (a)
$\textit{Re}_c=10$
and (b)
$\textit{Re}_c=250$
, with
$\kappa =0.2$
. The channel centreline is at
$Y = 0.5$
.
Figure 4. Long description
Panel A: A line graph showing particle trajectories in a constant-viscosity fluid and under viscosity gradients. The horizontal axis is labeled t and ranges from 0 to 1000. The vertical axis is labeled Y and ranges from 0.2 to 0.8. The graph includes four lines representing different values of epsilon (0, 0.5, 1, and -0.5). The lines show the migration of particles due to viscosity gradients, with some particles moving towards the upper half and others towards the lower half of the channel. Panel B: A line graph showing particle trajectories in a constant-viscosity fluid and under viscosity gradients. The horizontal axis is labeled t and ranges from 0 to 500. The vertical axis is labeled Y and ranges from 0.2 to 0.8. The graph includes four lines representing different values of epsilon (0, 0.5, 1, and -0.5). The lines show the migration of particles due to viscosity gradients, with some particles moving towards the upper half and others towards the lower half of the channel.
To characterise the dynamics fully, figure 5(a) shows the equilibrium positions
$Y^{\mathrm{eq}}$
against
$\textit{Re}_c$
. At constant viscosity, the stable centreline equilibrium gives way to a pair of symmetric off-centre equilibria at
$\textit{Re}_c^{\mathrm{crit}}\approx 151$
for
$\kappa =0.2$
, close to the value
${\approx} 148$
obtained by Anand & Subramanian (Reference Anand and Subramanian2023) with a point-particle formulation valid for
$\textit{Re}_p\ll 1$
and
$\kappa \ll 1$
. A viscosity gradient shifts the equilibria monotonically towards the less-viscous region, consistent with figure 4. For a mild gradient,
$\epsilon =0.25$
(down-triangles), an additional pair of stable and unstable equilibria emerges at higher inertia, whereas for
$\epsilon =0.5$
and
$1$
(squares and up-triangles) a single equilibrium persists over the entire range
$\textit{Re}_c\in [10,250]$
. The gradient thus acts as a symmetry-breaking tilt that turns the supercritical pitchfork into an imperfect one, analogous to the effect of gravity on a non-neutrally buoyant sphere, which shifts the equilibria towards one wall in both horizontal and vertical channel orientations (Fox et al. Reference Fox, Schneider and Khair2021). The saddle node (fold) of this imperfect pitchfork defines a critical gradient
$\epsilon _{\textit{crit}}(\textit{Re}_c)$
: for
$\epsilon \lt \epsilon _{\textit{crit}}$
three equilibria coexist (two stable, one unstable), whereas for
$\epsilon \gt \epsilon _{\textit{crit}}$
only one stable equilibrium remains. The case
$\epsilon =-0.5$
(diamonds) is analogous, since, as shown below, negative gradients map exactly onto positive ones. Figure 5(b) shows that
$\epsilon _{\textit{crit}}$
increases monotonically with
$\textit{Re}_c$
: stronger inertia requires a stronger gradient to destroy the upper branch. At
$\textit{Re}_c=250$
(inset),
$\epsilon _{\textit{crit}}\approx 0.38$
, which explains figure 4(b): for
$\epsilon =0.5$
and
$1$
, only one stable equilibrium remains.
(a) Equilibrium position
$Y^{{eq}}$
against
$\textit{Re}_c$
for
$\epsilon =0,\ {0.25,}\ 0.5,\ 1$
and
${-}0.5$
. Filled (open) symbols denote stable (unstable) equilibria. (b) Critical gradient
$\epsilon _{\textit{crit}}(\textit{Re}_c)$
for
$\epsilon \gt 0$
, which exists only for
$\textit{Re}_c\gt \textit{Re}_c^{\textit{crit}}$
. Inset:
$Y^{{eq}}$
against
$ \epsilon$
at
$\textit{Re}_c=250$
, showing saddle-node annihilation of the high-viscosity and central branches at
$\epsilon _{\textit{crit}}\approx 0.38$
. Here
$\kappa =0.2$
.
We attribute the generation of these hydrodynamic transitions to the transverse lift force (
$F_y$
) distributions. Figure 6 plots the corresponding
$F_y$
profiles (scaled by
$\tilde {\rho } \tilde {G}^2 \tilde {R}^4$
) for a sphere held at successive transverse positions across the channel while free to translate streamwise and to rotate. In this representation, zero crossings are equilibria, stable where
$\partial _Y F_y\lt 0$
. At
$\textit{Re}_c=10$
(figure 6
a), the constant-viscosity profile is nearly linear and restoring, crossing zero exactly at
$Y=0.5$
. A positive gradient shifts the whole curve downwards, most strongly near the centreline, moving the stable equilibrium lower; a negative gradient shifts it upwards. At
$\textit{Re}_c=250$
(figure 6
b), the constant-viscosity curve crosses zero three times, reflecting the two off-centre stable equilibria and the unstable centreline. For
$\epsilon =0.25\lt \epsilon _{\textit{crit}}$
, the curve shifts downwards but retains three crossings. For
$\epsilon =0.5$
and
$1$
(
$\gt \epsilon _{\textit{crit}}$
), its upper local maximum falls below
$F_y=0$
, the force-level manifestation of the saddle-node annihilation:
$F_y\lt 0$
throughout the upper half, so any particle released there is driven to the lower equilibrium.
Lift force
$(F_y)$
acting on a sphere held at successive transverse positions
$(Y)$
, at (a)
$\textit{Re}_c=10$
and (b)
$\textit{Re}_c=250$
, with
$\kappa =0.2$
.
Figure 6. Long description
The image contains two line graphs labeled (a) and (b), each showing the lift force (Fy) acting on a sphere held at successive transverse positions (Y). Both graphs have the same x-axis labeled Y, ranging from 0.3 to 0.7, and the same y-axis labeled Fy, ranging from -2.0 to 2.0 in graph (a) and from -0.5 to 0.3 in graph (b). The graphs are divided into two regions: Lower half and Upper half, separated by a dashed vertical line at Y = 0.5. Each graph contains multiple lines representing different values of epsilon (ε): 0, 0.25, 0.5, 1.0, and -0.5. In graph (a), the lift force decreases linearly from positive to negative values as Y increases. In graph (b), the lift force shows a more complex behavior, initially decreasing, then increasing slightly, and finally decreasing again. An inset in graph (b) zooms in on the region around Fy = 0, providing a detailed view of the lift force near this value.
In what follows, we rationalise these transitions with a reduced-order model. We model the transverse motion as an overdamped system,
$\zeta \dot \delta = F_y(\delta ; \textit{Re}_c, \epsilon )$
, with displacement
$\delta \equiv Y - 1/2$
and a dimensionless drag coefficient
$\zeta \gt 0$
. At constant viscosity, the system possesses a natural Couette-flow inversion symmetry: it is invariant under a
$180^\circ$
spatial rotation about the spanwise axis, mapping
$(x,y)\mapsto (-x,1-y)$
, combined with a Galilean shift. This exact symmetry demands that
$F_y$
is strictly odd in
$\delta$
, governed by a supercritical pitchfork normal form:
$\zeta \dot \delta =a_1(\textit{Re}_c)\,\delta -\beta \, \delta ^{3}+\mathcal{O}(\delta ^{5})$
, where
$a_1=\alpha (\textit{Re}_c-\textit{Re}_c^{\textit{crit}})$
and
$\alpha ,\beta \gt 0$
. Physically, for
$\textit{Re}_c\lt \textit{Re}_c^{\textit{crit}}$
(
$a_1\lt 0$
), the wall interactions exert a restoring lift force (Ho & Leal Reference Ho and Leal1974; Vasseur & Cox Reference Vasseur and Cox1976). Increased inertia flattens this lift profile (Anand & Subramanian Reference Anand and Subramanian2023) until the centreline slope reverses at
$\textit{Re}_c^{\textit{crit}}$
, producing two stable off-centre roots,
$\delta ^{{eq}}= \pm \sqrt {a_1/\beta }$
. The viscosity gradient introduces a cross-stream force
$f_v(\textit{Re}_c)\gt 0$
that acts as a symmetry-breaking tilt:
$\zeta \dot \delta =a_1\delta -\beta \delta ^{3}-\epsilon f_v$
. For
$a_1\gt 0$
, this imperfect pitchfork supports three exact roots, given by
$\delta _k = 2\sqrt { {a_1}/{3\beta }}\cos [ ({1}/{3})\arccos (- {3\sqrt {3}\beta ^{1/2}\epsilon f_v}/({2a_1^{3/2}} )) - {2\pi k}/{3} ],$
where
$k=0,1,2$
trace the upper stable, central unstable and lower stable branches, respectively. As the tilt
$\epsilon f_v$
increases, the upper and central roots coalesce at
$\delta ^*=\sqrt {a_1/(3\beta )}$
and mutually annihilate. This saddle-node limit defines the critical gradient:
\begin{equation} \epsilon _{\textit{crit}}(\textit{Re}_c) =\frac {2}{3\sqrt 3} \,\frac {\bigl [\alpha (\textit{Re}_c-\textit{Re}_c^{\textit{crit}})\bigr ]^{3/2}}{\beta ^{1/2}f_v(\textit{Re}_c)} . \end{equation}
Near the onset, where
$f_v$
varies slowly, (4.3) predicts
$\epsilon _{\textit{crit}}\propto (\textit{Re}_c-\textit{Re}_c^{\textit{crit}})^{3/2}$
. The data in figure 5(b) support this scaling: a power-law fit with
$\textit{Re}_c^{\textit{crit}}=151$
yields an exponent of
$1.5\pm 0.1$
, the uncertainty reflecting the numerical accuracy of
$\textit{Re}_c^{\textit{crit}}$
and the choice of fitting range. The behaviour for negative gradients follows from an exact discrete symmetry
$(y,\epsilon ,\textit{Re}_c) \mapsto \bigl (1-y,\,-\epsilon /(1+\epsilon ),\,\textit{Re}_c/(1+\epsilon )\bigr )$
, obtained by a
$180^\circ$
spatial rotation and Galilean shift. A negative-gradient state thus equates to a positive-gradient counterpart at a higher rescaled inertia; for example,
$\epsilon =-0.5$
at
$\textit{Re}_c$
maps onto
$\epsilon =1$
at
$2\textit{Re}_c$
, so the three equilibria observed for
$\epsilon =-0.5$
in figure 5(a) imply
$\epsilon _{\textit{crit}}(2\textit{Re}_c)\gt 1$
, consistent with the monotonic growth of the fold locus. The same reasoning applies to the particle-induced perturbation
$\eta '$
of the viscosity field introduced in § 4.1. In the non-zero-Reynolds-number regime,
$Pe_p=\textit{Re}_p\,\tilde \nu /\tilde D$
is not small (at
$\textit{Re}_c=10$
,
$Pe_p\approx 3$
for thermal transport with
$\tilde D\approx 1.4\times 10^{-7}\text{m}^{2}\,\text{s}^{-1}$
, and
${\approx} 400$
for a small-molecule solute in water with
$\tilde D\approx 10^{-9}\text{ m}^{2}\,\text{s}^{-1}$
), so the scalar is advected around the particle and a quantitative treatment would require the coupled flow-transport problem. Nevertheless, since
$\eta '$
vanishes for
$\epsilon =0$
, respects the symmetries above and enters at
$\mathcal{O}(\epsilon )$
, it can only renormalise the tilt
$f_v$
; because an imperfect pitchfork is the generic unfolding of a symmetric pitchfork under any such perturbation (Golubitsky & Schaeffer Reference Golubitsky and Schaeffer1984), the bifurcation structure is preserved, and only the value of
$\epsilon _{\textit{crit}}$
may change. Finally, to assess the contributions of weak-inertia lifts, we consider a sphere with velocity slip
$\tilde {\boldsymbol{V}}$
, relative spin
$\tilde {\boldsymbol{\varOmega }}_{\textit{rel}}$
and kinematic viscosity
$\tilde \nu$
. It experiences a Saffman lift (Saffman Reference Saffman1965)
$\tilde F_S=6.46\,\tilde \eta _0\tilde R^{2}\tilde V(\tilde G/\tilde \nu )^{1/2}$
directed towards lower-velocity streamlines for a leading particle (
$\tilde V\gt 0$
) and towards higher-velocity streamlines for a lagging one (formally valid for
$\tilde R\tilde V/\tilde \nu \ll \textit{Re}_p^{1/2}\ll 1$
), and a Magnus lift (Rubinow & Keller Reference Rubinow and Keller1961)
$\tilde {\boldsymbol{F}}_M=\pi \tilde \rho \tilde R^{3}\, \tilde {\boldsymbol{\varOmega }}_{\textit{rel}}\times \tilde {\boldsymbol{V}}$
. Since weakly inertial lift is determined by the leading-order velocity field, we evaluate both contributions using the Stokes-limit kinematics,
$V\equiv \tilde V/\tilde G\tilde R=\epsilon \kappa /3$
, obtaining
$F_S=2.15\,\epsilon \,\textit{Re}_c^{-1/2}$
and
$F_M=\pi (\tilde \varOmega _{\textit{rel}}/\tilde G)\, \epsilon \kappa /3$
upon non-dimensionalisation by
$\tilde \rho \tilde G^{2}\tilde R^{4}$
. The relative spin vanishes at
$\mathcal{O}(\epsilon )$
, while the velocity slip is
$\mathcal{O}(\epsilon )$
; hence, the Magnus lift is at most
$\mathcal{O}(\epsilon ^{2})$
and is negligible. The Saffman lift, by contrast, is
$\mathcal{O}(\epsilon )$
, and is therefore comparable in magnitude and identical in sign to the measured tilt (
$F_S\approx 0.68\epsilon$
at
$\textit{Re}_c=10$
, against a centreline shift of
$F_y$
of
${\approx} 0.5\epsilon$
–
$0.6\epsilon$
in figure 6
a), although its formal validity conditions are only marginally satisfied at the lower end of our tested range (
$\textit{Re}_c \approx 10$
,
$\kappa = 0.2$
). We confirm that the trends reported throughout this study qualitatively persist across the tested confinements,
$\kappa \in [0.1,0.25]$
.
5. Conclusions
We have investigated the motion of a rigid sphere in viscosity-stratified plane Couette flow. In the Stokes limit, cross-stream migration is prohibited for all confinement ratios and gradient strengths; the particle attains a constant streamwise velocity slip, set by the competition between the curvature-induced Faxén lag and the gradient–disturbance coupling, and rotates at half the local vorticity. At non-zero Reynolds numbers, the gradient acts as a symmetry-breaking perturbation of the pitchfork bifurcation. Below the pitchfork threshold, it shifts the centreline equilibrium into the less-viscous half, yielding a single off-centre stable equilibrium. Above it, we identify a critical gradient beyond which a saddle-node bifurcation eliminates the high-viscosity equilibrium, producing deterministic focusing in the low-viscosity region regardless of the initial position. The Magnus lift is negligible, whereas a slip-induced Saffman lift, whose estimate is indicative only at the lower end of the range considered (
$\textit{Re}_c\approx 10$
,
$\kappa =0.2$
), has the same sign as, and a magnitude comparable to, the viscosity-gradient-induced tilt; the low-
$\textit{Re}_c$
displacement thus likely reflects both this slip-induced lift and the gradient-induced modification of the wall-induced inertial lift. These findings provide a unified framework for understanding and controlling particle migration in microfluidic flows with naturally occurring or externally imposed viscosity gradients.
Acknowledgements
The authors acknowledge PARAM Pravega HPC system at SERC, IISc.
Funding
This research is supported by ANRF (grant no. ANRF/ARG/2025/003186/ENS).
Declaration of interests
The authors declare no conflict of interest.
η~(x~)
κ=0.2
Rec=250
ϵ=1
κ=0.2
O(ϵ)
Ux
Ωz
Rep=0.01
κ=0.1
Rec=1
Y=0.55
ϵ=0
ϵ=±0.5,1
Rec=10
Rec=250
κ=0.2
Y=0.5
Yeq
Rec
ϵ=0, 0.25, 0.5, 1
−0.5
ϵcrit(Rec)
ϵ>0
Rec>Reccrit
Yeq
ϵ
Rec=250
ϵcrit≈0.38
κ=0.2
(Fy)
(Y)
Rec=10
Rec=250
κ=0.2
Open access
