Inertial focusing of neutrally buoyant particles in heterogenous suspensions

Abstract The modelling-based design of microfluidic devices leads to highly efficient process intensification, which provides insights into different temporal and spatial scales at which processes in various fields of application could be performed. This requires not only an understanding of the underlying mechanisms of different processes at the micro scale, but also the development of relevant computational tools. The macroscopic models are often unable to produce conclusive evidence for a given mechanism in systems with the complexity characterizing almost all chemical and biochemical processes. By contrast, mesoscale methods possess the unique ability to model relatively large physical systems, and, at the same time, effectively capture the essential features of the micro- and nanoscale structure, architecture, and relevant interactions. We demonstrate the feasibility and usefulness of this novel tool by considering a movement of neutrally buoyant particles in straight microchannels. The two-dimensional lattice Boltzmann method with immersed boundary conditions was used to study the influence of Reynolds number and particle diameter ratio on formation of particle trains. It was shown that an increase in particle diameter ratio leads to a less stable final particle configuration. An increase in Reynolds number was not found to significantly influence the train stability in the tested range.


Introduction
The development and widespread use of micro-and nanofluidic technologies has significantly influenced process design ranging from laboratory to industrial scales. Important applications include (but are not limited to) single cell analysis [1], biocatalysts [2], fuel cell technologies [3] and biosensors [4]. The aforementioned micro-and nanofluidic approaches offer process intensification with more efficient heat and mass transport and in many cases also environmentally friendly technologies. In this context reliable numerical simulations present an invaluable tool which on one hand helps to unravel the principles behind the complex fundamental phenomena and process optimization on the other. Common simulation approaches include standard CFD techniques [5], lattice Boltzmann method (LBM) [6][7][8][9][10][11], dissipative particle dynamics (DPD) [12] and smoothed particle hydrodynamics [13].
A promising technology which today is mainly restricted to the micro scale is inertial microfluidics. It is based on the inertial lift effect first reported in Nature by Segré and Silberberg [14] in 1961. Their key observation was that randomly dispersed particles in a macroscopic pipe flow migrate to their specific equilibrium position of 0.6R (R is the pipe radius) when the Reynolds number is sufficiently low and the particle size is not negligible compared to the pipe diameter. A similar behaviour was observed for rectangular channels where particles additionally focus across the channel cross section [15]. In the latter case the number of equilibrium positions together with focusing pattern can be tuned by varying the Reynolds number and channel aspect ratio.
Since its first discovery in the 1960s, several experimental and theoretical/computational studies have been devoted to understanding of particle focusing in channels of various geometries and flow conditions e. g. [16][17][18][19]. Besides lateral and cross-sectional focusing particles can also exhibit longitudinal ordering along the flow direction which results in the formation of staggered particle trains [20]. Until now the ability to form particle trains was mainly investigated in monodisperse suspensions. Apart from an experimental study performed by Gao et al. [21] the effect of bi-and polydispersity has not been thoroughly addressed.
To better understand the focusing ability of rigid, neutrally buoyant particles in bidisperse suspensions and to fill the gap we performed extensive 2D numerical simulations based on the lattice Boltzmann method coupled to immersed boundary conditions (IB-LBM). The influence of Reynolds number and particle diameter ratio on formation of staggered trains was investigated. In the present study we demonstrate that an increase in particle diameter ratio leads to a less stable final particle configuration. An increase in Reynolds number was not found to significantly influence the train stability in the tested range but did shorten the focusing times.

