A pair of particles in inertial microfluidics: effect of shape, softness, and position

Abstract

Lab-on-a-chip devices based on inertial microfluidics have emerged as a promising technique to manipulate particles in a precise way. Inertial microfluidics exploits internal hydrodynamic forces and the mechanical structure of particles to achieve separation and focusing. The article focuses on the hydrodynamic interaction of two particles. This will help to develop an understanding of the dynamics of particle trains in inertial microfluidics, which are typical structures in multi-particle systems. We perform three-dimensional lattice Boltzmann simulations combined with the immersed boundary method to unravel the dynamics of various mono- and bi-dispersed pairs in inertial microfluidics. We study the influence of different starting positions for mono- and bi-dispersed pairs. We also change their deformability from relatively soft to rigid and choose spherical and biconcave particle shapes. The observed two-particle motions in the present work can be categorized into four types: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair. We show that stable pairs become unstable when increasing the particle stiffness. Furthermore, a pair with both capsules in the same channel half is more prone to become unstable than a pair with capsules in opposite channel halves.

---

1 Introduction

Biomedical studies of biological cells play a crucial role in diagnosing several fatal diseases¹ such as malaria,² cancer,^3,4 and the human immunodeficiency virus (HIV) infection,⁵ to name a few. In recent years, microfluidic⁶ lab-on-a-chip devices⁷ have emerged as a promising technique to precisely manipulate and control biological cells needed on a commercial level.⁸ Lab-on-a-chip devices for separating particles and biological cells rely on external force fields or internal hydrodynamic forces.⁹ Lab-on-a-chip techniques using external force fields include optical tweezers,¹⁰ dielectrophoresis,¹¹ magnetophoresis,¹² and acoustophoresis.¹³ On the other hand, deterministic lateral displacements at low Reynolds number using appropriately placed pillars¹⁴ and inertial microfluidics^15–20 exploit internal hydrodynamic forces to achieve particle separation. In contrast to common microfluidic lab-on-a-chip devices in which fluid inertia is negligible, inertial microfluidics operates in an intermediate range between Stokes and turbulent regimes, where flow is still laminar. The possibility to attain high throughput makes inertial microfluidics an attractive option among other microfluidic particle-separation techniques. In this article, we investigate how different factors such as shape, softness, and position influence the motion of a pair of particles in a microfluidic channel.

Inertial focussing was first reported by Segré and Silberberg for solid particles with the well-known tubular pinch effect^21,22 and then spurred numerous theoretical,^23–26 computational,^19,27–30 and experimental^31–34 studies to gain a deeper understanding of this phenomenon. Inertial migration of a neutrally buoyant solid particle is essentially the result of a balance between shear-gradient lift force, directed towards the channel wall,³⁵ and a wall repulsion force, directed towards the channel center.³⁶ Applying an external force to the solid particle along the channel axis, the resulting Saffman force³⁷ also contributes significantly to the force balance in the cross-stream direction. For almost five decades from the observation of the Segré–Silberberg effect, inertial migration was mainly a blue-skies-research problem. In 2007 Di Carlo et al.³⁸ demonstrated ordering and separation of particles in a microchannel using inertial focussing, which gave rise to the field of inertial microfluidics. Subsequently, Hur and co-workers³⁹ experimentally showed inertia-driven separation and enrichment of deformable capsules. Soft capsules experience an additional lift force induced by the deformability of the cell, which drives them towards the channel center and which thereby competes with the inertial lift force.⁴⁰

For biomedical applications of inertial microfluidics it is crucial to strive for a deeper understanding of the inertial migration of deformable capsules. The dynamics of soft capsules is widely studied at low Reynolds numbers in canonical flows^41–45 and cellular blood flow simulations.^46–49 Only recently, inertial microfluidics of deformable cells has caught attention.^20,50–58 Using two-dimensional numerical simulations, Shin and Sung⁵¹ showed that the final equilibrium position for a given value of capsule elasticity varies with the Reynolds number. In their following work⁵⁴ they identified different types of capsule dynamics such as tank-treading, swinging, and tumbling. According to Kilimnik et al.⁵⁰ the steady-state position of a deformable capsule in the lateral direction is independent of the flow rate and it only depends on the softness. This was later verified by Schaaf and Stark²⁰ using 3D lattice Boltzmann simulations. They also proposed that the inertial lift force in the vicinity of the channel center scales as a cube of capsule radius. More recently, Alghalibi and colleagues⁵⁸ investigated cross-stream inertial migration of a soft particle in pipe flow. Interestingly, they found that irrespective of the particle's softness and the flow rate, the particle always ends up at the channel center for an initial particle diameter of 0.2 times the pipe diameter.

Practical applications require manipulation of several soft particles at a time and, hence, it is important to study how particles interact in flow. There are a few systematic studies in literature to understand the interaction of solid and soft particles in inertial microfluidics. Yan et al.⁵⁹ investigated the dynamics of two solid particles in linear shear flow generated by two plates moving against each other. They observed that particles either settle down on the centerline where flow velocity is zero or exhibit oscillatory motion when released from symmetric initial positions. Schaaf et al.⁶⁰ reported that for two closely placed solid particles in a channel flow, the inertial lift forces acting on the leading and lagging particles differ significantly. In their subsequent work⁶¹ they extended the understanding of a particle pair to study the behavior of particle trains⁶² in inertial microfluidics. At low Reynolds numbers, Lac et al.⁶³ and Aouane et al.⁶⁴ investigated the interaction of two soft capsules in shear and Poiseuille flows, respectively. The inertial counterparts of these two problems were also studied. Doddi and Bagchi⁶⁵ observed a transition from irregular self-diffusion⁶⁶ to different types of spiraling motions (outward, fixed orbit, and inward) with increasing Reynolds number. Lan and Khismatullin⁶⁷ reported that a pair of deformable cells in a rectangular microchannel undergoes either swapping or passing trajectories, while deformable cells in a series experience damped oscillations. Like Lan and Khismatullin,⁶⁷ Schaaf et al.⁶⁰ also observed these swapping and passing trajectories for a pair of solid particles. Krüger et al.⁶⁸ looked into the suspension of deformable capsules in Poiseuille flow at finite Reynolds numbers. They discovered that particle–particle interactions enhance the migration towards the channel center for sufficiently large Reynolds and capillary numbers.

In the present work we perform three-dimensional lattice-Boltzmann simulations combined with the immersed boundary method^20,60,69 to investigate the dynamics of mono- and bi-dispersed particle pairs. They flow through a square microchannel in the inertial regime. We systematically investigate the effect of different starting positions for mono- and bi-dispersed pairs. We also vary their deformability from relatively soft to rigid and choose spherical and biconcave particle shapes. Reynolds number, particle radius, and axial particle distance are kept the same in all the cases. The observed dynamics of the different mono- and bi-dispersed particle pairs can be categorized into four types: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair.

