Conditional stability of particle alignment in
finite-Reynolds-number channel flow

Anupam Gupta, Pascale Magaud, Christine Lafforgue, Micheline Abbas

► To cite this version:

Anupam Gupta, Pascale Magaud, Christine Lafforgue, Micheline Abbas. Conditional stability of particle alignment in finite-Reynolds-number channel flow. Physical Review Fluids, 2018, 3 (11), 10.1103/PhysRevFluids.3.114302. [arXiv:2006.02497] [hal-02145949]

HAL Id: hal-02145949

https://hal.insa-toulouse.fr/hal-02145949

Submitted on 4 Jun 2019

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Conditional stability of particle alignment in finite-Reynolds-number
channel flow

Anupam Gupta,<sup>1,2,*</sup> Pascale Magaud,<sup>3,4</sup> Christine Lafforgue,<sup>5</sup> and Micheline Abbas<sup>1,†</sup>

<sup>1</sup>Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INP, UPS, Toulouse, France

<sup>2</sup>FERMaT, Université de Toulouse, CNRS, Toulouse, France

<sup>3</sup>Institut Clément Ader (ICA), Université de Toulouse, Toulouse, France

<sup>4</sup>Université de Limoges, 33 rue François Mitterrand, 87032 Limoges, France

<sup>5</sup>Laboratoire d'Ingénierie des Systèmes Biologiques et des Procédés, Université de Toulouse,
CNRS, INRA Toulouse, France

<figure><img src="image_2.png" /></figure>

(Received 6 March 2018; published 20 November 2018)

Finite-size neutrally buoyant particles in a channel flow are known to accumulate at specific equilibrium positions or spots in the channel cross-section, if the flow inertia is finite at the particle scale. Experiments in different conduit geometries have shown that while reaching equilibrium locations, particles tend also to align regularly in the streamwise direction. In this paper, the force coupling method was used to numerically investigate the inertia-induced particle alignment, using square-channel geometry. The method was first shown to be suitable to capture the quasisteady lift force that leads to particle cross-streamline migration in channel flow. Then the particle alignment in the flow direction was investigated by calculating the particle relative trajectories as a function of flow inertia and of the ratio between the particle size and channel hydraulic diameter. The flow streamlines were examined around the freely rotating particles at equilibrium, revealing stable small-scale vortices between aligned particles. The streamwise interparticle spacing between aligned particles at equilibrium was calculated and compared to available experimental data in square-channel flow [Gao et al., Microfluid. Nanofluid. 21, 154 (2017)]. The new result highlighted by our numerical simulations is that the interparticle spacing is unconditionally stable only for a limited number of aligned particles in a single train, the threshold number being dependent on the confinement (particle-to-channel size ratio) and on the Reynolds number. For instance, when the particle Reynolds number is ≈1 and the particle-to-channel height size ratio is ≈0.1, the maximum number of stable aligned particles per train is equal to 3. This agrees with statistics realized on the experiments of [Gao et al., Microfluid. Nanofluid. 21, 154 (2017)]. It is argued that when several particles are hydrodynamically connected moving as a unique structure (the train) with a steady streamwise velocity, large-scale hydrodynamic perturbations induced at the train scale prohibit small-scale vortex connection between the leading and second particles, forcing the leading particle to leave the train.

DOI: 10.1103/PhysRevFluids.3.114302

I. INTRODUCTION

The experiments of Segre and Silberberg [1] shed light on the fact that neutrally buoyant particles
experience cross-streamline migration in a parabolic flow if the flow inertia is finite at the particle

<footnote><sup>*</sup>Present address: Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA; anupam.gupta@ensiacet.fr</footnote>

<footnote>†Corresponding author: micheline.abbas@ensiacet.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Scheme showing the possible focusing positions of particles on four spots in the cross-section of a square channel (light blue spheres) and their alignment in the streamwise direction (darker blue). Z indicates the flow direction. X and Y are the directions parallel to the walls. In the simulations particles are placed in the symmetry plane with respect to the X direction.</figcaption></figure>

scale. The dipole interaction with the quadratic part of the flow is responsible of the particle migration. Theoretical computation of the resulting lift force and its dependence on the flow inertia has progressed slowly over decades [2–5]. Understanding this phenomenon opened a new field of applications with the development of microfluidics, where separation or detection of microparticles is operated by hydrodynamic focusing like flow cytometry [6], single-cell encapsulation [7], and cell diagnostics [8]. It can be especially practical in the sense that external fields (like electrical, magnetical) or membranes are avoided.

In more recent experiments, particles were found to accumulate preferentially at equilibrium positions that depend on the conduit cross-section. The accumulation region consists of a ring in a tube flow and of spots at the center of channel faces in square or rectangular ducts as recently reported [9–12]. It has been also observed that, in addition to the existence of equilibrium positions in the cross-section, particles tend to become ordered or evenly spaced in the streamwise direction (so-called trains are formed) [8,13–17]. These observations were obtained in several flow geometries. A sketch of particles assembled in the form of a streamwise train is illustrated in Fig. 1 in the case of square-channel flow.

These particle assemblies originate from the interaction, in shear flow, of particle pairs at finite-flow inertia in the presence of the walls. The experimental observations (usually by optical techniques) of particle trains suggest that at the end of pair interactions, an equilibrium interparticle (streamwise) spacing is reached. This spacing varies like Re_p⁻¹/² (Re_p being the particle Reynolds number defined at the end of the introduction), as it was obtained in tube and later in square-channel flows [13,17]. Neutrally buoyant particles transported by shear flow induce local streamline reversal at finite inertia [18]. As the interparticle spacing in the train structures decreases with the flow inertia, it was first suggested by Matas *et al.* [13] that the train formation is related to the flow induced by one particle in finite-inertia shear, as a particle causes the reversal of streamline direction, but a second particle following such a streamline is cut off from receding by the wall. The two-dimensional (2D) pair dynamics was later investigated by Yan *et al.* [19] in wall-bounded shear (linear) flow. The authors revealed that the particle pair can reach a stable equilibrium or limit cycles at finite inertia, depending on the streamwise boundary conditions.

Nevertheless Lee *et al.* [20] have measured interparticle spacings in channel flow, at different downstream positions of the channel and plotted histograms. Interestingly, the peak in interparticle spacing seemed to continuously shift to larger distances further downstream. The authors noted that this shift becomes noticeable after particles travel long distances of order hundred times the channel height and attributed this to residual viscous repulsive interactions.
We show in this paper that particles assembled in the streamwise direction due to finite-flow inertia reach stable interparticle spacings if a small number of particles is involved. However the apparently long-lived trains become unstable if a large number of particles are aligned, in which case the leading particle leaves the train. The corresponding dynamics seems to be very slow. This observation is made possible by simulating the full dynamics of a few particles aligned along the flow direction on a single spot in the square channel and very long simulation domains to avoid the effect of periodic boundaries in the flow direction or the distant interaction between different trains at different spots. We also show that the maximum size of a stable train depends on the operating conditions that can be gathered under two dimensionless numbers: the particle confinement and the Reynolds number. The particle confinement is defined as the ratio between the particle diameter $d_p$ and the channel hydraulic diameter $H$. The Reynolds number describes the competition between inertial and viscous forces, either at the channel scale, $Re = UH/\nu$ (the so-called channel Reynolds number), or at the particle scale $Re_p = Re * (d_p/H)^2$ (particle Reynolds number). $U$ is the average channel flow velocity and $\nu$ is the kinematic viscosity.