Methods
Lattice Boltzmann method (LBM) is based on the discretized version of the Boltzmann equation [6,22]. A central quantity of interest is the density distribution function f(x,c, t) where x and c are position and velocity at time t respectively. In a discretized form an appropriate velocity set is needed. We employ a standard D2Q9 velocity set where velocity vectors are conveniently represented as: Using this notation, we can associate a density distribution function f i (x, t) to particle population at the position x and time t moving in the direction c i . From the density distribution function macroscopic density and velocity fields are recovered using: In its simplest form the Boltzmann equation can be discretized using the Bhatnagar-Gross-Krook collision operator [23], which yields: Left-hand side of Eq. 4 is known as streaming, while the right-hand side is commonly referred to as collision. Any external forces are introduced via F i ext . Here τ represents the relaxation parameter which is related to the fluid viscosity by: In all our simulations the relaxation parameter is set to τ = 1.0 and the speed of sound c s to c s ¼ 1= ffiffiffi 3 p . The equilibrium distribution function is given by: For particulate flows it is especially powerful to combine the lattice Boltzmann method to immersed boundary conditions (IB-LBM) [24]. Within the IB-LBM framework two coordinate systems are defined. A Lagrangian system is where particles are represented as a set of points referred to as Lagrangian markers while the Eulerian grid is where the standard LBM process takes place and therefore represents the surrounding fluid. To establish a communication between the two coordinate systems we employed an algorithm named multi direct forcing (MDF) [25]. A central idea of MDF is to iteratively satisfy the no-slip condition at the particle-fluid boundary through steps called velocity interpolation and force spreading. For a detailed algorithm and comparison of various schemes we refer the reader to the excellent review given by Kang and Hassan [26]. Here we only highlight points relevant to our model. For velocity interpolation and spreading of forces back to Eulerian grid we used the simplest form of the discretized delta function Δ(x) = ϕ(x)ϕ(y) as first proposed by Peskin in his doctoral thesis [27]: and the Lagrangian correction force taking form of: Correction forces were incorporated through an external force term as: This type of setup has been successfully used before to simulate inertial focusing of single neutrally buoyant particle in stratified flows [28].
One of the drawbacks associated with the IB-LBM is the fact that the liquid is contained in the particle itself. It is therefore necessary to use a suitable time integrator which prevents the internal fluid from affecting the motion of a particle significantly. According to Feng in Michaelides [29] particle translational and angular velocities should be updated using: Here m p and I p represent the particle mass and mass moment of inertia respectively. Note that for our case of a neutrally buoyant particle term containing gravitational acceleration g can be omitted. Particle positions are updated using the first order forward Euler scheme: and the velocities at the Lagrangian markers are adjusted at each timestep according to: where the subscript k denotes the k-th Lagrangian marker. To prevent particles from interpenetrating each other we used a short distance repulsive model as proposed by Glowinski et al. [30]: where M is the particle mass, U is the particle translational velocity, R its radius and ε = 10 −4 . We set d min = 2R and Δr = 2Δx, respectively.
Our simulation setup is summarized as follows. At the beginning of each simulation we randomly distribute 6 smaller and 6 larger cylindrical particles in a channel with dimensions 2000 × 300 Eulerian points. To save computation time, particles were not allowed to be placed at the channel centreline initially. The distance between neighbouring Lagrangian markers is set to d = 0.48Δx and kept constant throughout the study. Initially particles and fluid are at rest. The flow is invoked and maintained through application of a constant body force which is in- . At the channel inlet and outlet, we implemented periodic boundary conditions while at the top and bottom half-way bounce back is used. The use of periodic boundaries does not exactly coincide with an experimental setup. However, periodic boundary conditions are routinely used in such simulations and do not influence the results provided that the channel is long enough (see for example [28]). We performed a total of 4 × 10 6 simulation steps for each set of conditions. Simulations were performed using in-house developed computer code written in Fortran 90.

Results and discussion
In the following subsections we provide results for particle distributions perpendicular to the main flow direction, some examples of simulated velocity fields, estimated focusing times and trace the train instabilities which occur at high particle diameter ratios. In our simulations Reynolds number ranges from Re = 50 to Re = 130. Note that we define the Reynolds number for a channel without particles as: The combinations of radii used in our simulations along with their corresponding volume fractions are given in Table 1. All our results are given in lattice units. No conversion to physical units has been performed during the study.

Train stability depends on the particle diameter ratio
When the steady state is reached, particles generally focus on specific lateral positions symmetrically with respect to the channel centreline due to inertial forces. As a first step to investigating the train stability it is therefore instructive to compute particle distributions perpendicular to the main flow direction. Here the distribution width serves as a loose criterion for the train stability and helps to identify cases which deviate from the generally accepted trends observed for the focusing of a single particle or train formation in homogenous suspensions. For our analysis we sampled particle configurations from the second half of each simulation.
As can be deduced from the distributions presented in Fig. 2, larger particles are always shifted closer to the channel centreline. This type of behaviour has been reported for single particle dynamics as well as for ordering of homogenous suspensions [28,31]. We therefore   conclude that heterogeneity does not influence the steady-state particle lateral position. At this point we only make qualitative comparisons with experimental data, since a detailed comparison of 2D simulations and 3D physical world is not fully justifiable.
Distribution widths also reveal that as the particle diameter ratio R 2 / R 1 is close to 1.0, suspensions mimic behaviour of homogenous systems, i. e. particles are narrowly focused at their equilibrium positions. We display two examples of simulated velocity fields for such conditions on Fig. 3. As the diameter ratio increases beyond 2.0 along with volume fraction increase, however, the distributions of smaller particles start do display a significant broadening. An analysis of lateral positions over time, which we provide in the Supplementary information, reveals that a set of conditions where R 1 = 10, R 2 = 40 and Re = 50 should be excluded from this debate, since in that particular case, broadening is the result of slow focusing as a consequence of low Reynolds number and motion through crowded media. In most other cases where R 2 /R 1 is larger than 2.0 and the volume fraction is greater than 5%, smaller particles tend to oscillate in a regular fashion around their equilibrium positions (Fig. 4). We display two examples of simulated velocity fields for high diameter ratio cases on Fig. 5. These findings are concisely represented on Fig. 6.
We stress that this type of behaviour is not simply a consequence of increased volume fraction, since for example a combination R 1 = 30 and R 2 = 40 (7.9 vol%) does not display such oscillations, whereas a set R 1 = 10 and R 2 = 40 (5.3%) does. Simulated velocity fields for every set of conditions are given in the Supplementary material.
A change in the Reynolds number did not result in any appreciable difference in the tested range.
As we have briefly shown in the previous paragraphs, the reason for broadening of particle distribution functions lies in oscillations which develop when both the concentration and particle diameter ratio are high enough. To verify that the oscillations are not an artifact of boundary conditions, we repeated some simulations using periodic length of L = 3000 instead of L = 2000 while keeping volume fractions intact. Keeping in mind that each starting configuration is random, resulting oscillations were comparable. On Fig. 6, the distribution widths are expressed as standard deviations from mean particle positions. However, concentrations and particle diameter ratio being high enough are not sufficient for oscillations to develop. On Fig. 5 we have shown that such motion is a consequence of particles forming pairs or even triples which collectively move throughout the channel. Therefore   after an initial transient period particles have to "meet" in order to hydrodynamically interact and start their oscillatory behaviour which does not decay over time (Fig. 7). A similar behaviour has been investigated for a pair of particles under pressure driven flow [32]. In other words, if particles do not meet in order to interact, oscillations will not develop and formed train will be classified as stable even though concentration and diameter ratio might be high enough. Given that we started our simulations from random particle configuration, the results are unbiased with respect to pair formation. We also emphasize that in all cases where oscillations did develop, they were observed strictly for smaller particles. We have never observed pair formation between a larger and smaller particle or between two large particles. In the tested range, Reynolds number was not found to influence pair or triple formation.