The outline of the article is as follows. In Section 2 we discuss the computational setup and simulation parameters of the problem and explain the lattice Boltzmann simulations. Results from the numerical investigations are presented in Section 3 and we close with concluding remarks in Section 4.

2 Computational setup and numerical methodology

2.1 Microfluidic setup for a flowing pair of particles

Fig. 1 shows the schematic of our microfluidic setup. We study the hydrodynamic interaction and inertial migration of a flowing pair of soft particles in a quadratic microchannel. We start from a fully-developed Poiseuille flow and place two particles in the midplane at y = 0. We apply periodic boundary conditions in the z-direction along the channel axis and the no-slip boundary condition at the channel walls in the remaining directions. The length of the microchannel is large enough so that particles do not interact with their periodic images unless the axial separation becomes as large as the channel length. The pressure gradient or force density α required to pump the fluid through the channel is calculated from the analytical expression for the velocity field of a Poiseuille flow:⁷⁰where U_z, W, and η are the axial velocity, channel half-width and dynamic viscosity of the fluid, respectively.

Fig. 1 Computational setup for a pair of neutrally buoyant flowing particles in a microchannel. Two particles of non-dimensional radius a and initial axial separation Δz are placed in a fully developed Poiseuille flow in the midplane at y = 0.

2.2 Coupled lattice Boltzmann finite element immersed boundary method (LBFEIBM)

2.2.1 Lattice Boltzmann method. In the last two decades the lattice Boltzmann method^71,72 has emerged as an efficient alternative to conventional techniques of computational fluid dynamics. In particular, it was applied to multiphysics problems involving soft matter and fluid flow.⁷³ The lattice Boltzmann method operates at the mesoscale and provides a solution of the discretized Boltzmann equation. It is based on the lattice Boltzmann equation that governs the advection and collision of populations of particles. They are described by distribution functions f_i(x,t) for a discrete set of velocity vectors c_i. In time step Δt these distributions evolve according towherestands for the local equilibrium distribution of velocity vectors c_i relative to flow velocity u. Here, we implement the D3Q19 lattice with 19 discrete velocity vectors chosen from c_i = {(0, 0, 0), cyc(±1, 0, 0), cyc(±1, ±1, 0)} with weights w_i = {1/3, 1/18, 1/36}. Furthermore, c_s in eqn (3) is the speed of sound and the relaxation time τ in eqn (2) is related to the kinematic viscosity , where in lattice units. The last term in eqn (2) is due to the forcing scheme proposed by Guo et al.⁷⁴ It takes into account the force density driving the fluid and forces exerted by particles. Based on this scheme, a body force density f is related to F_i as

In order to obtain the macroscopic picture of the flow field, fluid density ρ and momentum density ρu at each lattice node can be calculated using the zeroth and first moments of the populations f_i:

We enforce the no-slip condition at channel walls using the regularized method of Latt and Chopard.^75,76 The numerical accuracy of our simulations is maintained by setting the relaxation time τ ≤ 1. Please refer to ref. 20 for further details. The numerical solution of the lattice Boltzmann equation is obtained using an open-source parallel lattice Boltzmann solver from the Palabos project,⁷⁷ while implementation of the immersed boundary method and particle dynamics uses in-house libraries.

2.2.2 Modeling of soft capsules. We employ the finite element method to discretize the surface of the particle with elements,⁷⁸i.e., linear triangular elements. Tessellation of the particle surface is achieved by dividing each triangle of an icosahedron into four triangles (for more details refer to ref. 69) and the newly created nodes are then radially shifted to the circumsphere. The process is repeated until the distance between two neighboring nodes approaches the lattice spacing. This procedure of generating mesh from a highly symmetric solid (icosahedron in the present case) enables us to achieve superior mesh quality in terms of isotropy and homogeneity.⁶⁹ The same procedure is used for the solid particles as well.

The capsule dynamics is determined by a combination of strain, bending, and volume energies. For the strain energy we rely on the Skalak model⁷⁹ valid for an isotropic and homogeneous capsule. It is formulated as a function of the strain invariants of the Green strain tensor and given as

where k_s is the shear modulus and k_a is the area dilation modulus. The strain invariants I₁ and I₂ can be expressed in terms of in-plane principal extension ratios λ₁ and λ₂:

I₁ = λ₁² + λ₂² − 2, (8)

I₂ = λ₁²λ₂² − 1. (9)

The discretized version of the bending energy proposed by Helfrich^80,81 takes the form

where k_b is the bending modulus, θ_ij is the angle between the area normal vectors of two neighboring surface elements, and θ⁰_ij is the reference angle from the undeformed (initial) state. Implementation of the bending energy prevents cusp formation and buckling of soft capsules.

In this work, we consider soft particles with an impermeable membrane. To implement the constraint of constant volume, at least approximately, we choose a standard method and incorporate the volume energy⁸²

where k_v, V and V₀ are the volume modulus, actual capsule volume, and reference volume, respectively. To keep the volume changes below 1%, we choose k_v = 75 000ρν²/W². Note, these changes do not affect the dynamics of the particle.

Following Krüger et al.⁶⁸ and Schaaf and Stark,²⁰ we set k_a/k_s = 2 and k_b/(k_sa²) = 2.87 × 10⁻³. We can compare our parameters with the ones for red blood cells. To mimic their the mechanical behavior in blood flow simulations, the ratio k_b/(k_sa²) is set to 2.36 × 10⁻³,⁴⁷ which is very close to the value in present simulations. However, the ratio k_a/k_s in current simulations is significantly smaller than that used for a healthy red blood cell, which allows our spherical capsules to deform significantly.

The total force acting on node i can be obtained from the principle of virtual work;⁸³ given as F_i = ∂[E_s + E_b + E_v]/∂x_i. The derivatives of strain, bending, and volume energies were derived analytically and then inserted directly into the solver to save computational costs. The analytical form of the derivatives can be found in ref. 84.

2.2.3 Modeling of solid particles. The motion of a neutrally buoyant solid particle is governed by Newton's law and its equivalent for the angular momentum:

S(t + Δt) = S(t) + U(t), (12)

MU(t + Δt) = MU(t) + F_fluid + F_Feng, (13)

Iω(t + Δt) = Iω(t) + T_fluid + T_Feng, (14)

where, S, M, I, U, and ω are the particle position, mass, moment of inertia, linear velocity, and angular velocity, respectively. Force F_fluid and torque T_fluid exerted by the fluid on the particle are computed using the iterative procedure of Inamuro,⁸⁵ which ensures the no-slip condition at the boundary of the solid particle. Suzuki and Inamuro⁸⁶ demonstrated that for simulating the dynamics of solid particles at finite Reynolds numbers correctly, the effect of a fictitious fluid inside the solid particle must be taken into account when using the immersed boundary method. Therefore, one needs to introduce Feng's rigid body approximation⁸⁷ with the help of additional forces and torques, F_Feng and T_Feng, respectively.