The paper is organized as following. The numerical method is described and validated in Secs. II and III. These two sections are included to assess the relevance of the force coupling method for the simulation of inertia-induced particle migration and alignment. The reader can skip these two sections if not interested in numerical details. In Sec. IV, stable particle assemblies are investigated close to equilibrium. The train formation process and the stable train properties are described as a function of the Reynolds number and of the number of aligned particles. Instability of particle alignment is observed as soon as a large number of particles are aligned in the flow direction. The paper is ended with a discussion (Sec. V) on the possible driving mechanism.

Note that regarding the particle size, for notation convenience and comparison with other theoretical and experimental frameworks, the particle radius $a$ is used in Sec. II, in some places in Sec. III, and in the Appendix. The particle diameter $d_p$ is exclusively used starting from Sec. IV.

## II. NUMERICAL METHOD FORMULATION

The description of the numerical method can be found in Ref. [21]. It is rewritten in this paper before the validation section, for the sake of completeness. Direct numerical simulations of single-phase flows are performed by using the code JADIM for an incompressible Newtonian fluid [22]. The unsteady three-dimensional (3D) Navier-Stokes equations discretized on a staggered grid are integrated in space using the finite volume method. All terms involved in the balance equations are written in a conservative form and are discretized using second-order centered schemes in space. The solution is advanced in time by a second-order semi-implicit Runge-Kutta/Crank Nicholson time-stepping procedure, and incompressibility is achieved by correcting the pressure contribution, which is the solution of the Poisson equation.

Numerical simulations of particle trajectories and suspension flow dynamics are based on multipole expansion of momentum source terms added to the Navier-Stokes equations (namely, force-coupling method as formulated in Refs. [23,24]). Flow equations are dynamically coupled to Lagrangian tracking of particles. The fluid is assumed to fill the entire simulation domain, including the particle volume. The fluid velocity and pressure fields are solutions of continuity Eq. (1) and momentum balance Eqs. (2) and (3):

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \rho \frac{Du}{Dt} = -\nabla p + \mu \nabla^2 \mathbf{u} + f(\mathbf{x}, t), \qquad (2) $$

$$ f_i(\mathbf{x}, t) = \sum_{n=1}^{N_p} F_i^n \Delta[\mathbf{x} - \mathbf{Y}^n(t)] + G_{ij}^n \frac{\partial}{\partial x_j} \Delta'[\mathbf{x} - \mathbf{Y}^n(t)], \qquad (3) $$
where **u** is the fluid velocity, and $\rho$ and $\mu$ are, respectively, the density and dynamic viscosity of the fluid. The body force distribution $f(x, t)$ in the momentum balance Eq. (3) accounts for the presence of particles in the flow. It is written as a multipole expansion truncated after the second term. The first term of the expansion called the monopole represents the force $F^n$ that the particle exerts on the fluid, due to particle inertia, external forcing, or particle-to-particle contact forces [Eq. (4)]. The second term, called dipole, is based on a tensor $G^n$ sum of two contributions: an antisymmetric part is related to external torques applied on the particle, and a symmetric part that accounts for the resistance of a rigid particle to deformation by ensuring zero average strain-rate inside the particle volume, Eq. (5):

$$F^n = (m_p - m_f) \left( g - \frac{dV^n}{dt} \right) + F_{\text{ext}}^n, \quad (4)$$

$$S_{ij}^n(t) = \frac{1}{2} \int \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \Delta'[\mathbf{x} - \mathbf{Y}^n(t)] d^3 x = 0, \quad (5)$$