Focusing times
In the present work we define the focusing time as an average time needed for a particle of certain radius to reach its equilibrium lateral position. In microfluidics it is a common practice to measure the entrance length rather than the focusing time. Due to periodic boundary conditions, however, entrance lenght is an ill-defined quantity. Therefore we provide focusing time which is measured in simulation steps. In order to estimate the focusing time, trajectories of particles of certain radius at a given Reynolds number were averaged in order to remove the influence of initial height on the focusing time. Oscillations which develop in some cases make the estimation of focusing times relatively problematic, resulting in larger error associated with its focusing time estimate.
The analysis of Fig. 8 reveals that an increase in particle radius leads to a shorter focusing time. This trend is retained over the range of Reynolds numbers at which simulations were performed. An increase in Reynolds number generally shortens the focusing time. Both effects are more pronounced when particles are smaller and vanish with an increase in radius. Both trends correspond well with dynamics observed in single particle studies or studies which were concerned with homogenous suspensions. We therefore conclude that heterogeneity does not influence individiual particle focusing times. As long as no hydrodynamic interactions which would result in pair or triple formation take place, particles will behave as individuals.

Conclusions
In our 2D computational study we have thoroughly investigated hydrodynamic focusing of neutrally buoyant particles in heterogenous suspensions. In our simulation the channel Reynolds number spans from Re = 50 up to Re = 130, which is a range which may find relevance in microfluidic applications. To the best of our knowledge there exists only one experimental study where authors addresed effects of heterogeneity on the aforementioned phenomena. However, due to the fact that the physical world is 3D and our simulations are 2D we did not make any direct comparisons with their results since it is well known that soft matter often behaves differently when confined to 2D. Additionally, volume fractions in their study are for at least an order smaller than those we have used.
Our main conclusions are summarized as follows. In heterogenous suspensions, train stability is critically influenced by volume fraction of particles and particle diameter ratio. When the diameter ratio is close to 1.0, suspensions will mimic the behaviour of homogenous systems. As the diameter ratio is increased beyond approximately 2.0, smaller particles often tend to form pairs or even triples which oscillate and degrade the train stability. If, on the other hand, particles never come close enough in order to engage in hydrodynamic interaction, the train will be stable even above before established treshold. In the tested range we have not found any correlation between the train stability and the Reynolds number. An increase in Reynolds number and particle size did shorten the focusing times, which agrees with results reported for single particle dynamics and homogenous suspensions. In the future, studies are needed which will address higher Reynolds numbers and higher degrees of heterogeneity. Further and more detailed studies are also required to elucidate mechanisms of particle oscillations in such heterogenous suspensions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.   7. Evolution of pairwise distances in steady state for a pair (left) and a triple (right) of interacting particles. As the number of interacting particle grows, the complexity of oscillations grows as well. Clearly, oscillations do not decay over time. An interacting pair was observed in a suspension where R 1 = 20, R 2 = 50 and a triple was observed in a suspension where R 1 = 10, R 2 = 40. In both cases the Reynolds number was set to Re = 100.