2.2.4 Immersed boundary method. The hydrodynamics of soft and solid particles in inertial microfluidics is essentially a fluid–structure interaction (FSI) problem on the micron scale. In the present case, the motion of the bulk fluid in a microchannel is responsible for the movement of particles; in turn, particles exert forces on the fluid via their surfaces. We implement this interaction between fluid and particles using the immersed boundary method (IBM) of Peskin.^88,89 The immersed boundary method is one of the fluid–structure interaction methods that belong to the mixed Eulerian-Lagrangian framework. There also exist fully Eulerian and Lagrangian FSI methods (see Fig. 1 of ref. 90 for more details). The classical IBM treats vertices of the [Doublestruck P]₁ elements of the discretized particle surfaces as Lagrangian markers. The advection velocity ẋ_i of Lagrangian markers is obtained using an interpolation of the velocity values u from neighboring Eulerian points, which in our case are the lattice nodes of the lattice Boltzmann method:In contrast, the body force density f̃ exerted by Lagrangian markers on neighboring Eulerian points is computed using the spreading operation:Following ref. 88 the discretized delta function δ(X − x_i(t)) in eqn (15) and (16) are represented as δ(X − x_i(t)) = ϕ(X − x_i(t))ϕ(Y − y_i(t))ϕ(Z − z_i(t)), where

The interaction between soft capsules and the bulk fluid is resolved using the classical IBM. However, for solid particles we employ interpolation and spreading operations in an iterative manner^85,86 to obtain the force F_fluid and torque T_fluid directly from the flow field. Then, advection velocity U is obtained using eqn (13). This particular approach is known as the multi-direct-forcing variation of the immersed boundary method.^91,92

2.3 Simulation parameters

In our microfluidic setup we simulate two types of particle pairs, mono- and bi-dispersed (see Fig. 2). Mono-dispersed pairs consist of spherical particles of the same stiffness. Bi-dispersed pairs are classified into two subtypes; first, a pair with particles of different shapes but the same stiffness, and second, a pair with particles of varying stiffness but the same shape. We consider the combination of a spherical (leading) and biconcave (lagging) particle for the first type, and for the second type, two spherical particles. The initial orientation of a biconcave particle in the flow field is shown in a magnified view in Fig. 2. We apply the transformation
to the vertices of a spherical particle of radius a to model the biconcave shape.^93,94

Fig. 2 Classification of different types of particle pairs with radius a. Each pair type is simulated for the same initial distance Δz, for different values of softness La_lead, La_lag, and five sets of different initial positions x_lead, x_lag. Alphanumeric labels M1–M3 and B1–B4 are the identification tags given to each case. Note that the mono-dispersed pair with La = ∞ corresponds to the case of solid particles.

In all the cases we set the non-dimensional particle radius a and the initial axial separation Δz equal to 0.4 and 1, respectively. For the biconcave particle the radius characterizes the spherical particle from which it is generated using eqn (18). The value of the Reynolds number Re = 2WU_max/ν is fixed to 10 in all the simulations. Here, U_max is the maximum fluid velocity in the channel center and ν is the kinematic viscosity. Note, in a typical inertial microfluidics setup, the Reynolds number varies from 1 to 100¹⁷. Following Schaaf and Stark,²⁰ we use the Laplace number to define the softness of the particle; it is given as

La = k_sa/ρν². (19)

In our simulations we consider cases with relatively soft (La = 10) and stiff (La = 100) particles (please refer to Fig. 2). We also include a mono-dispersed pair of solid particles (case with La = ∞ in Fig. 2) to compare the dynamics of soft and solid particles of spherical shape. For comparison we note that for a red blood cell suspended in plasma, one has a = 4 × 10⁻⁶ m and k_s = 5.3 × 10⁻⁶ N m⁻¹,⁴⁷ and the Laplace number turns out to be 0.0212.

In practice, there might exist a viscosity contrast between the fluid within soft particles (cytoplasm in biological cells) and the bulk liquid (η_cyto/η_bulk > 1). For example, the cytoplasmic viscosity of the red blood cell is five times the plasma viscosity η_plasma= 1 mPa s.⁴⁷ Moreover, tumour cells such as A549 can have cytoplasmic viscosity in the range of 144.8 to 1390.7 Pa s,⁹⁵ which is several orders of magnitude higher than the plasma viscosity. We shortly summarize what influence this will have. An increase in the viscosity of the interior fluid or cytoplasm attenuates the deformation of soft particles and biological cells.⁵⁰ Thus, while for a very soft particle a viscosity contrast close to one should not significantly alter the migration dynamics, a viscosity contrast significantly larger than one makes particles less deformable and modifies their equilibrium positions. In the case of relatively stiff particles, the viscosity contrast will have a negligible influence on the particle dynamics. Moreover, for biconcave soft particles in inertial flows with Re ∼ [scr O, script letter O](1), there exists a critical value of viscosity contrast above which the transition from tank-treading to tumbling motion takes place.⁵² Later, we will see that biconcave particles in the present work already exhibit tumbling motion. With this in mind, we set the density and viscosity of the soft particle identical to the values of the bulk fluid.

In each simulation we fix the leading particle position in the x-direction to 0.2 and vary the position of the lagging particle from −0.4 to 0.4 in steps of 0.2. The leading particle's starting position is kept sufficiently far from the stable equilibrium point of isolated particles (see Fig. 5 and 9) so that disturbances from the lagging particle can have a significant effect on the motion of the leading particle. The implemented LBFEIBM solver has been employed in a series of works by our group^19,20,60,61 (see Appendix A for benchmarking). According to our previous work we discretize the channel width with 90 lattice nodes for the case of soft particles and 75 lattice nodes for solid particles. In our simulations, the necessity of a high grid resolution arises to ensure the proper coupling between particles and the bulk fluid. Moreover, relatively more lattice nodes are required in the case of soft particles to capture the particle deformation and lift force accurately. Lastly, the surface of soft and solid particles is discretized using 20 480 elements.

3 Results and discussion

In this section, we first explain how the overall dynamics of a particle pair is categorized in the present work for the different types of particle pairs as illustrated in Fig. 2. A detailed explanation of the hydrodynamic behavior is provided later in Sections 3.1 and 3.2 for the different mono- and bi-dispersed pairs, respectively.

In Poiseuille flow the inertial lift force migrates a particle towards a specific position between channel center and wall, where softer particles are located closer to the center. Depending on the direction of migration, the particle either accelerates or slows down. For two closely placed particles the lift force acting on each particle is modified compared to the single-particle case⁶⁰ and the modified cross-streamline motion causes particles to approach or move away from each other. Notably, the initial interaction and positions of particles determine their dynamics at a later stage.