where $m_p$ and $m_f$ are, respectively, the mass of the particle and that of the fluid in the region occupied by the particle. The finite-size particle is accounted for by spreading the momentum source terms around the particle center $\mathbf{Y}^n$ using a Gaussian spherical envelope, one for the monopole $\Delta(\mathbf{x}) = (2\pi\sigma^2)^{-3/2}e^{(-|\mathbf{x}|/2\sigma^2)}$ and another one for the dipole $\Delta'(\mathbf{x}) = (2\pi\sigma'^2)^{-3/2}e^{(-|\mathbf{x}|/2\sigma'^2)}$. The widths of the Gaussian envelopes, $\sigma$ and $\sigma'$, are set with respect to the particle radius $a$ such that the settling velocity and the hydrodynamic perturbation generated by a particle in a shear flow are both exactly matched to Stokes solutions [$\sigma = a/\sqrt{\pi}$ and $\sigma' = a/(6\sqrt{\pi})^{1/3}$] for a single particle.

The particle translation and rotation velocities are obtained from a local weighted average of the volumetric fluid velocity (respectively, rotational velocity) field over the region occupied by the particle [Eqs. (6) and (7)]:

$$\mathbf{V}^n(t) = \int \mathbf{u}(\mathbf{x}, t) \Delta[\mathbf{x} - \mathbf{Y}^n(t)] d^3 \mathbf{x}, \quad (6)$$

$$\Omega^n(t) = \frac{1}{2} \int [\nabla \times \mathbf{u}(\mathbf{x}, t)] \Delta'[\mathbf{x} - \mathbf{Y}^n(t)] d^3 \mathbf{x}. \quad (7)$$

Particle trajectories are then obtained from numerical integration of the equation of motion as in Eq. (8):

$$\frac{d\mathbf{Y}^n}{dt} = \mathbf{V}^n. \quad (8)$$

This modeling approach allows calculating the hydrodynamic interactions with a moderate computational cost. For a neutrally buoyant particle, the monopole and the antisymmetric contribution to the dipole are strictly zero. Only the symmetric part of the dipole (Stresslet) allows us to account for the interaction between the particle and the shear flow. Eight grid points per particle diameter are usually sufficient to correctly capture this interaction.

The method has been validated in the limit of vanishing particle Reynolds number [23,24]. It has later been extended to the case of finite-flow inertia at the particle scale, i.e., $\text{Re}_p = O(1)$ [25,26]. Loisel et al. [15] have shown that the Stresslet components of a single particle placed in a linear flow compare very well with DNS measurements, up to particle Reynolds number equal to 5 [27]. Additional validation tests with a single particle in quadratic flow are presented in the next section and in the Appendix.

As for the interaction between two spheres in a linear flow, Yeo and Maxey [25] have shown that the FCM gives the right relative particle trajectories at $\text{Re}_p = O(1)$. When two particles are initially placed in the shear plane (perpendicular to the vorticity direction), their relative trajectory remains in-plane, and it is open or reversed depending on the initial shift in the shear direction ($\delta y$),
<figure><img src="image_4.png" /><figcaption>FIG. 2. Slip velocity (left panel) and migration velocity (right panel) vs. particle position $Y/H$ for $Re = 13$ (triangles) and $Re = 39$ (squares) with particle size $d_p/H = 0.06$. The solid black curve shows the law proposed by Ref. [31]. The vertical dotted line indicates a distance from the wall equal to the particle radius. The horizontal dashed line in left panel represents the Faxen correction.</figcaption></figure>

of the lagging particle with respect to the leading one. The bifurcation between the two types of trajectories is close to the one found in LBM simulations [28]. The off-plane spiraling interaction is less well captured when the gap between particle surfaces is smaller than $0.1d_p$; however, the amplitude of the relative velocity is very small in that case, and it does not play a significant role in the system studied in this paper (particles in the same shear plane).

### III. VALIDATION OF THE NUMERICAL METHOD

At very low Reynolds number, a small neutrally buoyant spherical particle follows the flow streamlines. Near a wall, both the translational and rotational particle velocities are smaller than the local fluid flow velocities [29]. However, the particle does not experience a wall-normal motion for reversibility reasons. If the flow is slightly inertial at the particle scale, the neutrally buoyant particle experiences lift perpendicular to the flow streamlines, in the presence of shear, the intensity and direction of the lift depending on the flow configuration, and on whether the particle is free to rotate or not. In channel flow, the interaction of the particle Stresslet with the curved background flow profile induces a lift force oriented toward the channel walls when the particle is located near the central region [2]. This force is enhanced by flow inertia. When the particle is very close to the wall, the particle slip is large. The particle slip in the presence of shear near a wall leads to a lift force oriented toward the high velocity region (as computed for instance by Cherukat and Mclaughlin [30]). Hence, there is an equilibrium position, between the flow center and the walls, where the particle is transported force-free. The equilibrium position is closer to the channel walls when the flow inertia increases as it was demonstrated theoretically in channel flow, first by Ref. [4] up to $Re = O(100)$ and later by Ref. [5] up to $Re = O(1000)$, assuming pointlike particles.

The validation tests shown here were realized in square-channel flow. Periodic boundary conditions were used in the flow direction (Z) and no slip at the walls (in X and Y directions). The ratio of the particle diameter to channel height was $d_p/H = 0.06$ and $0.11$. The channel length in the streamwise direction was equal to $28.8d_p$, where $d_p$ is the particle diameter. The grid distribution was set to ensure 8 grid points per particle diameter. The fluid flow was initially set to the steady solution of square-channel flow, and a constant pressure gradient was applied in the z direction. The particle was seeded at different Y locations in the midplane ($X = H/2$).

#### A. Particle freely moving in square-channel flow

In the first test, the particle was moving freely during approximately $10a^2/v$, before its streamwise and wall-normal velocities were recorded (a is the particle radius). The streamwise slip and wall-normal particle velocities are shown in Fig. 2. The two velocity components are compared,
at channel Reynolds number $Re = 13$ and 39, to theoretical expressions for a pointlike particle in 2D Poiseuille flow (see the summary on this in Asmolov *et al.* [31]). For the smallest particle size, the effect of the flow three-dimensionality on the particle motion is expected to be relatively small. The slip velocity is normalized by $aG_m$, where the shear rate $G_m = 4U_m/H$ is calculated from the maximum velocity in the channel center $U_m$. The particle slip is not impacted by the flow inertia in this range of Reynolds numbers. The agreement with the theoretical velocity as derived from Goldman *et al.* [29] is acceptable near the wall. However, the slip does not vanish in the channel center because of the flow curvature (this effect on the particle motion was written formally in Faxén laws at low Reynolds numbers). The slip magnitude is roughly twice of the Faxén correction $4U_m a^2/(3H^2)$. The same observations were reported in the studies of Loisel *et al.* [15] and Asmolov *et al.* [31], realized in 2D Poiseuille flow, using numerical simulations based on FCM and Lattice Boltzmann, respectively. As for the migration velocity (scaled by $aG_m$ and $Re_p/6\pi$), its trend is similar to the prediction based on a pointlike particle at low but finite Reynolds number [3]. The shift to the left of the numerical points is a joint consequence of the flow being 3D instead of 2D and to the underestimation of the hydrodynamic interaction between the particle and the wall. We note that this underestimation seems to be only effective on the wall-normal direction and not on the slip parallel to the wall.

## B. Lift force computation

The force coupling formulation allows solving the mobility problem, i.e., the particle is displaced and rotated under a given forcing. A neutrally buoyant particle is not subject to any external forcing. The direct calculation of the force pushing the particle to move across the flow streamlines is not possible, because unlike other particle-resolved methods, the force coupling method does not guarantee the no-slip boundary condition on the particle surface, and therefore the surface traction cannot be directly calculated. Instead, an iterative algorithm is set to compute the wall-normal force that should be applied on a particle placed in a shear flow, to prohibit the particle motion in the wall-normal direction.

After recording the velocities of the freely moving particle in the previous test, a force was then applied on the particle, only in the wall-normal direction, in a way that ensures zero migration velocity. This force was applied to the particle motion through $F_{ext}^n$ in Eq. (4). Its initial value was set to zero. The force was then updated at the iteration $k$ from the value at iteration $k-1$, according to a penalty method: $F_{ext}(k) = F_{ext}(k-1) - \lambda[6\pi \mu a V(k-1)]$. The iterations were stopped when $V(k)$ became very close to zero. $\lambda$ is an arbitrary constant, which should be chosen not very low to reduce the time needed for convergence and not very high to avoid numerical instability. Note that the computation of this force was first realized in the case of a particle placed near a wall in a linear flow. The details of this test are written in the Appendix. The similarity between this force applied to prohibit the particle wall-normal motion and the theoretical predictions of the quasi-steady lift force, led us to call it “lift force” in this paper.

Figure 3 shows the evolution of the lift force calculated $F_l/(\rho U_m^2 a^4/H^2)$ as a function of the particle position, in the midplane ($X = H/2$) of the square-channel flow. The particle radius $a$ is used in the force scaling. The negative sign indicates a force pushing the particle away from the wall. The numerical results, at different particle diameters and channel Reynolds numbers, are compared to the theoretical work of Hood *et al.* [32]. Their work was developed in square-channel flow geometry, assuming that the wall falls in the inner layer perturbed by the particle (weak inertial stresses compared to viscous stresses). The inertial lift force in the x and y directions was shown to depend on the particle radius in the form $F_l/(\rho U_m^2 a^4/H^2) = C_4 + C_5 a/H$, where $C_4$ and $C_5$ are constants that depend only on the location of the particle. For the lowest Reynolds number ($Re = 13$) and smallest particle size, the numerical force obtained by the FCM is in good agreement with the profile established in Ref. [32]. Scaling the force by $\rho U_m^2 a^4/H^2$ is consistent in the channel center, but not near the channel wall. Still at $Re = 13$, the dimensionless lift is lower when larger particle size is used. Note that, when the flow conditions are unchanged, the larger the particle the stronger
<figure><img src="image_5.png" /></figure>

FIG. 3. Lift force scaled by $ρU_m^2 a^4 / H^2$, acting on a particle in a square-channel flow versus the position of the particle in the y direction (at $x = H/2$) for different Re. The red and blue symbols are the lift force from our simulation for particle diameter $d_p/H = 0.06$ and 0.11, respectively. The symbols are for $Re = 120$ (circles), 38 (squares), and 13 (triangles). The corresponding $Re_p$ is between 0.05 and 1.5. The red dashed and blue solid lines are obtained from the solution of Ref. [32] for particle diameters $d_p/H = 0.06$ and 0.11, respectively.

is the impact of truncating higher order terms in the multipole expansion (quadrupole, sextupole...). When the Reynolds number is increased, the force calculated by the FCM deviates with respect to this scaling (it becomes lower), in a way coherent with the theoretical analysis based on matched asymptotic expansions [4,5].

## IV. PARTICLE ALIGNMENT IN CHANNEL FLOW

Here and in the following sections, the particles are freely transported by a square-channel flow, unless otherwise stated. We considered particularly the square-channel flow configuration to discuss the stability of particle alignment, because in the conditions of this paper, the equilibrium positions are stable and well established (at the midplane of the four channel walls) and trains have been well characterized by the experiments of Gao et al. [17]. The trajectory of a single particle migrating toward an equilibrium spot (as sketched in Fig. 1), using the same numerical method, can be found in Ref. [33]. When particles are randomly seeded in the simulation domain, they experience first a lateral motion, i.e., perpendicular to the velocity isocontours, then cross-lateral migration, i.e., parallel to the closest wall, till they reach equilibrium positions [33]. Both processes are slow, the former stage being faster than the following one. The establishment length scale of particle migration is quite large [$O(1000H)$ at $Re = O(100)$], and the lower the Reynolds number the larger the establishment length is. To focus on the streamwise ordering process, the lateral and cross-lateral migration stages are bypassed by initially placing the particles near their equilibrium spots (in the midplane $x = H/2$). During the simulations, the particles were observed to remain in this symmetry plane.

The operating conditions consist of a channel Reynolds number $O(100)$, a particle confinement $d_p/H$ in the range 0.077–0.14, and a solid volume fraction less than 1%. The flow was resolved using a uniform mesh grid with 78 × 78 grid points in the square cross-section, to ensure 8 grid points per particle diameter for a reasonable numerical accuracy. We carefully verified that the box length $L$ does not impact the results shown here ($29 \le L/d_p \le 60$, $d_p$ being the particle diameter).

The streamlines around a single particle are first shown, since they contribute to the alignment process. Then particle relative trajectories are used to show stable assemblies when small number
<figure><img src="image_6.png" /><figcaption>FIG. 4. Zoom on the flow streamlines in the (YZ) frame of a particle at equilibrium in a square-channel flow. The black region illustrates the particle (stretched for convenience). Re<sub>p</sub> = 0.5 (left), 1.5 (center), and 3.0 (right). The two stagnation points (indicated with crosses) get closer to the particle surface as Re<sub>p</sub> increases. The fore-aft asymmetry increases with Re<sub>p</sub>. Red arrows show the flow direction.</figcaption></figure>

of particles is involved (two or three) and unstable assemblies as soon as the number of particles
becomes “large.”

A. Streamlines around a single particle at equilibrium

Figure 4 shows the flow streamlines in the (YZ) frame attached to a single particle located at an equilibrium spot, for different Reynolds numbers at particle confinement $d_p/H = 0.11$ (the images are stretched and zoomed for format convenience). The saddle points both in front of and behind the particle take place in the presence of shear as soon as the flow inertia is finite at the particle scale [34]. Whereas the reverse streamlines induced by a particle in a channel flow are open under Stokes flow condition, they are of spiralling nature at finite inertia [20]. The stable spirals act as attractor regions. The centers of these forward and backward spirals are closer to the particle surface as Re<sub>p</sub> increases and the size of this attractive region becomes wider. The two horizontal lines added to these figures show that the gap (in the y direction) between the forward and backward attractors also increases with Re<sub>p</sub>, in relation with the symmetry breakup at finite-flow inertia.

B. Stable particle assembly

When two particles are found close to each other near an equilibrium spot [as illustrated in Fig. (1)], they are located in-plane, i.e., in the plane parallel to the flow and to one wall-normal direction. Studies dedicated to the interaction of a particle pair in a linear flow at finite-flow inertia [19,28], show that the relative trajectory of the lagging particle with respect to the leading one is open (respectively, reversed) when the shift $\delta Y$ in the position of the particle centers in the wall-normal direction is large (respectively, small).

Figure 5 shows the relative trajectory of a leading particle (particle 2) with respect to the lagging one (particle 1) in the square-channel flow at $d_p/H = 0.11$. The trajectory of particle 2 is trapped in a basin of attraction, nearby the forward attractor of particle 1, following a spiralling motion. At $Re_p = 0.5$, the trajectory of particle 2 with respect to particle 1 is similar to the streamlines around the freely rotating particle 1 (when isolated). The first part of the relative trajectory (in Fig. 5) is of reversing nature since $\delta Y = Y_{P2} - Y_{P1}$ is initially small and negative. However, after reversing, particle 2 does not go off to infinity. The inertia-induced lift induces cross-streamline relative motion, which is coupled to the trajectory reversal leading to an equilibrium spacing between the particles. This type of interaction of two finite $Re_p$ spheres, is an additional aspect compared to the open and reversed trajectories in linear flow. It involves the quadratic nature of the flow and the proximity of the particles to the wall. When $Re_p$ increases, particle 2 converges faster toward equilibrium, the convergence pathway depending on the initial position of particle 2 with respect to particle 1. Figure 5(b) ($Re_p = 1.5$) shows that if the relative position is chosen carefully, the
<figure><img src="image_7.png" /></figure>

FIG. 5. Panels (a), (b), and (c) contain the pair-dynamics in a two-particle train for Re<sub>p</sub> = 0.5, 1.5, and 3.0, respectively. The trajectories of the leading particle with respect to the lagging one are overlaid on the streamlines around a single particle. The initial position of the center of the leading particle is shown with asterisks. Red arrows show the flow direction. (d) Flow streamlines (in the frame of leading particle) around a stable pair of particles obtained for Re<sub>p</sub> = 3, showing a stable small recirculation connecting both particles.

leading particle converges toward the basin of attraction even if the initial distance is as far as 9d<sub>p</sub>.
At equilibrium, the vortex in front of the lagging particle and the one behind the leading particle
connect to form one close vortex as shown is Fig. 5(d).

We realized the same type of simulations with three particles. Figure 6(a) shows the particle
trains at equilibrium for three different Re<sub>p</sub> at d<sub>p</sub>/H = 0.11. The streamwise position of the lagging
particle is set arbitrarily to zero for this figure. The interparticle spacing and average train distance
from the wall both decrease when Re<sub>p</sub> is increased. These trains are not perfectly aligned in the
flow direction, but they are relatively inclined as shown in Fig. 6(b). This has been also observed
experimentally by Matas et al. [13] at higher Re<sub>p</sub> in a tube flow, while the inclination is absent at
smaller Re<sub>p</sub>. The evolution of the train inclination with the Reynolds number is coherent with the
gap between forward and backward stagnation points shown in Fig. 4, which increases with Re<sub>p</sub>.
In addition, it can be noted that the spacing between the leading and second particle is 10% greater
than between the second and third (lagging) one in all cases.

The distance between the train barycenter and the closest wall $y_T$, as well as the average
interparticle distance at equilibrium $l$ are plotted in Fig. 7. Most of the simulations were realized
with square-channel flow and $d_p/H = 0.11$. When $Re_p$ is increased, the train gets closer to the
channel wall (similarly to the single particle) and the average interparticle distance decreases. The
decrease of the average distance with $Re_p$, observed similarly in the experiments [13,16,17], is
consistent with the fact that the basin of attraction is closer to the particle surface, when the particle
Reynolds number is increased (Fig. 4). Figure 7 contains also information on particle assembly
when the number of particles per train is increased. At a given $Re_p$, the train gets slightly closer
to the wall when the number of particles per train increases. The average interparticle distance
seems to slightly decrease as well. The train statistics are compared to the experimental ones of
Gao et al. [17], which were realized in similar conditions. The trend of the trains statistics with
<figure><img src="image_8.png" /><figcaption>FIG. 6. (a) Steady configuration of three-particle trains at different Reynolds numbers in the (YZ) plane. Particles in the same train have the same color and are linked with lines for eye guidance. Re<sub>p</sub> = 0.5 (red), 1.5 (blue), and 3 (black). The z coordinate of the lagging particle in a train is set arbitrarily to 0. The velocity profile is represented on the right of the figure. (b) Plots representing the center positions of the three-particle trains at different Reynolds numbers [same colors as in (a)].</figcaption></figure>

respect to Re<sub>p</sub> is similar. There is a uniform shift between experimental and numerical results. The discrepancy of train positions at equilibrium is almost suppressed when the mesh resolution is twice finer. However, the shift persists for the average interparticle spacing. This issue deserves to be further delved into in the future if precise information on the stability of particle alignment is needed. It requires the computation of the interaction between several particles in channel flow, at identical operating conditions, using different available (at least numerical) methods.

### C. Unstable assembly

When a fourth particle is seeded close to the three-particle train of Fig. 6 (Re<sub>p</sub> = 1.5), either in front of or behind the train, the front particle is lifted up, leaving a stable three-particle train behind it. The relative trajectory of a fourth particle placed in front of a three-particle stable train is shown in Fig. 8(a). Several initial configurations led to the same result. Even if the leading particle tends to follow the reversed streamlines in a first stage, it does not converge toward the attractor. The departure of the front particle can be found in Fig. 8(b) from the evolution of the spacing between the front and second particles (in black lines). A video sequence of particle positions in the channel for this case is shown in the supplementary material [35]. The particle that leaves the train reaches an equilibrium position y<sub>P</sub> located slightly further from the wall than the position y<sub>T</sub> of the train left
<figure><img src="image_9.png" /></figure>

FIG. 7. Properties of stable particle trains. (a) shows the distance between the train barycentre and the closest wall $y_T$. (b) shows the average interparticle distance $l$. Red Squares, upward-pointing triangles, and black plus symbols are for 2-, 3-, and 4-particle assembly in square channel with $\frac{d_p}{H} = 0.11$. Blue downward-pointing triangles and black circles are for 3- and 4-particle assembly in square channels with $\frac{d_p}{H} = 0.14$ and 0.08, respectively. Black cross is for 4-particle in rectangular channel ($\frac{d_p}{W} = 0.09$). Black stars are from the experiments of [17] realized in square channel flow with $\frac{d_p}{H} = 0.11$. The magenta square and cyan triangle are obtained for 2- and 3-particle trains, using a twice finer numerical resolution in a square channel with $d_p/H = 0.11$.

behind and is thus slightly faster. Since we used periodic boundary conditions in the flow direction, the particle that leaves the train from the front joins it from behind, the new leading particle leaves in turn, and so on. The same observations were noted for trains with the larger number of particles.

Note that this instability starts to be discernible after the particles travel a long distance downstream (i.e., ≈30–40*H* in the case of Fig. 8). It is probably for this reason that [20] have observed a shift of the distribution of interparticle distances toward larger values in the measurements at different distances from the channel inlet, without observing any change in interparticle spacings within images of dimensions of *O*(10*H*). To detect the eventual departure
<figure><img src="image_10.png" /><figcaption>FIG. 8. (a) Trajectories of the fourth particle placed in front of a stable three-particle train at Re<sub>p</sub> = 1.5 and d<sub>p</sub>/H = 0.11 for different initial positions. These trajectories are relative to the (new) leading particle of the remaining stable three-particle train (at Re<sub>p</sub> = 1.5). The trajectories are overlaid on the streamlines, in (YZ) plane. Panel (b) contains the evolution of the relative particle spacing in the train (ΔZ). Black, red, and blue lines correspond to the distance between front-second, second-third, and third-fourth. The Z coordinate corresponds to the streamwise position of the front, second, and third particles for each curve, respectively.</figcaption></figure>

of the leading particle from the train structure by optical techniques, it would be required to use either two synchronized cameras at different streamwise positions or a camera frame following a long-lived train that contains a large number of particles.

# V. DISCUSSION ON THE DESTABILIZING MECHANISM

The maximum number of stable particles in a train can be tuned by changing the particle confinement and/or the fluid inertia, as summarized in Fig. 9. Note that to change d<sub>p</sub>/H while keeping constant the particle Reynolds number, the channel Reynolds number should be changed accordingly. It is clear from Fig. 9 that both Re<sub>p</sub> and confinement play an important role. For rectangular channel cross-section, the confinement is defined as the ratio between the particle diameter and the channel height or width, whichever is smaller (setting the largest velocity gradients). The number of particles stable in a train increases when the flow inertia is increased and when the particle size is decreased. It is striking to note that the maximum length of stable train structures is approximately equal to the channel hydraulic diameter (H for square channels).

The numerical results encouraged us to revisit the statistics of Gao et al. [17] realized on experiments with particle size d<sub>p</sub>/H = 0.11 in square-channel flow. The histograms of the number of aligned particles in a single train exhibit a very sharp peak at N<sub>p</sub> = 3, for Re<sub>p</sub> ranging from 0.1 to 3 and for particle concentration between 0.02 and 1%. The percentage of trains constituted
<figure><img src="image_11.png" /><figcaption>FIG. 9. A diagram showing the maximum number of particles in a stable train for different particle Reynolds number $Re_p$ and confinement $d_p/H$ in a square-channel flow, as predicted from the FCM simulations.</figcaption></figure>

of three particles is shown in Fig. 10 as a function of the suspension concentration (defined as the solid volumetric fraction). This figure shows that most of the trains are constituted of three aligned particles at low concentration ($\phi = 0.08\%$). A longer train detected in the experiments (with imaging techniques [17]) as such, might be the result of a transient alignment. As the concentration increases, the percentage of three-particle trains decreases but remains dominant. The concentration has the dual effect of both increasing the number of particles available for alignment and the dispersive hydrodynamic interactions between them.

From these results, it seems that the particle assembly under finite inertia is a weakly coupled system. The origin of particle alignment seems to result from a favorable vortex connection between consecutive particles as illustrated in Fig. 5. The vortex generated behind the front particle interacts with the vortex induced in front of the lagging particle, minimizing by that the fluctuating kinetic energy (in a similar way to particles interacting in oscillatory fluid flows [36], as discussed by Ref. [37]). However, the vortex connection does not seem to occur when the train exceeds a

<figure><img src="image_12.png" /><figcaption>FIG. 10. Probability of finding trains with three aligned particles, obtained from post-processing the experiments of Ref. [17], versus the suspension concentration. The probability is the highest at low concentration. Different symbols are for different $Re_p$: triangles, stars, and circles are for $Re_p = 0.12$, 0.7, and 1.8, respectively. These results are for $d_p/H = 0.11$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Plots of velocity streamlines relative to the front particle in the microchannel in the XZ plane for (a) a single particle, (b) two-particle train, and (c) three-particle train at $Re_p = 1.5$. Colored contours show the streamwise velocity perturbation (subtracting the unperturbed channel flow from the instantaneous velocity field) scaled by the average flow velocity. Flow direction is from left to right.</figcaption></figure>

certain number of aligned particles. Velocity perturbation induced by the individual particle Stresslet in bounded shear flow decays as $1/r^2$ at a distance from the particle center $r < H$. Since the train morphology does not fundamentally change when the number of particles increases, i.e., interparticle spacing does not decrease significantly when the train length increases, hydrodynamic repulsion between pairs cannot be the driving mechanism. Visualizations of the flow perturbation at the channel scale (Fig. 11) reveal that the assembled particles move like a unique coherent structure, with a perturbed outer region that expands as the number of aligned particles increases. Sequences of snapshots for the velocity perturbations reveal that the vortex connection starts first between the lagging and the second to last particle, etc., until reaching the front of the train. It is likely that the hydrodynamic perturbation induced by the large structure when its length reaches the channel scale, prohibits the vortex connection between the leading particle and the second one, pushing the front particle to move forward.

Nevertheless, the observed departure of the leading particle, for instance in the simulations corresponding to Fig. 8, does not depend drastically on the accuracy of the computed hydrodynamic perturbation induced by each particle. Since the Stresslet terms are the essential ingredients to capture the hydrodynamic interaction of the neutrally buoyant particles with the shear flow, we tuned the Stresslets terms to test whether this has an impact on the train stability. When we realized simulations with constant-value Stresslets (obtained from the converged three-particle train), instead of updating the Stresslets to maintain the zero-average strain rate inside the particle volume [Eq. (5)], very similar relative trajectories were observed. This suggests that the conditional stability described in this paper is a robust phenomenon that does not exclusively depend on the accuracy of the interactions at small (particle) scales.

A stability analysis, hard to realize on this system, would probably help to rationalize the effect of increasing the particle number on the stability of particle alignment. Here we limit our argument to the energy budget. At a given Reynolds number, when the number of particles assembled in a train increases, the slip velocity between the train and the ambient fluid flow increases. Although the neutrally buoyant particles move force-free in the flow, the slip velocity induces energy dissipation. The ratio between the energy dissipated by the train and the energy of the flow pushing forward the leading particle is plotted in Fig. 12. This ratio increases with the particle number and confinement, and it decreases with the Reynolds number. Figure 12 suggests that the particle assembly should not cost above a threshold (around 2.5% of the flow energy on the particle scale) for the system to remain stable.

## VI. CONCLUDING REMARKS

After the validation of the numerical method, we gave some insight on the dynamics of a pair of neutrally buoyant particles that tend to align in the streamwise direction and form trains in channel
<figure><img src="image_14.png" /><figcaption>FIG. 12. The ratio between the drag force on particle train and forward pressure force applied by the fluid on the front particle. The symbol legend is idem to Fig. 7.</figcaption></figure>

flows. All the results were obtained after the inertial migration stages were accomplished, where the particles were located close to stable equilibrium spots and in one shear plane. The simulations were realized in domains long enough to eliminate any influence from the streamwise periodic boundaries. Trains of particles revealed to be slightly closer to the channel walls than a single particle at equilibrium, and therefore they had a slower streamwise velocity. The trains were slightly inclined with respect to the flow direction (lifted forward) when the Reynolds number increases, as already observed in experiments realized in tube flow [13]. The trains were unconditionally stable only in a limited range of Reynolds numbers and particle diameter-to-channel height ratios. When the train length increases, the hydrodynamic perturbation induced by the train structure, is likely stronger than the perturbation induced by the Stresslets at the individual particle scale, pushing.

Our numerical results, obtained using a truncated multipole expansion, agree qualitatively well with the experiments [17] realized within the same range of operating conditions. Future investigation on the interaction between one and several pairs of spheres near channel walls are required to assess quantitatively the bifurcation between stable and unstable alignment. Moreover, the conclusions on the particle assembly are valid for moderate particle size and when the solid volumetric concentration is low. When the particle size is almost half of the channel height, additional sets of equilibrium positions take place alternating on opposite walls, like in Refs. [16,38]. This situation could not be examined by the Force Coupling Method as implemented here, mainly because of the truncation of the multipole expansion used in Eqs. (1)–(3). When the suspension volumetric concentration is not negligible ($\phi \gtrsim 0.5\%$) hydrodynamic dispersion is expected to decrease the alignment potentiality in a way complex to predict.

## ACKNOWLEDGMENTS

This work was granted access to the HPC resources of the supercomputer center CALMIP under Allocation No. 2017-P1002 and of GENCI under the Allocations No. x20162b6942 and No. A0012B06942. The authors would like to thank Dr. Y. Gao for revisiting the statistics in his experimental measurements, and his contribution to Fig. 10.
<figure><img src="image_15.png" /></figure>

FIG. 13. Left Panel: Numerical calculation of the wall-normal force (blue open circles) applied on a particle near a wall in linear flow, and comparison with different theoretical predictions. The force is scaled by $\mu a V_{\text{slip}}$. The different theories are from Saffman [39] (blue asterisks), Cherukat and Mclaughlin [30] (red upward-pointing triangle), Lovalenti in the Appendix of Ref. [30] (green square) and Magnaudet [40] (magenta right-pointing triangle). Rigth Panel: streamwise particle slip velocity versus $Re_p$. The red circles are from FCM simulations. The lines are the prediction of the slip velocity from Goldman *et al.* [29] at $y_p = 1.54$ and 2.35a, respectively.

APPENDIX: PARTICLE NEAR A WALL IN A LINEAR FLOW

The validation of the wall-normal force calculation was realized first by placing a neutrally-buoyant particle of radius $a$ near the bottom wall of a plane Couette flow. The domain size was $10.6a$ in the flow $x$ and wall-normal $y$ directions, and $8.1a$ in the spanwise $z$ direction. The computational grid was uniform, and the mesh size was $\Delta x = a/4$. Periodic boundary conditions were used in the $x$ and $z$ directions. The bottom wall was stationary and the top wall moving with $V_W = \dot{\gamma}H$, where $H$ is the distance between the walls and $\dot{\gamma}$ is the shear rate. The particle was placed at a given position, and iterations were realized to find the force required to prohibit particle motion in the wall-normal direction.

An example of the results is shown in Fig. 13(a) obtained by placing the particle at $y_{P0} = 1.66a$ near the bottom wall. This figure shows the increase of the wall-normal force with the Reynolds number defined in linear flow as $Re_p = \dot{\gamma}a^2/\nu$. The force is scaled by the viscous drag, $\mu a V_{\text{slip}}$, where $V_{\text{slip}}$ is the particle slip velocity (in the streamwise direction) with respect to the unperturbed local fluid flow, and $\mu$ is the dynamic fluid viscosity. Note that $V_{\text{slip}}$ is not known *a priori* but calculated from the simulation result at equilibrium, upon completion of the iterative procedure used to obtain the force. Figure 13(b) shows $V_{\text{slip}}$ scaled by the wall velocity as a function of $Re_p$. It is compared in the same graph with the Stokes flow limit from [29] at the two wall-normal positions $y_p = 1.54$ and $2.35a$ reported in Table 2 of their paper. Note that in our simulations the particle position $y_p$, initially equal to $1.5a$ in the simulations at different $Re_p$ is found between $1.66a$ and $1.7a$ at the end of the iterative procedure. The calculated slip is close to the Stokes flow prediction. Its amplitude decreases slightly with the Reynolds number. This is not an inertial effect. It is rather related to the fact that the steady particle position is not the same.

Several theoretical works allowed determining the lift force applied on a particle in a linear flow assuming no fluid acceleration at the particle scale. In this paper, we compare the numerical results with that of Refs. [30,39,40]. The expressions of the lift force obtained by the different works are listed below. All of them are scaled with $\mu a V_{\text{slip}}$ and take into account the proximity of the particle to the wall, except the expression of Saffman [39] obtained in unbounded shear flow. Figure 13 shows that the numerical results agree very well with the theories that take into account the wall presence.

The expressions of the lift force resulting on a finite-size particle in a linear flow are given in this Appendix, from different sources in the literature. The first one does not take into account the presence of the wall. The last three contain the parameter $\kappa$ which is the ratio between the particle position with respect to the closest wall and the particle radius:
(1) Eq. (3.11) of Saffman 1965 [39], with a correction of $\frac{1}{4\pi}$ published in the erratum:

$$F_l = \frac{81.2}{2\pi} a \sqrt{\dot{\gamma}/\nu} = 6.46 \sqrt{\text{Re}_p}. \quad (A1)$$

(2) Eq. (4.2) of Cherukat and McLaughlin 1994 [30]. This work also accounts for the distance between wall and particle; $\kappa = a/\nu_{p0}$ and valid in the regime $\nu_{p0} \ll \min(L_s, L_\gamma)$, where Stokes length $L_s = \nu/V_{slip}$, Saffman length $L_\gamma = (\nu/\dot{\gamma})^{1/2}$ and dimensionless parameter $\Lambda_\gamma = \dot{\gamma}a/V_{slip}$:

$$F_l = \text{Re}_s I = \text{Re}_p (A \Lambda_\gamma^{-1} - B + C \Lambda_\gamma), \quad (A2)$$

where $\text{Re}_s = aV_{slip}/\nu$; $A = 1.7631 + 0.3561\kappa - 1.1837\kappa^2 + 0.84516\kappa^3$, $B = 3.24139/\kappa + 2.676 + 0.8248\kappa - 0.4616\kappa^2$, and $C = 1.8081 + 0.879585\kappa - 1.9009\kappa^2 + 0.98149\kappa^3$.

(3) Eq. (A 14) of Cherukat and McLaughlin 1994 [30], which is derived by Lovalenti. We nondimensionalized Eq. (A 14) by $(\rho v a V_{slip})$:

$$F_l = \frac{18\pi}{32}\text{Re}_s + \frac{7}{8}\pi\text{Re}_p - \frac{66\pi}{64}\text{Re}_p\left(\frac{1}{\kappa} + \frac{27}{16}\right)\left(\frac{6\pi 55}{144 \times 4}\right)\text{Re}_p\left(\frac{a\dot{\gamma}}{V_{slip}}\right) \\ = \pi\text{Re}_p\left[\frac{18}{32}\frac{1}{\Lambda_\gamma} + \frac{7}{8} - \frac{66}{64}\left(\frac{1}{\kappa} + \frac{27}{16}\right)\left(\frac{6 \times 55}{144 \times 4}\right)\Lambda_\gamma\right]$$

(4) Eq. (17) of Magnaudet 2003 [40], which is valid for spherical bubble, droplet and particle. It is valid for the object close to the wall. For particles, in case of viscosity ratio $\lambda \to \infty$, the nondimensionalized lift force is

$$F_l = \frac{\pi}{4}\text{Re}_s(3/2)^2 \left[ \left(1 + 11/8\Lambda_\gamma^2 5/3\right) - \frac{11}{6}\Lambda_\gamma(\kappa^{-1} + 0.84) \right] \\ = \frac{9\pi}{16}\text{Re}_p \left[ \left(\Lambda_\gamma^{-1} + 55/54\Lambda_\gamma\right) - \frac{11}{6}(\kappa^{-1} + 0.84) \right].$$

[1] G. Segré and A. Silberberg, Behaviour of macroscopic rigid spheres in poiseuille flow part 2. experimental results and interpretation, J. Fluid Mech. 14, 136 (1962).

[2] B. P. Ho and L. G. Leal, Inertial migration of rigid spheres in two-dimensional unidirectional flows, J. Fluid Mech. 65, 365 (1974).

[3] P. Vasseur and R. G. Cox, The lateral migration of a spherical particle in two-dimensional shear flows, J. Fluid Mech. 78, 385 (1976).

[4] J. A. Schonberg and E. J. Hinch, Inertial migration of a sphere in Poiseuille flow, J. Fluid Mech. 203, 517 (1989).

[5] E. S. Asmolov, The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number, J. Fluid Mech. 381, 63 (1999).

[6] J. Oakey, R. W. Applegate Jr., E. Arellano, D. Di Carlo, S. W. Graves, and M. Toner, Particle focusing in staged inertial microfluidic devices for flow cytometry, Anal. Chem. 82, 3862 (2010).

[7] J. F. Edd, D. Di Carlo, K. J. Humphry, S. Köster, D. Irimia, D. A. Weitz, and M. Toner, Controlled encapsulation of single cells into monodisperse picolitre drops, Lab Chip 8, 1262 (2008).

[8] S. C. Hur, H. Tat K. Tse, and D. Di Carlo, Sheathless inertial cell ordering for extreme throughput flow cytometry, Lab Chip 10, 274 (2010).

[9] J.-P. Matas, J. F. Morris, and É. Guazzelli, Inertial migration of rigid spherical particles in poiseuille flow, J. Fluid Mech. 515, 171 (2004).

[10] D. Di Carlo, D. Irimia, R. G. Tompkins, and M. Toner, Continuous inertial focusing, ordering, and separation of particles in microchannels, Proc. Natl. Acad. Sci. USA 104, 18892 (2007).
[11] A. A. S. Bhagat, S. S. Kuntaegowdanahalli, and I. Papautsky, Enhanced particle filtration in straight microchannels using shear-modulated inertial migration, *Phys. Fluids* **20**, 101702 (2008).

[12] Y. S. Choi, K. W. Seo, and S. J. Lee, Lateral and cross-lateral focusing of spherical particles in a square microchannel, *Lab Chip* **11**, 460 (2011).

[13] J. P. Matas, V. Glezer, É. Guazzelli, and J. F. Morris, Trains of particles in finite-Reynolds-number pipe flow, *Phys. Fluids* **16**, 4192 (2004).

[14] B. Chun and A. J. C. Ladd, Inertial migration of neutrally buoyant particles in a square duct: An investigation of multiple equilibrium positions, *Phys. Fluids* **18**, 031704 (2006).

[15] V. Loisel, M. Abbas, O. Masbernat, and E. Climent, Inertia-driven particle migration and mixing in a wall-bounded laminar suspension flow, *Phys. Fluids* **27**, 123304 (2015).

[16] S. Kahkeshani, H. Haddadi, and D. Di Carlo, Preferred interparticle spacings in trains of particles in inertial microchannel flows, *J. Fluid Mech* **786**, R3 (2016).

[17] Y. Gao, P. Magaud, L. Baldas, C. Lafforgue, M. Abbas, and S. Colin, Self-ordered particle trains in inertial microchannel flows, *Microfluid. Nanofluid*. **21**, 154 (2017).

[18] G. G. Poe and A. Acrivos, Closed-streamline flows past rotating single cylinders and spheres: Inertia effects, *J. Fluid Mech*. **72**, 605 (1975).

[19] Y. Yan, J. F. Morris, and J. Koplik, Hydrodynamic interaction of two particles in confined linear shear flow at finite Reynolds number, *Phys. Fluids* **19**, 113305 (2007).

[20] W. Lee, H. Amini, H. A. Stone, and D. Di Carlo, Dynamic self-assembly and control of microfluidic particle crystals, *Proc. Natl. Acad. Sci. USA* **107**, 22413 (2010).

[21] G. Wang, M. Abbas, and E. Climent, Modulation of large-scale structures by neutrally buoyant and inertial finite-size particles in turbulent Couette flow, *Phys. Rev. Fluids* **2**, 084302 (2017).

[22] I. Calmet and J. Magnaudet, Large-eddy simulation of high-schmidt number mass transfer in a turbulent channel flow, *Phys. Fluids* **9**, 438 (1997).

[23] M. R. Maxey and B. K. Patel, Localized force representations for particles sedimenting in stokes flow, *Int. J. Multiphase Flow* **27**, 1603 (2001).

[24] S. Lomholt and M. R. Maxey, Force-coupling method for particulate two-phase flow: Stokes flow, *J. Comput. Phys.* **184**, 381 (2003).

[25] K. M. Yeo and M. R. Maxey, Dynamics and rheology of concentrated, finite-Reynolds-number suspensions in a homogeneous shear flow, *Phys. Fluids* **25**, 053303 (2013).

[26] K. Yeo, S. Dong, E. Climent, and M. R. Maxey, Modulation of homogeneous turbulence seeded with finite-size bubbles or particles, *Int. J. Multiphase Flow* **36**, 221 (2010).

[27] D. R. Mikulencak and J. F. Morris, Stationary shear flow around fixed and free bodies at finite Reynolds number, *J. Fluid Mech.* **520**, 215 (2004).

[28] H. Haddadi and J. F. Morris, Topology of pair-sphere trajectories in finite inertia suspension shear flow and its effects on microstructure and rheology, *Phys. Fluids* **2**, 1 (2015).

[29] A. J. Goldman, R. G. Cox, and H. Brenner, Slow viscous motion of a sphere parallel to a plane wall-II Couette flow, *Chem. Eng. Sci.* **22**, 653 (1967).

[30] P. Cherukat and J. B. Mclaughlin, The inertial lift on a rigid sphere in a linear shear flow field near a flat wall, *J. Fluid Mech.* **263**, 1 (1994).

[31] E. S. Asmolov, A. L. Dubov, T. V. Nizkaya, J. Harting, and O. I. Vinogradova, Inertial focusing of finite-size particles in microchannels, *J. Fluid Mech.* **840**, 613 (2018).

[32] K. Hood, S. Lee, and M. Roper, Inertial migration of a rigid sphere in three-dimensional poiseuille flow, *J. Fluid Mech.* **765**, 452 (2015).

[33] M. Abbas, P. Magaud, Y. Gao, and S. Geoffroy, Migration of finite-sized particles in a laminar square-channel flow from low to high Reynolds numbers, *Phys. Fluids* **26**, 123301 (2014).

[34] G. Subramanian and D. Koch, Inertial effects on the transfer of heat or mass from neutrally buoyant spheres in a steady linear velocity field, *Phys. Fluids* **18**, 1 (2006).

[35] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevFluids.3.114302 for sequence of particle positions in the channel.

[36] D. Klotsa, M. R. Swift, R. M. Bowley, and P. J. King, Interaction of spheres in oscillatory fluid flows, *Phys. Rev. E* **76**, 056314 (2007).
[37] K. Hood and M. Roper, Pairwise interactions in inertially driven one-dimensional microfluidic crystals, *Phys. Rev. Fluids* **3**, 094201 (2018).

[38] K. J. Humphry, P. M. Kulkarni, D. A. Weitz, J. F. Morris, and H. A. Stone, Axial and lateral particle ordering in finite-Reynolds-number channel flows, *Phys. Fluids* **22**, 081703 (2010).

[39] P. G. T. Saffman, The lift on a small sphere in a slow shear flow, *J. Fluid Mech.* **22**, 385 (1965).

[40] J. Magnaudet, Small inertial effects on a spherical bubble, drop or particle moving near a wall in a time-dependent linear flow, *J. Fluid Mech.* **485**, 115 (2003).