In Fig. 3 we use the time evolution of the axial particle distance Δz to identify four different types in the dynamics of a flowing particle pair: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair. In our notation the stable particle pair either directly approaches a constant particle separation at long times (see video M1) or the particle distance oscillates. In the second case we find both, damped oscillations ( see video M2) and ongoing bounded oscillations (see video B1). The source of damped and bounded oscillations is different as we discuss below in Sections 3.1 and 3.2. Finally, particles of an unstable pair continuously drift apart until the particle distance approaches the channel length ( see video B3).

Fig. 3 Characterization of the dynamics of a particle pair based on the temporal evolution of the particle distance Δz. For all stable pairs the axial separation eventually equilibrates to a constant value or undergoes bounded oscillations. In contrast, the leading and lagging particles of an unstable pair drift apart until Δz approaches the channel length.

The diagram in Fig. 4 summarizes the dynamics of the particle pairs, which we observed in our simulations for different mono- and bi-dispersed pairs when the initial position of the lagging particle is varied. The mono-dispersed pair M1 with soft capsules is always stable irrespective of the starting position x_lag. It is situated in the channel center as we will see below. On increasing the particle stiffness from La = 10 to La = 100 (label M2), the pair is stable when the leading and lagging particles initially occupy different channel halves and becomes unstable when the particles start in the same channel half. The transition between both cases occurs through damped oscillations. For the mono-dispersed pair M3 composed of rigid particles (La = ∞), we observe that only the symmetric pair with (x_lead, x_lag) = (0.2, −0.2) is stable, reminiscent of the observations in ref. 60.

Fig. 4 Influence of the starting position x_lag on the overall dynamics of different mono- (M1–M3) and bi-dispersed (B1–B4) pairs with the leading particle starting at x_lead = 0.2 (see Fig. 2 for further details on labels M1–M3 and B1–B4).

The dynamics of the bi-dispersed pairs B1 and B2 is qualitatively the same as that of pairs M1 and M2, however, with the difference that now the B1 pair performs bounded oscillations. Lastly, the symmetric placement of leading and lagging particles in the B3 and B4 pairs gives a stable pair, which we also observe for the B3 pair with x_lag = 0. In summary, irrespective of particle-pair type when the leading and lagging particles start from the symmetric positions (x_lead, x_lag) = (0.2, −0.2) the pair is stable. Also, stable pairs become unstable when increasing the particle stiffness. Finally, a pair with both particles in the same channel half is more prone to become unstable than a pair with particles starting in the opposite channel halves.

Before proceeding, we add a comment. Prohm and Stark¹⁹ demonstrated that a freely placed solid particle in Poiseuille flow might either migrate to the main axis or the diagonal depending on the particle size and cross section of a microchannel. Furthermore, these zero lift-force positions can either be stable, saddle points, or unstable. Particles located on unstable positions migrate to stable positions in response to a small perturbation. In the work of Schaaf and Stark,²⁰ it was observed that the majority of implemented deformable capsules migrated to the channel diagonal. In our case we initialize all particles in the midplane at y = 0. We observe that in all cases they essentially stay in the midplane during migration since for symmetry reasons the lift force component normal to the midplane is zero.

3.1 Mono-dispersed pairs

3.1.1 Single particles. To understand similarities and differences in the dynamics of a mono-dispersed pair compared to a single particle, we briefly look at single particles with the same softness as in the pairs M1–M3. Fig. 5 shows the inertial migration of a spherical particle for different values of softness. Particle radius and Reynolds number are fixed to 0.4 and 10, respectively. To make position x and time t non-dimensional in Fig. 5, the channel half-width W and viscous time scale W²/ν are used, respectively. While rigid particles show the typical inertial focusing, the particle size in the present case is large enough so that the softness corresponding to La = 10 is sufficient for the lift forces due to deformability to dominate over the inertial lift forces. Hence, the particle settles at the channel center independent of the starting position. Note that smaller particles need smaller La or larger deformability to occupy the channel center.²⁰ We also show in the inset of Fig. 5 how the particle with rigidity La = 10 transforms to a bullet-shaped particle in steady-state, while the particle with rigidity La = 100 stays nearly spherical.

Fig. 5 Lateral position plotted versus time of a spherical particle moving in the midplane of the channel at y = 0 for different rigidities La and starting positions. Inset: The equilibrium shape of soft particles for different La is shown.

3.1.2 Particle pairs. For the mono-dispersed pairs M1–M3 the lateral positions of the leading and lagging particles and their axial distance are plotted versus time in Fig. 6. The softest particles (pair M1) migrate to the channel center (see video M1), however, the time needed by the leading and lagging particles to reach the center varies with the initial position of the lagging particle. For x_lag = 0 and −0.2 the leading and lagging particles in the M1 pair migrate with nearly constant axial distance in time. Furthermore, when the leading and lagging particles initially move on the same streamline, x_lead = x_lag = 0.2, the leading particle experiences a very strong push towards the channel center, which significantly enhances its axial migration velocity. However, once it has come close to the channel center in a very short time, it takes a longer time to fully reach the channel center. Finally, we note that, initially, the axial particle distance for the M1 pair with x_lag = 0.4 and −0.4 grows at the same rate. Nevertheless, the steady-state value of Δz for x_lag = 0.4 is twice as large as that for x_lag = −0.4. To quantify the hydrodynamic interaction of soft particles in the M1 pair, we compare the migration time t^m the leading and lagging particles need to reach the channel centerline with the corresponding time of a single particle. Such a comparison allows to demonstrate how strongly the total lift force is influenced by the neighboring particle in a pair. Thus, in Fig. 7, we plot the ratio t^m_pair/t^m_single for different values of the starting position x_lag. We find that the presence of the lagging particle with x_lag = ±0.4 slows down the migration of the leading particle to the channel center, while for the remaining three starting positions the migration of the leading particle is enhanced. Moreover, the lagging particle migrates faster when both particles are in different channel halves initially and vice versa.

Fig. 6 Lateral positions of the leading and lagging particles and their axial distances Δz plotted versus time for different initial positions x_lag and for fixed x_lead = 0.2. The mono-dispersed pairs M1–M3 with different rigidity La moving in the channel midplane at y = 0 are considered. Once the axial particle distance Δz is beyond 10 and close to the channel length, the simulations are stopped. Insets in the middle column show the behavior of x_lead, x_lag, and Δz at later times.

Fig. 7 Migration time t^m of particles in the mono-dispersed pair M1 relative to a single particle plotted versus the initial position x_lag for fixed x_lead. t^m is the time needed to reach the channel centerline with a distance smaller than 10⁻³. Note that it is not possible to include the point x_lag = 0 for the lagging particle since it is already at the equilibrium position for La = 10.

We now focus our attention on the dynamics of particles in the M2 pair. For x_lag = 0.2, similar to the M1 pair, the lagging particle initially pushes the leading particle towards the channel center. Thus, leading and lagging particles drift at a different axial velocity, which continuously increases the axial distance and makes the pair unstable. Note that we stop the simulations once the axial particle distance Δz is beyond 10 and therefore close to the channel length, where a particle interacts with the image of the other particle. Continuing the simulations of the leading particle from the momentary position at x = 0.055 without the other particle present shows that it migrates towards the center instead of drifting back to the known equilibrium position at x ≈ 0.3. Thus, we identify a second type of stable position at x = 0, which we attribute to the strong wall repulsion for the large particles, we use in our simulations, and the fact that the inertial lift force is small close to the center. The M2 pair with x_lag = 0.4 also turns out to be unstable since already initially, they drift with very different axial velocities. Interestingly, for x_lag = −0.4 the particle pair is stable although with a larger axial distance compared to the corresponding M1 pair with La = 10, while for the symmetric initial positions with x_lag = −0.2 the axial distance stays close to the starting value similar to the case La = 10.

A new feature is observed for the initial position x_lag = 0. The M2 pair is stable but reaches a constant axial distance close to the initial value only after extended damped oscillations in both the lateral positions and Δz (see video M2), as the insets in Fig. 6 show. Similar oscillations in the axial distance were reported earlier in numerical simulations of a pair of solid particles⁶⁰ and a train of soft particles⁶⁷ in a rectangular microchannel. Experimental evidence on damped oscillations of solid particles in a microchannel can be found in Lee et al.⁹⁶ Recently, based on theoretical calculations, Hood and Roper⁹⁷ suggested that in such oscillations, inertia and viscosity switch their roles, i.e., viscous interactions of a particle with the channel wall and other particles initiate oscillations while inertial forces dampen them.

We further investigated the nature of these damped oscillations for x_lag = −0.05, −0.1, and −0.15. We find that as the initial position of the lagging particle moves away from the channel center, the maximum amplitude of oscillations decreases [see Fig. 8(a)]. We assign this to the reduced viscous coupling, when the particles are further away from each other and an increased damping due to the (inertial) lift forces, which force the particle to a specific lateral position.⁶⁰ In all cases, particles equilibrate with an axial separation very close to its initial value. In the inset of Fig. 8(b) we visualize the trajectories of the leading and lagging particles in the center-of-mass frame for x_lag = 0. They reveal the inward spiraling motion of the leading and lagging particles. For decreasing x_lag the extent of the spirals shrink because of the reduction in the value of max(Δz). As the main plot in Fig. 8(b) shows, the damping rate of the oscillations in Δz first does not vary significantly for decreasing x_lag but then increases sharply.

Fig. 8 Damped oscillations in the M2 pair. Plotted as a function of x_lag, the graphs show in (a) the maximum amplitude max(Δz) and the steady-state value of the axial separation Δz_eq and in (b) the damping rate ζ of the oscillations. Inset in (b): trajectories of the leading and lagging particles in the center-of-mass frame for x_lag = 0.

Finally, we discuss the case of a pair of solid particles (M3) in Fig. 6. Only the symmetric placement with x_lag = −0.2 and x_lead = 0.2 provides a stable flowing pair, reminiscent of the findings in ref. 60. Unlike soft particles with x_lag = −0.2, the initial distance between two solid particles no longer remains constant over time. It grows until the particles have reached their final lateral positions. Earlier, we saw that in mono-dispersed pairs M1 and M2 with x_lag = ±0.4, the axial separation grows at the same rate only at early times. Now, the temporal increase of the axial distance Δz is the same for the whole time interval for these two cases since the lagging particles already start close to their equilibrium positions (see Fig. 5) and only the leading particle has to relax.

3.2 Bi-dispersed pairs

In this subsection we concentrate on the dynamics of different bi-dispersed pairs. We employ an approach similar to that in the previous discussion on mono-dispersed pairs. We first look at the migration of a biconcave particle so that later it can be compared with bi-dispersed pairs.

3.2.1 Single particles. Fig. 9 shows the temporal evolution of the lateral position of a biconcave particle with rigidities La = 10 and 100. It migrates to a stable equilibrium point very close to that shown in Fig. 5 for a spherical particle of the same rigidity. However, for La = 10 the inertial migration is significantly slower compared to the spherical particle of identical La value.

Fig. 9 Lateral position of a biconcave particle (left inset) plotted versus time when moving in the midplane of the channel at y = 0 for two values of rigidity and starting position. Right inset: Lift-off in the y-direction plotted versus time for La = 10.

At low Reynolds numbers a biconcave particle in shear flow undergoes either tank-treading, tumbling, or swinging motion depending on the shear rate, softness, and viscosity ratio.⁹⁸ Here, for more rigid particles at La = 100 we observe a solid-body like tumbling motion in the y = 0 plane, while for the softer biconcave particle at La = 10 a deformation-driven tumbling motion is observed as illustrated in Fig. 10. Initially [Fig. 10(b and c)], the particle rotates because the lower part being closer to the center moves faster. Later, when the particle aligns itself with the flow direction [Fig. 10(d)], it transforms into a parachute-like shape [Fig. 10(f)]. Eventually, the transition to a convex disc-like shape [Fig. 10(h)] completes one tumbling cycle, during which the particle has moved closer to the channel center. This also increases the period of the tumbling motion.

Fig. 10 Deformation-driven tumbling motion of a biconcave particle of rigidity La = 10 starting at x₀ = 0.2 in the midplane of the microchannel when looked from the top. Flow direction is from left to right and the cell tumbles counterclockwise. The snapshots are taken from the trajectory in Fig. 9.

Note that the biconcave particle with La = 10 does not fully reach the channel center at x = 0 and also leaves the midplane at y = 0 [inset of Fig. 9]. Therefore, it does not have a symmetric shape and performs a tumbling motion. We call it wobbling since the tumble axis varies in time, as we explain further below in Fig. 12. This is qualitatively similar to the observation of Fedosov et al.⁹⁹ for the tumbling-to-slipper transition of a red blood cell in cylindrical microchannels. A slipperlike form of the biconcave particle occurs in Poiseuille flow at low Reynolds numbers,¹⁰⁰ which in steady state undergoes tank-treading motion. Thus, we observe a transition from deformation-driven tumbling to wobbling when the particle approaches the center. It also occurs in the bi-dispersed pair B1 as discussed below.

3.2.2 Particle pairs. In Fig. 11 for the bi-dispersed pairs B1 and B2 we plot the lateral positions of the leading and lagging particles and their axial distance versus time. Similar to the M1 pair with La = 10 also the B1 pair with the same softness moves to the channel center. The axial distance in steady state is larger compared to the M1 pair. This is mainly because the biconcave particle migrates slower than spherical particles (compare Fig. 5 and 9). Due to the tumbling of the lagging biconcave particle, bounded oscillations in the axial distance and x_lead occur. Except for x_lag = − 0.2, all the other biconcave particles in the B1 pairs exhibit the transition from deformation-driven tumbling to wobbling close to the channel center as described earlier. We illustrate this transition for x_lag = 0.2 in Fig. 12 using the trajectory of a Lagrangian marker initialized on the periphery of the biconcave particle in the midplane at y = 0. Initially, the deformation-driven tumbling is indicated by the zig-zag trajectory, which then transitions to a helical trajectory due to wobbling (see also video B1). The different tumbling motions are highlighted using blue and red planes in Fig. 12. An exception is observed for the biconcave particle with initial position at x_lag = −0.2. It migrates very quickly towards the channel center in contrast to the single particle in Fig. 9 and performs directly the wobbling motion.

Fig. 11 Lateral positions of the leading and lagging particles and their axial distances Δz plotted versus time for different initial positions x_lag and for fixed x_lead = 0.2. The bi-dispersed pairs B1 and B2 with different rigidity as well as B3 and B4 consisting of leading and lagging particles with different rigidity La moving in the channel midplane at y = 0 are considered. Inset: The lift-off in the y-direction for the biconcave particle in the B1 pair is shown.

Fig. 12 A transition from the deformation-driven tumbling (blue plane) to wobbling (red plane) motion occurs when a biconcave particle approaches the channel center. The particle has rigidity La = 10 and moves in the B1 pair with x_lag = 0.2. The distance from the midplane is shown.

At stiffness La = 100 stable particle pairs only occur for x_lag= −0.2 and −0.4. The pronounced bounded oscillations have vanished although modulations with small amplitudes are visible. All other B2 pairs are unstable and axially move away from each other while reaching their equilibrium positions with finite distance from the center. An exception is the case with x_lag = 0. The leading particle is strongly dragged to the center by the biconcave particle and ultimately reaches the center, which we already identified as an additional stable lateral position when discussing the dynamics of the M2 pairs.

We now discuss the bi-dispersed pairs B3 and B4 where leading and lagging particles are spherical and have different stiffnesses. In Fig. 11 we plot lateral positions and axial distances versus time. When the leading particle is softer (La = 10) than the lagging particle (La = 100), the B3 pairs with x_lag = 0 and −0.2 form stable pairs. The leading soft particle succeeds to drag the stiffer particle to the channel center, where it occupies the additional stable position so that the pair can stay bounded. In the remaining cases the lagging particle migrates towards the off-center stable position, which ultimately makes the B3 pairs unstable.

In the B4 pairs the rigidity values are switched between the particles. As a result both leading and lagging particles move towards the channel center. A special case occurs when both particles start from the symmetric lateral positions x_lead = −x_lag = 0.2. They obstruct each other from aligning at the channel center, which results in oscillations in the lateral positions. However, the axial distance between the particles only fluctuates little about its initial value and the pair is stable. In the rest of the cases the leading and lagging particles successfully align at the channel center. They drift apart at different axial velocities due to their different equilibrium shapes (see Fig. 5), which directly affect the drag acting on each particle in the axial direction. Hence, even if the leading and lagging particles have the same lateral position, their axial distance does not remain bounded. As shown in Fig. 11, for the unstable pairs the axial distance increases at the same rate once both particles are at the channel center. We did not simulate the unstable pairs for longer times because it is computationally expensive as the axial distance grows very slowly. On extrapolating the trend of the curve for x_lag = 0.2, one would need ≈400 time units for Δz to reach Δz = 10. We also simulated the B4 pairs with the starting position of the leading particle x_lead = 0. However, we did not find any stable configurations in these cases.

4 Concluding remarks

Inertia driven separation and focusing of soft particles in microfluidics has direct relevance for biomedical applications. Inertial microfluidics has recently emerged as a robust and precise technique for particle manipulation since it relies entirely on internal hydrodynamic forces. Moreover, inertial microfluidics functions at finite Reynolds numbers, which allows faster processing of samples. In the present article we investigated the hydrodynamics of a particle pair in inertial microfluidics using 3D numerical simulations. We combined the lattice Boltzmann method to compute the fluid flow with the finite element method to model particle dynamics. To couple fluid flow and particle motion, the immersed boundary method was used.

Fig. 13 The equilibrium axial distance Δz_eq between the leading and lagging particle is plotted versus x_lag for stable particle pairs.

In this work we considered mono- (M1–M3) and bi-dispersed (B1–B4) particle pairs of different softness and particle shape (see Fig. 2). For each pair we investigated the influence of different starting positions of the lagging particle and identified four types: stable pair, stable pair with damped oscillations, stable pair with bounded oscillations, and unstable pair.

We found that the starting position of the lagging particle does not affect the stability of the mono-dispersed pairs with stiffness La = 10 (M1). On increasing the stiffness to La = 100 (M2), the mono-dispersed pairs become unstable when the leading and lagging particles start in the same channel half. Moreover, when the lagging particle starts at x_lag = 0, the pair undergoes damped oscillations. The oscillations arise due to hydrodynamic interactions between particles, while inertial forces damp them. The maximum amplitude of the damped oscillations decreases as we shift the starting position of the lagging particle away from the leading particle. Finally, from mono-dispersed pairs of perfectly solid particles (M3) only the pair with initially symmetric lateral positions about the centerline retains its stability.

When we change the shape of the lagging particle from spherical to biconcave, bi-dispersed pairs with particle stiffness La = 10 are always stable. They show bounded oscillations in the axial distance due to the tumbling motion of the biconcave particle. For La = 100 (B2) the biconcave shape does not change the qualitative behavior when compared to the mono-dispersed pairs (M2) of the same particle stiffness.

Finally, the soft leading (La = 10) and stiff lagging (La = 100) particles in the bi-dispersed pair B3 assemble at the channel center. The pair is stable when the lagging particle starts at x_lag = 0 or −0.2, while the remaining starting positions result in unstable pairs. When the stiffness of the leading and lagging particles are interchanged (pair B4), both particles migrate towards the channel center irrespective of the starting position and they ultimately drift apart due to their different shapes. Only for symmetric starting positions about the channel center a stable pair with lateral oscillations occur.

Our work explores and classifies different two-particle motions, which should be observable in the flow of mono- and bi-dispersed particle suspensions in the regime of inertial microfluidic. To explicitly test the different types of pair dynamics, one could use microchannels with inlets, which place both particles at specific lateral positions.

Finally, in inertial microfluidics the multi-particle states consist of linear trains either moving on the same streamline or in a staggered configuration as observed in experiments.^96,101 As we demonstrated in ref. 61 for rigid particles, knowing the bounded two-particle states one can directly infer the lattice constant of the staggered trains, while the lattice constant of linear trains on the same streamline typically have a value determined by the range of the two-particle interaction. In this context, the steady-state value of the axial distance between the leading and lagging particles is summarized in Fig. 13 for mono- and bi-dispersed pairs, when they are stable. Thus, the present study will help to understand the multi-particle train configurations occurring in inertial microfluidics for soft particles and their mixtures, which will be more involved compared to the case of ref. 61, where only rigid particles of the same size were studied.

Conflicts of interest

There are no conflicts of interest to declare.

Appendix A: benchmarking of the present LBFEIBM solver

We briefly present the validation of the present coupled lattice Boltzmann finite element immersed boundary method (LBFEIBM) for soft and solid particles. To benchmark the modeling of soft particles, we look at the inertial migration of a deformable capsule of radius a = 0.2 and softness La = 10 in a quadratic microchannel at Reynolds number Re = 100. In the present case, the deformable capsule is initialized on the channel diagonal. We use 120 lattice nodes to discretize the channel width, while the surface of the capsule is discretized using 5120 elements. Time-wise variation of the capsule position along the channel diagonal is shown in Fig. 14. The dashed black line in Fig. 14 indicates the equilibrium position reported by Schaaf and Stark.²⁰

Fig. 14 Inertial migration of a deformable capsule with La = 10 along the diagonal of a square microchannel at Reynolds number Re = 100. The dashed line shows the equilibrium position reported in ref. 20. The inset shows the equilibrium shape of the capsule.

Lastly, to validate the implementation of solid particles, we demonstrate damped oscillations of the axial distance in a pair of solid particles. For this, we consider a pair of solid particles in the midplane of a rectangular microchannel with an aspect ratio 0.5. We fix the particle radius a, initial position (x_lead, x_lag), initial axial distance Δz, and Reynolds number Re to 0.4, (0.24, −0.2), 2, and 10, respectively. We use 75 lattice nodes to discretize the channel width. The surface of the solid particles is approximated using 20480 [Doublestruck P]₁ elements. The temporal evolution of the ratio Δz/a is plotted in Fig. 15, along with its steady-state value reported by Schaaf et al.⁶⁰

Fig. 15 Damped oscillations of the axial distance Δz in a pair of solid particles with the initial position (x_lead, x_lag) = (0.24, −0.2). The dashed line shows the equilibrium distance reported in ref. 60.

Acknowledgements

We thank C. Schaaf for sharing the source code. We acknowledge support from the Deutsche Forschungsgemeinschaft in the framework of the Collaborative Research Center SFB 910.

References

1. G. Y. Lee and C. T. Lim, Trends Biotechnol., 2007, 25, 111–118.
2. S. Suresh, J. Spatz, J. Mills, A. Micoulet, M. Dao, C. Lim, M. Beil and T. Seufferlein, Acta Biomater., 2005, 1, 15–30.
3. S. Suresh, Acta Biomater., 2007, 3, 413–438.
4. A. van de Stolpe, K. Pantel, S. Sleijfer, L. W. Terstappen and J. M. den Toonder, Cancer Res., 2011, 71, 5955–5960.
5. X. Cheng, D. Irimia, M. Dixon, K. Sekine, U. Demirci, L. Zamir, R. G. Tompkins, W. Rodriguez and M. Toner, Lab Chip, 2007, 7, 170–178.
6. G. M. Whitesides, Nature, 2006, 442, 368–373.
7. E. K. Sackmann, A. L. Fulton and D. J. Beebe, Nature, 2014, 507, 181–189.
8. C. Haber, Lab Chip, 2006, 6, 1118–1121.
9. A. Karimi, S. Yazdi and A. M. Ardekani, Biomicrofluidics, 2013, 7, 021501.
10. D. G. Grier, Nature, 2003, 424, 810–816.
11. B. Çetin and D. Li, Electrophoresis, 2011, 32, 2410–2427.
12. T. P. Forbes and S. P. Forry, Lab Chip, 2012, 12, 1471.
13. M. Evander, L. Johansson, T. Lilliehorn, J. Piskur, M. Lindvall, S. Johansson, M. Almqvist, T. Laurell and J. Nilsson, Anal. Chem., 2007, 79, 2984–2991.
14. A. Hochstetter, R. Vernekar, R. H. Austin, H. Becker, J. P. Beech, D. A. Fedosov, G. Gompper, S.-C. Kim, J. T. Smith, G. Stolovitzky, J. O. Tegenfeldt, B. H. Wunsch, K. K. Zeming, T. Krüger and D. W. Inglis, ACS Nano, 2020, 14, 10784–10795.
15. D. D. Carlo, Lab Chip, 2009, 9, 3038.
16. J. M. Martel and M. Toner, Annu. Rev. Biomed. Eng., 2014, 16, 371–396.
17. J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M. E. Warkiani and W. Li, Lab Chip, 2016, 16, 10–34.
18. S. Razavi Bazaz, A. Mashhadian, A. Ehsani, S. C. Saha, T. Krüger and M. Ebrahimi Warkiani, Lab Chip, 2020, 20, 1023–1048.
19. C. Prohm and H. Stark, Lab Chip, 2014, 14, 2115.
20. C. Schaaf and H. Stark, Soft Matter, 2017, 13, 3544–3555.
21. G. Segré and A. Silberberg, Nature, 1961, 189, 209–210.
22. G. Segré and A. Silberberg, J. Fluid Mech., 1962, 14, 136–157.
23. B. P. Ho and L. G. Leal, J. Fluid Mech., 1974, 65, 365–400.
24. E. S. Asmolov, J. Fluid Mech., 1999, 381, 63–87.
25. J.-P. Matas, J. F. Morris and É. Guazzelli, J. Fluid Mech., 2009, 621, 59–67.
26. K. Hood, S. Lee and M. Roper, J. Fluid Mech., 2015, 765, 452–479.
27. D. D. Carlo, J. F. Edd, K. J. Humphry, H. A. Stone and M. Toner, Phys. Rev. Lett., 2009, 102, 094503.
28. X. Shao, Z. Yu and B. Sun, Phys. Fluids, 2008, 20, 103307.
29. I. Lashgari, M. N. Ardekani, I. Banerjee, A. Russom and L. Brandt, J. Fluid Mech., 2017, 819, 540–561.
30. H. T. Kazerooni, W. Fornari, J. Hussong and L. Brandt, Phys. Rev. Fluids, 2017, 2, 084301.
31. J.-P. Matas, J. F. Morris and É. Guazzelli, J. Fluid Mech., 2004, 515, 171–195.
32. A. A. S. Bhagat, S. S. Kuntaegowdanahalli and I. Papautsky, Phys. Fluids, 2008, 20, 101702.
33. K. Miura, T. Itano and M. Sugihara-Seki, J. Fluid Mech., 2014, 749, 320–330.
34. C. Liu, G. Hu, X. Jiang and J. Sun, Lab Chip, 2015, 15, 1168–1177.
35. J. Feng, H. H. Hu and D. D. Joseph, J. Fluid Mech., 1994, 277, 271–301.
36. H. Brenner, Chem. Eng. Sci., 1961, 16, 242–251.
37. P. G. Saffman, J. Fluid Mech., 1965, 22, 385–400.
38. D. D. Carlo, D. Irimia, R. G. Tompkins and M. Toner, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 18892–18897.
39. S. C. Hur, N. K. Henderson-MacLennan, E. R. B. McCabe and D. D. Carlo, Lab Chip, 2011, 11, 912.
40. C. K. W. Tam and W. A. Hyman, J. Fluid Mech., 1973, 59, 177–185.
41. F. Risso, F. Collé-Paillot and M. Zagzoule, J. Fluid Mech., 2006, 547, 149.
42. S. K. Doddi and P. Bagchi, Int. J. Multiphase Flow, 2008, 34, 966–986.
43. G. Coupier, B. Kaoui, T. Podgorski and C. Misbah, Phys. Fluids, 2008, 20, 111702.
44. L. Shi, T.-W. Pan and R. Glowinski, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 056308.
45. D. S. Hariprasad and T. W. Secomb, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 033008.
46. P. Bagchi, Biophys. J., 2007, 92, 1858–1877.
47. T. Krüger, Rheol. Acta, 2015, 55, 511–526.
48. K. Vahidkhah, P. Balogh and P. Bagchi, Sci. Rep., 2016, 6, 28194.
49. P. Balogh and P. Bagchi, Physiol. Rep., 2019, 7, e14067.
50. A. Kilimnik, W. Mao and A. Alexeev, Phys. Fluids, 2011, 23, 123302.
51. S. J. Shin and H. J. Sung, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 046321.
52. Y. Kim and M.-C. Lai, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 066321.
53. D. Salac and M. J. Miksis, J. Fluid Mech., 2012, 711, 122–146.
54. S. J. Shin and H. J. Sung, Int. J. Heat Fluid Flow, 2012, 36, 167–177.
55. Z. Y. Luo, S. Q. Wang, L. He, F. Xu and B. F. Bai, Soft Matter, 2013, 9, 9651.
56. Y.-L. Chen, RSC Adv., 2014, 4, 17908–17916.
57. Z. Wang, Y. Sui, A.-V. Salsac, D. Barthès-Biesel and W. Wang, J. Fluid Mech., 2016, 806, 603–626.
58. D. Alghalibi, M. E. Rosti and L. Brandt, Phys. Rev. Fluids, 2019, 4, 104201.
59. Y. Yan, J. F. Morris and J. Koplik, Phys. Fluids, 2007, 19, 113305.
60. C. Schaaf, F. Rühle and H. Stark, Soft Matter, 2019, 15, 1988–1998.
61. C. Schaaf and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2020, 43, 50.
62. J. F. Morris, J. Fluid Mech., 2016, 792, 1–4.
63. E. Lac, A. Morel and D. Barthès-Biesel, J. Fluid Mech., 2007, 573, 149–169.
64. O. Aouane, A. Farutin, M. Thiébaud, A. Benyoussef, C. Wagner and C. Misbah, Phys. Rev. Fluids, 2017, 2, 063102.
65. S. K. Doddi and P. Bagchi, Int. J. Multiphase Flow, 2008, 34, 375–392.
66. D. Leighton and A. Acrivos, J. Fluid Mech., 1987, 177, 109–131.
67. H. Lan and D. B. Khismatullin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 012705.
68. T. Krüger, B. Kaoui and J. Harting, J. Fluid Mech., 2014, 751, 725–745.
69. T. Krüger, F. Varnik and D. Raabe, Comput. Math. Appl., 2011, 61, 3485–3505.
70. H. Bruus, Theoretical microfluidics, Oxford University Press, Oxford New York, 2008.
71. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Clarendon Press Oxford University Press, Oxford New York, 2001.
72. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer, Berlin, Heidelberg, 2016.
73. T. Krüger, S. Frijters, F. Günther, B. Kaoui and J. Harting, Eur. Phys. J.-Spec. Top., 2013, 222, 177–198.
74. Z. Guo, C. Zheng and B. Shi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 046308.
75. J. Latt, PhD thesis, University of Geneva, 2007.
76. J. Latt, B. Chopard, O. Malaspinas, M. Deville and A. Michler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 056703.
77. J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani, D. Lagrava, F. Brogi, M. B. Belgacem, Y. Thorimbert, S. Leclaire, S. Li, F. Marson, J. Lemus, C. Kotsalos, R. Conradin, C. Coreixas, R. Petkantchin, F. Raynaud, J. Beny and B. Chopard, Comput. Math. Appl., 2021, 81, 334–350.
78. H. Elman, D. Silvester and A. Wathen, Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics, Oxford University Press, 2014.
79. R. Skalak, A. Tozeren, R. Zarda and S. Chien, Biophys. J., 1973, 13, 245–264.
80. W. Helfrich, Z. Naturforsch., C: J. Biosci., 1973, 28, 693–703.
81. G. Gompper and D. M. Kroll, J. Phys. I, 1996, 6, 1305–1320.
82. U. Seifert, Adv. Phys., 1997, 46, 13–137.
83. J.-M. Charrier, S. Shrivastava and R. Wu, J. Strain Anal. Eng. Des., 1989, 24, 55–74.
84. T. Krüger, Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear, Springer Spektrum, Wiesbaden, 2012.
85. T. Inamuro, Fluid Dyn., 2012, 44, 024001.
86. K. Suzuki and T. Inamuro, Comput. Fluids, 2011, 49, 173–187.
87. Z.-G. Feng and E. E. Michaelides, Comput. Fluids, 2009, 38, 370–381.
88. C. S. Peskin, Acta Numer., 2002, 11, 479–517.
89. R. Mittal and G. Iaccarino, Annu. Rev. Fluid Mech., 2005, 37, 239–261.
90. S. S. Jain, K. Kamrin and A. Mani, J. Comput. Phys., 2019, 399, 108922.
91. Z. Wang, J. Fan and K. Luo, Int. J. Multiphase Flow, 2008, 34, 283–302.
92. S. K. Kang and Y. A. Hassan, Int. J. Numer. Methods Fluids, 2010, 66, 1132–1158.
93. E. Evans and Y.-C. Fung, Microvasc. Res., 1972, 4, 335–347.
94. S. Ramanujan and C. Pozrikidis, J. Fluid Mech., 1998, 361, 117–143.
95. K. Wang, X. H. Sun, Y. Zhang, T. Zhang, Y. Zheng, Y. C. Wei, P. Zhao, D. Y. Chen, H. A. Wu, W. H. Wang, R. Long, J. B. Wang and J. Chen, R. Soc. Open Sci., 2019, 6, 181707.
96. W. Lee, H. Amini, H. A. Stone and D. D. Carlo, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 22413–22418.
97. K. Hood and M. Roper, Phys. Rev. Fluids, 2018, 3, 094201.
98. A. Z. K. Yazdani and P. Bagchi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 026314.
99. D. A. Fedosov, M. Peltomäki and G. Gompper, Soft Matter, 2014, 10, 4258–4267.
100. B. Kaoui, G. Biros and C. Misbah, Phys. Rev. Lett., 2009, 103, 188101.
101. S. Kahkeshani, H. Haddadi and D. D. Carlo, J. Fluid Mech., 2016, 786, R3.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sm00276g

---