Numerical investigation of the formation and stability of homogeneous pairs of soft particles in inertial microfluidics

Abstract

We investigate the formation and stability of a pair of identical soft capsules in channel flow under mild inertia. We employ a combination of the lattice Boltzmann, finite element and immersed boundary methods to simulate the elastic particles in flow. Validation tests show excellent agreement with numerical results obtained by other research groups. Our results reveal new trajectory types that have not been observed for pairs of rigid particles. While particle softness increases the likelihood of a stable pair forming, the pair stability is determined by the lateral position of the particles. A key finding is that stabilisation of the axial distance occurs after lateral migration of the particles. During the later phase of pair formation, particles undergo damped oscillations that are independent of initial conditions. These damped oscillations are driven by a strong hydrodynamic coupling of the particle dynamics, particle inertia and viscous dissipation. While the frequency and damping coefficient of the oscillations depend on particle softness, the pair formation time is largely determined by the initial particle positions: the time to form a stable pair grows exponentially with the initial axial distance. Our results demonstrate that particle softness has a strong impact on the behaviour of particle pairs. The findings could have significant ramifications for microfluidic applications where a constant and reliable axial distance between particles is required, such as flow cytometry.

1. Introduction

Microfluidic devices play an increasingly important role in disease diagnostics. Due to their small footprint, high portability, relatively low cost and ever-improving manufacturing techniques, microfluidic devices have the potential to revolutionise point-of-care applications. Microfluidic devices exploit physical effects, such as flow and cell dynamics, at a smaller scale than conventionally sized devices. For example, cells can be manipulated into forming regularly spaced pairs and trains that are required for applications such as flow cytometry (Hur, Tse & Di Carlo 2010) and cell encapsulation (Moon et al. 2018).

Inertial microfluidics (IMF) is a relatively new research field that emerged in the late 2000s (Di Carlo et al. 2007; Russom et al. 2009). While conventional microfluidic devices are mostly operated in the limit of small Reynolds number, the Reynolds number in IMF is typically of the order of $10$–$100$ due to relatively high flow speeds. Besides an increased throughput, IMF features additional physical effects that can be exploited for particle manipulation, focusing in the channel cross-section, and separation. At finite inertia, particles and cells experience additional lift forces which lead to cross-streamline migration of particles that would not migrate in the Stokes limit (Ho & Leal 1974; Schonberg & Hinch 1989; Asmolov 1999; Matas, Morris & Guazzelli 2009). Inertial effects can also lead to axial ordering of particles (Matas et al. 2004; Lee et al. 2010), an effect that can be exploited to overcome the Poisson statistics that is often plaguing applications in the non-inertial regime (Lagus & Edd 2013). The ultimate aim is to bring IMF to maturity and routinely use it for diagnostic applications, such as recovery of rare cells from blood (Tanaka et al. 2012), separation of particles by deformability (for instance, diseased red blood cells from healthy ones) or search for sepsis markers (Gossett et al. 2012).

When the Reynolds number is of the order of $10$–$100$, rigid spherical particles in tube flow migrate to a radial equilibrium position approximately 60 % away from the tube centreline. This so-called Segré–Silberberg effect (Segre & Silberberg 1961) is caused by a balance of inertial shear-gradient lift forces, pushing the particles towards the wall, and wall repulsion forces caused by an increased pressure between the particles and the wall (Ho & Leal 1974; Schonberg & Hinch 1989; Asmolov 1999; Matas et al. 2009). The resulting lateral motion of the particles towards their equilibrium position is termed inertial migration. Altering the channel cross-section (Kim et al. 2016) or curving the channel into serpentine (Zhang et al. 2014) or spiral (Warkiani et al. 2014) geometries changes the number and location of lateral equilibrium positions and can accelerate lateral focusing of particles (Martel & Toner 2014).

Soft particles give rise to more complex behaviour than rigid particles since particle softness alone leads to cross-streamline migration, even in the non-inertial limit (Chen 2014). Additionally, the analysis of soft particle dynamics is more difficult due to the changing particle shape. The inertial migration of a single soft particle has been investigated extensively in a variety of different flow conditions, including shear flows (Ma et al. 2019), channel flows (Coclite et al. 2020), stratified flows (Jyothi, Renganathan & Pushpavanam 2019) and viscoelastic flows (Ni & Jiang 2020). Hur et al. (2011) identified experimentally that softer particles migrate to equilibrium positions closer to the channel centre. Kilimnik, Mao & Alexeev (2011) also demonstrated that the equilibrium position of soft particles is essentially independent of Reynolds number, a result confirmed by Schaaf & Stark (2017).

The situation gets even more interesting when particles are sufficiently close and interact hydrodynamically. In IMF stable pairs and trains of particles can form under some circumstances, which brings up important implications for particle manipulation, focusing and sorting. Lee et al. (2010) first identified the self-assembly of particle pairs in inertial flows through the mechanism of reversing streamlines. Since then, a number of studies have investigated the pair/train formation further, both numerically and experimentally. For stable pairs, axial distance has been shown to be independent of initial position for rigid particle pairs (Humphry et al. 2010; Schaaf & Stark 2020), while different trajectory types have been identified, depending on the interaction between the particles. Of these identified trajectory types, only damped oscillation trajectories have been shown to lead to stable pairs (Lan & Khismatullin 2014; Schaaf, Rühle & Stark 2019). Kahkeshani, Haddadi & Di Carlo (2016) observed the existence of two different axial equilibrium distances between rigid pairs, with the preference to each being Reynolds number dependent. Importantly, Schaaf & Stark (2020) identified that during the formation of particle trains, stable particle pairs first form and then group together to build particle trains. Patel & Stark (2021) investigated the effect of pair softness and shape for mono- and bi-disperse pairs, finding that increased particle softness leads to increased pair stability. However, the formation mechanisms of particle pairs and the conditions leading to stable pairs are still not well understood. Moreover, the distinct properties of soft particle pairs, compared with rigid particle pairs, need to be explored further to improve the design of inertial microfluidic devices that are used for soft cells.

In this paper we investigate the dynamics of a pair of identical soft capsules in straight channel flow at moderate inertia via immersed-boundary-lattice-Boltzmann-finite-element (IB-LB-FE) simulations (§ 2). We validate the model by simulating the interaction of a pair of soft capsules in a simple shear flow and the lateral migration of a single soft capsule in channel flow (§ 3). Our results (§ 4) show that pairs of identical soft capsules exhibit six different trajectory types, two of which have not been observed for rigid pairs. We demonstrate that the lateral equilibrium positions of the particles largely determine the stability of the pair and that the stabilisation of the axial distance between both particles occurs only after the lateral migration phase. These observations lead to the hypothesis that certain mechanisms contribute more strongly to the overall system: flow development as zeroth-order, single particle lateral migration as first-order, and axial spacing as second-order effects. We also observe two distinct phases of pair formation, an early axial approach phase and a later spiralling convergence phase that is largely independent of the initial phase. During the spiralling phase, particles are tightly coupled through hydrodynamic interactions. The spiralling dynamics is characterised by frequency and damping coefficients that are determined by particle softness, while the pair formation time is mostly dependent on the initial positions of the particles. We argue that our findings have strong implications for the understanding of particle interactions in mildly inertial flows and the design of inertial microfluidic devices aiming at the formation of regularly spaced particle pairs and trains (§ 5).

2. Physical and numerical model

The physical and numerical models are briefly outlined in §§ 2.1 and 2.2, respectively.

2.1. Physical model

2.1.1. Governing equations and physical parameters

We consider a single or two soft capsules flowing in a simple planar shear flow or in a straight channel. We assume an incompressible Newtonian liquid. The suspended particles are hyperelastic and neutrally buoyant capsules which are filled with the same liquid and are spherical in their undeformed state. While the liquid is governed by the incompressible Navier–Stokes equations, we employ two different elastic models for the capsules, either the Skalak model (Skalak et al. 1973)

(2.1)\begin{equation} w_{s} = \frac{\kappa_{s}}{12} (I_1^2 + 2 I_1 - 2 I_2) + \frac{\kappa_\alpha}{12} I_2^2 \end{equation}

or the neo-Hookean model

(2.2)\begin{equation} w_{s} = \frac{\kappa_{s}}{6} \left(I_1 - 1 + \frac{1}{I_2 + 1}\right), \end{equation}

where $w_{s}$ is the areal energy density, $I_1$ and $I_2$ are the in-plane strain invariants (Krüger, Varnik & Raabe 2011), and $\kappa _{s}$ and $\kappa _\alpha$ are the elastic shear and area dilation moduli. We include a membrane bending energy

(2.3)\begin{equation} w_{b} = \frac{\kappa_{b}}{2} \left(H - H^{(0)}\right)^2, \end{equation}

where $H$ and $H^{(0)}$ are the trace of the surface curvature tensor and the spontaneous curvature, respectively, and $\kappa _{b}$ is the bending modulus.

The two flow scenarios considered are (i) simple shear flow between two flat and rigid plates and (ii) force-driven flow in a straight and rigid channel with a rectangular cross-section. The no-slip boundary condition is assumed at the surfaces of the channel and the particles.

The relevant parameters are liquid density $\rho$ and kinematic viscosity $\nu$; the radius $a$ of the undeformed capsules; the elastic shear modulus $\kappa _{s}$, bending modulus $\kappa _{b}$ and area dilation modulus $\kappa _\alpha$ of the capsules; the channel half-width $w$, half-height $h$ and length $L$; and either shear rate $\dot {\gamma }$ for simple shear flow or maximum velocity $U_{max}$ at the channel centre for force-driven flow.

2.1.2. Dimensionless groups

The particle Reynolds number ${Re}_{p}$ is used for the simple shear flow cases and as defined by Doddi & Bagchi (2008),

(2.4)\begin{equation} {Re}_{p} = \frac{\dot{\gamma} a^2}{\nu}. \end{equation}

For the channel flow cases, the channel Reynolds number ${Re}_{c}$ is used. We adopt the definition of Schaaf & Stark (2017), Schaaf et al. (2019),

(2.5)\begin{equation} {Re}_{c} = \frac{U_{max} w}{\nu}. \end{equation}

The capillary number $Ca$ is the ratio of the viscous stress of the liquid to the characteristic elastic shear stress of the capsule membrane,

(2.6)\begin{equation} {Ca} = \frac{\rho \nu \dot{\gamma} a}{\kappa_{s}}. \end{equation}

While the capillary number depends on the flow field, the Laplace number $La$ is a combination of material properties only and suitable to isolate the contribution of particle softness in inertial flows (Schaaf & Stark 2017). The Laplace number is defined as the ratio between the particle Reynolds number and the capillary number,

(2.7)\begin{equation} {La} =\frac{{Re}_{{p}}}{{Ca}} = \frac{\kappa_{s} a}{\rho \nu^2}. \end{equation}

Other dimensionless groups are the confinement ratio $\chi = a/h$ or $\chi = a/w$, depending on whether $h$ or $w$ is smaller, the channel aspect ratio $\alpha = w/h$, the reduced dilation modulus $\tilde {\kappa }_\alpha = \kappa _\alpha / \kappa _{s}$ and the reduced bending modulus $\tilde {\kappa }_{b} = \kappa _{b} / (\kappa _{s} a^2)$.

2.1.3. Characteristic scales

In order to non-dimensionalise reported distances (or positions) and times, we use characteristic length scales and a time scale. Depending on context, distances or locations are non-dimensionalised either by particle radius $a$, channel half-width $w$ or channel half-height $h$. For channel flow, time is non-dimensionalised by the advection time

(2.8)\begin{equation} t_{ad} = \frac{a}{U_{max}}. \end{equation}

2.2. Numerical model

The numerical model consists of a partitioned fluid-structure interaction solver in which the lattice Boltzmann (LB) method is used for the liquid, the finite element (FE) method for the capsule dynamics and the immersed boundary (IB) method for the fluid-structure interaction. This IB-LB-FE solver has previously been employed in the study of the dynamics of deformable red blood cells and capsules (Krüger et al. 2013; Krüger, Holmes & Coveney 2014). Here, we provide essential properties of the model, while comprehensive details are available elsewhere (Krüger et al. 2011).

For the LB method, we use the D3Q19 lattice (Qian, D'Humières & Lallemand 1992) and the BGK collision operator (Bhatnagar, Gross & Krook 1954) with relaxation time $\tau$. The viscosity of the liquid and the relaxation time satisfy

(2.9)\begin{equation} \nu = c_{s}^2 \left(\tau - \frac{\Delta t}{2}\right), \end{equation}

where $c_{s}$ is the lattice speed of sound and $\Delta t$ is the time step. For the D3Q19 lattice, $c_{s}^2 = \Delta x^2/(3 \Delta t^2)$ holds where $\Delta x$ is the lattice resolution. For the channel flow cases, flow is driven by a constant body force following the forcing method of Guo, Zheng & Shi (2002). This form of the LB method is widely used in the field of fluid dynamics, including in previous IMF studies (Schaaf & Stark 2017; Schaaf et al. 2019).

Each capsule is represented by a surface mesh with $N_{f}$ flat triangular faces (or elements) defined by three nodes (or vertices) each. At a given time step, the capsule mesh is generally deformed. The hyperelastic forces acting on each vertex are calculated as a function of the mesh deformation state through an explicit scheme. The shear and area dilation forces result from the deformation gradient tensor of each face, while the bending forces are related to the angles between normal vectors of pairs of neighbouring faces (Krüger et al. 2013).

We employ an IB method with a three-point stencil (Peskin 2002). The forces obtained from the FE scheme are spread from the Lagrangian mesh to the Eulerian lattice where they act on the surrounding fluid nodes through the LB algorithm. The updated fluid velocity is then interpolated at the location of each mesh node. The positions of the mesh nodes are updated using the forward-Euler method, assuming a massless membrane which is appropriate for neutrally buoyant capsules. This treatment recovers the no-slip boundary condition at the surface of the capsules and the momentum exchange between the liquid and the capsule membrane.

The no-slip boundary condition at the resting and moving walls for channel and shear flow, respectively, is realised by the standard half-way bounce-back condition (Ladd 1994). The flow is periodic along those directions that are not confined by walls (flow direction in the channel case, and shearing plane in the shear flow case). The channel length $L$ is chosen sufficiently long to avoid the interaction of capsules with their periodic images. In all simulations the hydrodynamic forces are sufficient to prevent contact between capsules – capsules do not come closer than approximately $2 \Delta x$ – and an artificial repulsion force between capsules is not required. Likewise, capsules always keep a large distance from the confining walls due to hydrodynamic lift, and an artificial capsule-wall repulsion force is not needed.

3. Benchmark tests

We test our model by comparing simulation results obtained from our solver with previously published results from other groups: trajectories of a soft particle pair in shear flow (§ 3.1) and lateral migration of a single soft capsule in channel flow (§ 3.2).

3.1. Soft particle pair in shear flow

Doddi & Bagchi (2008) numerically investigated the effect of inertia on the interaction between a soft particle pair in a shear flow, whereby the particles are slightly offset from the centre of the channel as shown in figure 1. By increasing the shear rate $\dot {\gamma }$, the trajectories of the pair switch from passing to reversing. Doddi & Bagchi employed the neo-Hookean model, (2.2), for the capsules and the front-tracking/IB method proposed by Unverdi & Tryggvason (1992). Here, we reproduce these results with our solver to test the accuracy of the IB-LB-FE method and its ability to capture hydrodynamic interactions of soft capsules in mildly inertial flows. Table 1 shows the relevant parameters of this benchmark case.

Figure 1. Schematic of the first benchmark case: soft particle pair in shear flow. (a) Both particles are located on a plane parallel to the shear direction. (b) The shear rate is defined by the speed of the moving walls and their separation: $\dot {\gamma } = 2 u_{w} / L$. Both particles have the same distance from the walls, where one particle is closer to the bottom and the other closer to the top wall. Particles are initially separated by $\delta x_0$ and $\delta z_0$ along the $x$- and $z$-axes, respectively. Simulation parameters are reported in table 1.

Table 1. Parameters of the first benchmark case: soft particle pair in shear flow. See figure 1 for an illustration of the set-up. The shear rate $\dot {\gamma }$ depends on ${Re}_{p}$ according to (2.4), and the shear elasticity $\kappa _{s}$ is obtained from (2.6). The liquid density is set to $1$ in simulation units.

While Doddi & Bagchi (2008) did not report whether they used a bending resistance, the work originally employing the numerical method included a finite bending resistance to avoid folding of the membrane (Eggleton & Popel 1998). Therefore, we ran two sets of simulations: one without bending resistance, the other with a finite bending resistance ($\tilde {\kappa }_{b} = 0.00287$). Figure 2(a) shows our simulated particle trajectories for different values of ${Re}_{p}$ in the plane highlighted in figure 1 in comparison with the results of Doddi & Bagchi (2008). Good agreement is seen between the original results and our simulations that include a bending resistance, while simulations without bending resistance show some differences in the later stages of the passing trajectories. Independent of the bending resistance, we see that, as ${Re}_{p}$ is increased, each particle moves closer to the mid-plane between the walls, resulting in the transition from passing to reversing trajectories. Using the IB-LB-FE model, we observed this transition between ${Re}_{p} = 0.375$ and $0.575$, in agreement with previous results.

Figure 2. Trajectories of particle pairs in simple shear flow for various Reynolds numbers. Blue: ${Re}_{p} = 0.125$; orange: ${Re}_{p} = 0.375$; green: ${Re}_{p} = 0.575$; red: ${Re}_{p} = 0.75$. Black circles indicate the initial position of each particle. See figure 1 for the geometry set-up and table 1 for simulation parameters. (a) Comparison of our results for the neo-Hookean model with and without bending resistance with previous data (Doddi & Bagchi 2008). Squares are data points extracted from figure 8 in Doddi & Bagchi (2008) using WebPlotDigitizer v4.4. (b) Comparison between neo-Hookean and Skalak models with bending resistance obtained from our IB-LB-FE solver. (a) Comparison of our simulations with original results. (b) Comparison between neo-Hookean and Skalak models.

Doddi & Bagchi (2008) employed the neo-Hookean model only. Since other elastic capsule models are often used, we have repeated the same test for capsules equipped with the Skalak model to investigate the sensitivity of the trajectories to the details of the capsule model. The reduced bending modulus is $\tilde {\kappa }_b = 0.00287$ for both elastic models, and $\tilde {\kappa }_\alpha = 2$ for the Skalak model (note that $\tilde {\kappa }_\alpha$ is not defined for the neo-Hooken model). Figure 2(b) shows the comparison of trajectories for the neo-Hookean and the Skalak models as obtained from the IB-LB-FE solver. Overall, details of the trajectories are not significantly altered by the model. The Skalak model leads to slightly larger lateral displacements. This difference is probably caused by the reduced deformation of the capsules due to the strain-hardening properties of the Skalak model compared with the strain-softening neo-Hookean model. While the difference between the constitutive models is not trivial, the trajectories are qualitatively similar, and the transition between passing and reversing trajectories is also between ${Re}_{p} = 0.375$ and $0.575$ for the Skalak model.

3.2. Lateral migration of a single soft particle in channel flow

Before investigating the interaction of a particle pair in channel flow with moderate inertia, it is crucial to ensure the inertial effects on a single particle are captured accurately. We consider the inertial migration of a single soft capsule in a pressure-driven Poiseuille flow through a square duct with edge length $2w$ and axial length $L = 8w$, as shown in figure 3. Periodic boundary conditions are used along the flow axis.

Figure 3. Schematic of the first benchmark case: soft particle migration in channel flow. (a) The arrow indicates the flow direction. The grey plane indicates the channel cross-section. (b) The particle is initially located away from the channel centreline. Simulation parameters are reported in table 2.

Table 2. Parameters of the second benchmark case: soft particle migration in channel flow. See figure 3 for an illustration of the set-up. The channel Reynolds number is varied by the body force and, therefore, $U_{max}$ via (2.5), and the Laplace number is controlled by the shear elasticity via (2.7). The liquid density is set to $1$ in simulation units.

Schaaf & Stark (2017) investigated the effect of Laplace number on the lateral equilibrium position of a soft capsule at a given channel Reynolds number in the same geometry. They demonstrated that the lateral equilibrium position moves away from the channel centre as $La$ increases. For the case ${Re}_{c} = 100$, a transition occurs where the equilibrium position switches from a diagonal location to a channel face centre. As a hallmark of lateral migration in a square duct, one can generally distinguish two phases of migration: a first phase during which the capsule quickly migrates in the radial direction, and a second phase defined by a slower circumferential migration along a heteroclinic orbit (Nakagawa et al. 2015). This behaviour is caused by the shear gradients of the velocity field being more pronounced in the radial than in the circumferential direction.

We have run simulations for ${Re}_{c} = 10$ and $100$, and for $La$ between $1$ and $100$. The Skalak model was employed for the capsule membrane. Table 2 summarises the simulation parameters. The initial position of the capsule is the same for each set of parameters considered. We compare cross-sectional particle trajectories obtained from our solver with those results reported by Schaaf & Stark (2017). Figure 4 shows that the capsule trajectories and lateral equilibrium positions generally agree well. In the cases where results do not agree well, the equilibrium positions recovered by our solver are located either on the channel diagonals or midway along channel faces, as expected from symmetry considerations, whereas the trajectories of Schaaf et al. appear to have stopped short. This discrepancy could be explained by the runtime of the simulations by Schaaf & Stark (2017) who mentioned that not all of the particles ‘reach their equilibrium position on the diagonal or the main axis during the simulations’.

Figure 4. Cross-sectional trajectories and lateral equilibrium positions of a single capsule in channel flow for various values of $La$ at (a) ${Re}_{c} = 10$ and (b) ${Re}_{c} = 100$. Our results are indicated by solid lines with circles, and dashed lines with squares show the results reported by Schaaf & Stark (2017). Grey symbols mark the initial capsule position, and colourful symbols indicate the final positions on the channel cross-section. The channel centreline is located at $(y,z) = (0, 0)$. Dashed grey lines are guides for the eyes.

We also noticed some oscillations in the results of Schaaf et al. at ${La} = 5$ and $10$ that appear to affect the trajectories during the second phase of migration. Our simulations did not show these oscillations, which could be due to a different bending model employed.

Overall, our IB-LB-FE solver produces results consistent with those of previously validated solvers and has been shown to be suitable for the investigation of inertial migration and the hydrodynamic interaction of soft capsules.

4. Results and discussion

Previous works have investigated the formation and stability of pairs (Gupta et al. 2018; Udono 2020) and trains (Kahkeshani et al. 2016; Hu et al. 2020) of rigid particles in channel flow. Schaaf et al. (2019) identified four different types of the trajectories of two rigid particles, depending on their initial positions in the channel. Patel & Stark (2021) also observed each of these trajectory types for soft particles in mono- and bi-disperse pairs. In the following we analyse the interaction of a pair of equally soft capsules for different Laplace numbers. In § 4.1 we define the cases investigated and parameters used, while § 4.2 presents the trajectory types that occur at a given Laplace number. Section 4.3 investigates the relationships of the lateral equilibrium position of each particle and the axial equilibrium distance between both particles with Laplace number and initial positions. Section 4.4 analyses the formation of stable pairs with a view to understanding how Laplace number and initial position effect focusing time and distance.

4.1. Case definition

Originally proposed by Schaaf et al. (2019), two particles are placed in a pressure-driven flow through a rectangular duct with width $2 w$ and height $2 h$ with aspect ratio $w/h = 2$ (figure 5). The length of the channel is $L = 10 h$, and the flow is periodic along the flow direction.

Figure 5. Schematic of particle pairs in a rectangular duct. Parameter values are given in § 4.1. (a) The flow is along the $x$-axis (blue arrow). Particles are initially located on the mid-plane with $y = \text {const.}$ (indicated by grey plane). (b) Depending on their initial position along the $x$-axis, we distinguish between the leading and the lagging particle.

Both capsules are initially located on the mid-plane between the side walls ($y = \text {const.}$), while the initial $x$- and $z$-coordinates are varied. We distinguish between the initially leading (farther downstream) and lagging (farther upstream) particles, according to their initial positions on the flow axis ($x$-axis). The limitation to the mid-plane is justified since particle equilibrium positions in high aspect-ratio channels at the Reynolds number investigated here are usually in the middle of the long channel edges (Prohm & Stark 2014). Particles initially located on this mid-plane will usually stay on this plane while moving along the $x$-axis and migrating along the $z$-axis. We have not observed particles leaving the mid-plane in any of our simulations. We expect that changing the channel aspect ratio $w/h$ will have only a minor influence on the results as long as the confinement by the side walls is lower than that by the top and bottom walls. Investigating the more general cases with arbitrary channel aspect ratios and initial particle positions not confined to the mid-plane would lead to an unmanageable number of free parameters and is beyond the scope of this paper.

In all following simulations we have employed the Skalak membrane model with a reduced dilation modulus $\tilde {\kappa }_\alpha = 2$ and reduced bending modulus $\tilde {\kappa }_{b} = 0.00287$. The channel Reynolds number is ${Re}_{c} = 10$ in all cases, and the numerical viscosity is kept at $\nu = 1/6\Delta x^2/\Delta t$. The channel dimensions are $2 w = 80 \Delta x$, $2 h = 160 \Delta x$ and $L = 560 \Delta x$. The undeformed capsule radius is $a = 16 \Delta x$. We focus on $La$, $\delta x_0$, $z_0^{lead}$ and $z_0^{lag}$ as free parameters (figure 5).

4.2. Interaction types of homogeneous soft particle pairs

We consider pairs of identical capsules (homogeneous pairs). First, we characterise the interaction types of the capsules before investigating the effects of Laplace number and initial particle positions in more detail.

In order to understand the general effect of softness and initial positions of the capsules, we simulated pairs at three different Laplace numbers: ${La} = 1$, $10$ and $100$. The investigated range of Laplace number was chosen to allow comparison with previous work (Schaaf & Stark 2017) while also ensuring numerical instabilities did not occur for more rigid particles. Softer particles were not considered since additional interesting behaviour was not observed beyond ${La} = 1$ for the considered confinement. The initial positions of both capsules along the height axis ($z$-axis) is varied with a total of 66 combinations (six initial positions of the leading particle: five in the top half of the channel and one at the centre; 11 initial positions of the lagging particle: five in the top half, five in the bottom half and one at the centre). Note that initial configurations mirrored at the $x$–$y$-plane lead to identical results. We also considered three different initial axial distances between the particle centres ($3a$, $5a$ and $7a$), leading to $3 \times 66 \times 3 = 594$ configurations in total. The actual number of simulations is $549$ since there are 45 symmetric cases when the leading particle is initially midway between the bottom and top walls.

Using rigid particles, Schaaf et al. (2019) identified four different particle interaction types: swap & scatter, pass & scatter, scatter and capture. Examples of each trajectory type are shown in figure 6, and brief descriptions are included in table 3. We changed the original names of the trajectory types that Schaaf et al. (2019) used in order to improve clarity when discussing the trajectories observed for soft particles. ‘Scatter’ refers to cases where the axial distance between particles eventually grows such that particles stop interacting with each other. ‘Capture’ indicates an interaction where the axial distance of the two particles is bound. Interestingly, we observe two capture subtypes: one where the axial (and also lateral) distance between the particles becomes constant after some time (‘stable’) and the other where the axial distance does not grow arbitrarily, but does not converge to a constant value either (‘partially stable’).

Figure 6. Typical snapshots at selected points in time and space for each trajectory type observed for soft particles. Cases (a–d) have been previously observed for rigid particles as well, while cases (e–f) have not. The initially leading/lagging particle is shown in red/blue, respectively. Note that the axial distances between different snapshots are not to scale. The six interaction types are characterised in more detail in table 3. (a) Swap & scatter. (b) Pass & scatter. (c) Capture. (d) Scatter. (e) Swap & capture. (f) Pass & capture.

Table 3. Descriptions of trajectory types of soft particle pairs with indication of whether this type is also observed for rigid particle pairs. Scatter means that the axial distance between particles grows until particles stop interacting with each other. Capture means that the axial distance is bound. Figure 6 visualises some example cases.

Using soft particles, we identified the same four previously observed trajectory types of rigid particles, and two new types that have not been found for rigid particles. We label the two new types as ‘swap & capture’ and ‘pass & capture’. Examples of the new types are shown in figure 6 and described in table 3. The newly observed partially stable particle pairs seem to occur only when particles are sufficiently soft, as will be discussed later.

Figure 7 shows the key results of our work: trajectory types resulting from all studied configurations (various initial lateral positions; initial axial particle distances of $3a$, $5a$ and $7a$; Laplace numbers ${La} = 1$, $10$ and $100$). On first inspection, a general trend can be seen that particle capture occurs under a wider range of initial conditions when the particles are softer (smaller $La$). Furthermore, particle capture is more likely when the initial axial distance $\delta x_0$ is smaller. For example, for the extreme case ${La} = 1$ and $\delta x_0 = 3a$, 90 % of all studied initial positions lead to capture. Only in cases where the leading particle is initially much closer to the channel centre (and, therefore, sufficiently faster) than the lagging particle, capturing does not occur.

Figure 7. Particle interaction types as a function of Laplace number, initial axial distance $\delta x_0$ and initial lateral positions $z_0$ for homogeneous pairs. The $x$- and $y$-axes of each panel indicate the lateral initial positions of the lagging and leading particles with respect to the channel centreline ($z_0 = 0$). The symbols indicate the interaction type as defined in table 3 and shown in figure 6. Symbols with a white background indicate a scattering trajectory. A dark grey background denotes a stable pair, while a light grey background indicates a partially stable pair. The dashed boxes define those configurations that are analysed in more detail in figure 8.

The increase of capture probability with decreasing $\delta x_0$ can be explained by the stronger hydrodynamic interaction of particles when they are initially closer. Additionally, particle softness is beneficial for capture since an initially lagging particle that is softer is able to increase its axial velocity more quickly than a more rigid particle by migrating closer towards the channel centre where the axial free-stream flow velocity is higher. Thus, a softer lagging particle has a higher chance of catching up with the leading particle before the leading particle moves away. We will take a closer look at the effect of Laplace number on particle dynamics in § 4.3.

Another notable effect of particle softness is that the diagrams in figure 7 become more symmetric with respect to $z_0 \leftrightarrow -z_0$ of the lagging particle. For ${La} = 100$ (more rigid particles), capture essentially only occurs when the leading and lagging particles are on different sides of the channel (positive $z_0$ for leading particle and negative $z_0$ for lagging particle). This effect can be explained by the known observation that pairs of rigid particles at high confinement $a/h$ are more stable when both particles are located on different sides of the channel centre (Patel & Stark 2021). Softer particles, however, migrate more quickly, and they migrate closer to the channel centre; therefore, the initial lateral position of softer particles is less important. Furthermore, the faster lateral migration of softer particles provides more opportunity for a soft lagging particle to remain within interaction range of the leading particle, resulting in a captured pair rather than being scattered.

Schaaf et al. (2019) did not observe the swap & capture and pass & capture types for rigid particle pairs. Figure 7 shows that both types become less common when $La$ increases. Thus, the results suggest that a critical Laplace number exists where these trajectory types disappear. Similarly, we found partially stable particle pairs only for ${La} < 100$.

Examination of figure 7 reveals that a small variation in the initial position can have a large effect on the trajectory type of the two particles. To investigate this effect further, two sets of configurations are selected where the Laplace number, the initial axial distance between the particles and the initial lateral position of the leading particle remain constant. The initial lateral position of the lagging particle is varied between $-0.5h$ and $+0.5h$. The selection of initial conditions is highlighted by the dashed areas in figure 7.

The time evolution of the axial distance between both particles is shown in figure 8. For a given value of $La$ and changing initial position of the lagging particle, various trajectory types can be seen, with stable pairs forming in some cases. Softer particle pairs (${La} = 10$) result in more trajectory types where captured particle pairs are created, but some similarities exist between the cases for ${La} = 10$ and $100$. For ${La} = 10$, the trajectory type transitions from capture when the lagging particle is located at the extreme of the opposite side of the channel centreline ($z_0^{lag} = -0.5 h$) to swap & capture as the lagging particle is released closer to the centreline. Similarly, for ${La} = 100$, the capture trajectory type also occurs when the lagging particle is at the extreme of the opposite side of the channel centreline. Once the initial position of the lagging particle is closer to the channel centreline, the early interaction for the stiffer particles (${La} = 100$) is similar to that of the softer particles (${La} = 10$) with some degree of swapping of lateral positions between the particles. However, differences begin to occur after this early interaction: the lagging particle moves away from the leading particle for ${La} = 100$, resulting in a swap & scatter type. Upon $z_0^{lag}$ becoming positive (i.e. both particles are initially on the same side of the channel), another trajectory type transition occurs: for ${La} = 10$, we find a transition to pass & capture, while the pair at ${La} = 100$ transitions to pass & scatter, marking a further difference between the cases for both Laplace numbers. The point of transition to a passing trajectory is also $La$-dependent; $z_0^{lag} = 0$ for ${La} = 100$ and $z_0^{lag} = 0.2 h$ for ${La} = 10$. Finally, for the largest studied values of $z_0^{lag}$, the trajectory types transition to scatter, irrespective of the Laplace number. All transitions described here are also visible in figure 7.

Figure 8. Time evolution of the axial distance $\delta x$ for various initial lateral positions of the lagging particle (shown in different colours) for (a) ${La} = 10$ and (b) ${La} = 100$. The line style denotes the resulting trajectory type. In all cases, $z_0^{lead} = 0.4h$ and $\delta x_0 = 3a$.

On a finer sweep of the initial position range of the lagging particle where $z_0^{lag}$ was varied by increments of $0.033 h$ (data not shown), no additional transitions were found, confirming that not all trajectory types exist for a given Laplace number.

For the same system properties (Reynolds number, Laplace number, channel aspect ratio, particle-channel confinement), different initial particle configurations can lead to entirely different outcomes: captured or scattered pairs. Captured pairs are observed within a wider range of initial positions when the particles are softer. Furthermore, depending on the initial configuration, the axial distance between particles in (partially) stable pairs can be different. Figure 8(a) shows the partially stable axial distance can vary between around $2.5 a$ and $3.5a$, a variation of nearly $50\,\%$. For the swap & capture and pass & capture types, the axial distances are different and also vary with the initial particle position. For the capture trajectory type at ${La} = 10$, the variation in final axial distance is smaller, however it remains dependent on the initial particle position. In contrast, the final axial distance of capture trajectory type pairs at ${La} = 100$ are independent of initial positions (figure 8b). We explore this behaviour further in § 4.3 since these are important observations given the formation of stable pairs is useful in many microfluidic applications when the particle spacing within a channel is required to be predictable and reproducible.

As we will elaborate next, the role of Laplace number and lateral equilibrium position are tightly connected with the stability of a captured particle pair.

4.3. Effect of softness on lateral and axial particle migration

After having identified different interaction types of soft particle pairs, we now turn our attention to the dynamics of particle capture. In their investigation of rigid particle pairs under the same flow conditions as this work, Schaaf et al. (2019) found that stable particle pairs form for certain initial lateral positions. Rigid particles show a damped oscillation of their relative distance before reaching their stable equilibrium configuration. We observed a similar oscillation for soft particles.

We investigated the oscillation of the relative distance between soft particles for $z_0^{lead} = 0.2 h$, $z_0^{lag} = -0.2 h$ and $\delta x_0 = 3a$ in the range ${La} = [5,125]$. All these cases have capture trajectories and lead to either stable or partially stable pairs. Figure 9(a) shows the time evolution of the lateral position of both particles as a function of $La$. In the early stages of migration, each particle undergoes a damped oscillation before reaching a $La$-dependent lateral equilibrium position. Trajectories of single particles at the same $La$ are included for comparison. Single particles follow similar migration paths to the particles within the pair of equal stiffness but without oscillation, demonstrating that the oscillations are a result of particle–particle interaction. While the softest particles ($La \leq 15$) migrate to the channel centre, stiffer particles reach equilibrium positions farther away from the channel centre as $La$ increases, converging to rigid limits extracted from Schaaf et al. (2019). This observation is in agreement with the results of Kilimnik et al. (2011).

Figure 9. Lateral motion of soft particles for $\delta x_0 = 3a$ and various Laplace numbers. (a) Time evolution of single, leading and lagging particles. Rigid limits extracted from Schaaf et al. (2019) for lateral equilibrium positions of rigid particles under the same flow conditions. (b) Zoomed area of (a) in the region of transition between off-centre and centreline lateral equilibrium positions. (a) Time evolution of lateral particle position. (b) Lateral equilibrium positions.

Figure 9(b) shows a zoomed area of figure 9(a), highlighting the transition between off-centre and centreline lateral equilibrium positions, also in comparison with the single-particle lateral trajectories. The point of transition between off-centre and centreline equilibrium positions for a single particle is termed the critical Laplace number and was found to be ${La}_{cr} \approx 18$ for a single particle (the full migration time at the critical Laplace number is large and not fully shown in figure 9b). However, figure 9(b) shows a particle pair with ${La} = 18$ migrating to an off-centre equilibrium position, as does a slightly softer pair with ${La} = 17$. These observations demonstrate that the formation of pairs slightly alters the inertial migration characteristics of soft particles.

The damped oscillation is also seen in the axial distance between both particles (figure 10). While particle pairs with ${La} > 15$ converge to a stable axial equilibrium distance, softer particles fail to do so. However, their axial distance remains small during the runtime of the simulations, indicating that these pairs are partially stable, rather than scattered. We find that pairs that form stable axial equilibrium distances have off-centre lateral equilibrium positions, while partially stable pairs have centreline lateral equilibrium positions in all cases. This observation can be attributed to the absence of shear-gradient-induced lift at the channel centreline. As a result, other forces have larger contributions to the overall dynamics of the system while the influence of numerical artefacts cannot be fully ruled out.

Figure 10. Time evolution of axial distance for various Laplace numbers. The initial conditions are the same as in figure 9.

We observed that, for the investigated parameter range, all pairs at the centreline are partially stable while all stable pairs are off-centre. We hypothesise that the determining factor for stability is lateral equilibrium position, rather than particle softness (Laplace number). For a given confinement and Reynolds number, the Laplace number determines the lateral equilibrium position. In particular, at $\chi = 0.4$ and ${Re} = 10$, we found ${La}_{cr} \approx 15$. To test the hypothesis, we investigated the behaviour of a particle pair at $\chi = 0.2$ and ${La} = 15$ for which we expect a single particle and particle pairs to assume an off-centre equilibrium position. Figure 11 shows the trajectories of particle pairs for two initial axial distances: $\delta x_0 = 3a$ to have the same $\delta x_0/a$ ratio and $\delta x_0 = 6a$ to have the same $\delta x_0/h$ ratio as the case with $\chi = 0.4$. The lateral migration of a single particle with the same softness, confinement and initial lateral position is included for comparison. Both sets of initial conditions lead to the pair migrating to an off-centre equilibrium position (figure 11a) and to stable pairs forming (figure 11b). Hence, the findings for $\chi = 0.2$ provide evidence that stable pairs require off-centre equilibrium positions. Applications relying on finely tuned axial distances between particles, therefore, might benefit from arrangements with sufficiently small $\chi$ and large $La$ for which particles do not migrate to the centreline.

Figure 11. Time evolution of a single particle and particle pair motion with particles of ${La} = 15$ and $\chi = 0.2$. Initial lateral positions of the particles are identical to figures 9 and 10 in non-dimensional units while the initial axial distance is equal to figures 9 and 10 in terms of channel height ($\delta x_0 = 6a$) and in terms of particle radius ($\delta x_0 = 3a$). (a) Time evolution of lateral particle positions. (b) Time evolution of axial particle distance.

The dependence of axial distance on initial configuration when particles are located at the centreline may explain some observations previously reported where the axial distance between particles is distributed within a certain range, with most of the measured distances clustered around one or several preferred distances (Lee et al. 2010; Kahkeshani et al. 2016). In realistic problems, a channel contains more than two particles, and each scatter interaction means that the leading particle may catch up with another particle that is farther downstream in the channel. This way, after some time, stable trains might form, even if individual pairs are not stable. Future studies should therefore investigate how the dynamics of single pairs is related to the dynamics of longer particle trains.

Figure 12(a) shows the absolute value of the lateral equilibrium position of the leading and lagging particles as a function of inverse $La$; the equilibrium position of a single particle under the same conditions is included for comparison. Also included for comparison is the equilibrium position of a rigid pair as reported by Schaaf et al. (2019) with which we obtain good agreement. There are two important findings beside the fact that softer particles end up closer to or at the channel centre. First, the absolute values of the lateral positions of both particles in a pair are virtually indistinguishable from each other. However, the lagging particle is slightly farther from the channel centre than the leading particle. This finding was also observed by Schaaf et al. (2019) for rigid particles. Second, both particles in a pair are consistently farther away from the channel centre than a single particle under the same conditions. This difference becomes smaller for increasing $La$, and our findings suggest that in the rigid limit ($1/{La} \to 0$), lateral equilibrium positions are the same for a single particle and particles in a pair. Indeed, Schaaf et al. (2019) found rigid particles having the same lateral equilibrium position when in a pair and in isolation. Interestingly, Gupta et al. (2018) found lateral equilibrium positions of a rigid particle train to be slightly closer to the channel wall than for a single particle. It is currently unclear how a train of soft particles under significant confinement would behave in relation to a single particle. Furthermore, we observed that the lateral equilibrium positions are independent of initial positions and the type of the particle interaction (data not shown).

Figure 12. (a) Lateral equilibrium positions for leading particle (blue), lagging particle (green) and a single particle under the same conditions (red). Note that leading and lagging particle have essentially the same lateral equilibrium positions and green markers are not always visible. Crosses indicate the difference between the lateral positions of the leading particle in a pair and the single particle. (b) Axial equilibrium distance. (c) Focusing times for the leading particle in a pair in the lateral direction and for the axial distance. The focusing time is defined as the time until the last occurrence of the position/distance being outside its equilibrium value $\pm$ the specified tolerance. Three different tolerances are included to highlight the general trends. (a) Lateral equilibrium positions. (b) Axial equilibrium distances. (c) Focusing time of particles with varying Laplace number for lateral equilibrium position and axial distance.

For the stable pairs, where the axial distance $\delta x$ converges to a constant value after some time, figure 12(b) shows the axial equilibrium distance as a function of $1/{La}$. We find that particles tend to form closer pairs when they are softer, and the lowest stable axial distance we observed is around $2.4a$, just above one particle diameter. Note that these particle pairs are also off-centre and on different sides of the channel centreline, so the actual distance between the particles is larger. We include the axial distance reported by Schaaf et al. (2019) for rigid pairs for comparison. The explanation of the gap between the rigid limit and the most rigid particle we have simulated is currently unclear. The axial distance could be strongly sensitive to mild particle deformation ($1/{La} < 0.01$), and future investigation of this parameter range is recommended.

The focusing time and distance, the time or distance until the last occurrence of a particle outside its equilibrium position within a given tolerance, is an important parameter in many inertial microfluidic applications. Figure 12(c) shows the focusing time of the leading particle in the lateral direction and the axial distance between the particles in a pair. Three tolerances have been selected for each direction as representative low, medium and high values. The magnitude of these tolerances are arbitrary and intended to demonstrate the general trend of lateral and axial focusing times. Note that the absolute tolerances in the lateral and axial directions are equal to reflect tolerances in practical applications. We find that the focusing time decreases slightly in the lateral direction for softer particles in the range ${La} = 50$–$125$. However, for soft particles below ${La} = 50$, the focusing time tends to increase, surpassing the focusing time of the most rigid particles investigated when ${La} = 20$. This general trend is repeated for the focusing time in the axial direction, with the softest particles having the longest focusing time. Given that the softest particles have lateral equilibrium positions closest to the channel centreline, where the axial velocity is larger, this also corresponds to the softest particles having the largest focusing distance.

Figure 12(c) shows that the lateral focusing occurs before the axial focusing. We will return to this observation in § 4.4 where we hypothesise that the particles within a pair must reach their lateral equilibrium position before being able to find their axial equilibrium distance.

Summarising the results so far, we have observed a number of ‘rules’ that all simulated particle pairs obey for the chosen value of Reynolds number (${Re} = 10$) and particle-channel confinement ($a/h = 0.4$) and the range of Laplace numbers and initial positions investigated.

1. (i) For all partially stable pairs, both particles are located at the channel centre.

2. (ii) For all stable pairs, both particles are off-centre and on different sides of the centreline and with essentially identical $\|z_{eq}\|$, i.e. $z_{eq}^{lag} \approx -z_{eq}^{lead}$.

3. (iii) A critical Laplace number exists for a given confinement ratio where particles with ${La} \leq {La}_{cr}$ always migrate to the channel centre and never form stable pairs.

These observations suggest that hydrodynamic particle interactions are fundamentally different for on-centre and off-centre pairs. Our findings, therefore, have important implications for microfluidic applications that rely on well-tuned axial distances between particles.

Finally, we observed that focusing of the lateral position occurs before stabilising the axial distance. We propose a high-level breakdown of the formation of stable particle pairs, shown schematically in figure 13. We distinguish between zeroth-, first- and second-order effects. Zeroth-order effects are those that occur without particles present, i.e. the flow field development. First-order effects denote the behaviour of single particles that cannot be explained by the unperturbed flow field alone, in particular the single-particle lateral migration. Second-order effects are caused by the interaction of two particles, including the oscillation of the lateral position and the stabilisation of the axial distance between particles. Our results reveal that second-order effects strongly depend on the first-order effects while the presence of a second particle has only minor consequences for the lateral motion of the other particle. The notion of second-order effects depending on first-order effects have been observed experimentally; Gao et al. (2017) found that particle trains only begin to form once particles reach their lateral equilibrium positions. As a result, microfluidic designers must consider the key parameters of the application carefully. Devices that are used for cytometry rely on axial particle ordering (second-order effects) and will require a different quality of consideration than devices for particle separation (mostly relying on first-order effects).

Figure 13. Contributions to the overall formation of stable particle pairs.

4.4. Effect of initial position on lateral and axial particle migration

To better understand how stable pairs form, we analyse the oscillations of captured pairs when initial positions are varied. We first investigate the effect of initial axial distance $\delta x_0$. Figure 14 shows the lateral positions of the leading and lagging particle and axial distance in time for initial distances in the range $[3a,11a]$. The Laplace number is set to ${La} = 36$ so that the lateral equilibrium position is equal to the initial position and the time development is mostly caused by axial rearrangement. Under these conditions, lateral oscillations are small (figure 14a); they arise from the flow field perturbations caused by the axially approaching particles (figure 14b).

Figure 14. Time evolution of (a) lateral particle positions and (b) axial distance between particles at ${La} = 36$ and $\chi = 0.4$. Initial lateral positions are set to $z_{eq}$. The initial axial distance is increased until a captured pair does not form before $t/t_{ad} \geq 10\,000$.

We identify the time $t_{ax}$ it takes for the axial distance to reach its first minimum as illustrated by vertical lines in figure 14(b). The resulting times $t_{ax}$ are shown in figure 15(a). It can be seen that the axial attraction time $t_{ax}$ has an exponential relationship with initial axial separation, with a behaviour $t_{ax} / t_{ad} \propto \exp (0.85 \delta x_0 / a)$. Thus, even a small increase in initial particle separation might lead to pairs not forming before the particles have reached the device outlet.

Figure 15. (a) Axial attraction time $t_{ax}$ at which the axial distance reaches its first minimum as denoted by vertical lines in figure 14(b). The solid line is an exponential $\propto \exp (0.85 \delta x_0/a)$. (b) Time evolution of lateral position of the lagging particle for different initial axial distances. (c) Time evolution of axial distance for different initial axial distances. (a) Axial attraction time as a function of initial axial distance. (b) Lateral position of the lagging particle vs shifted time $t' = t - t_{ax}$. (c) Axial distance vs shifted time $t' = t - t_{ax}$.

Figure 15(b,c) shows the lateral position of the lagging particle and the axial distance shifted in time such that the first minima of $\delta x(t)$ coincide ($t \to t' = t - t_{ax}$). We find that the amplitudes of both the lateral and axial oscillations increase with $\delta x_0$, converging to a constant magnitude when $\delta x_0 \approx 7a$. These observations suggest two phases exist in the axial migration of particle pairs. The first phase at $\delta x \geq 7a$ is slow, while the second phase at $\delta x < 7a$ is faster and includes stronger particle interactions leading to oscillations in both lateral position and axial distance. Importantly, the data shows that the second phase is not strongly affected by the first phase. Only for cases where particles are initially closer than $7a$, the system is initialised directly in the second phase, and the outcomes strongly depend on the initial distance $\delta x_0$.

Next we investigate the effect of initial lateral position $z_0$. The parameters are ${La} = 90$ and $\chi = 0.4$ for which the lateral equilibrium position is $z_{eq} \approx \pm 0.3 h$ and the axial equilibrium distance is $\delta x_{eq} \approx 3 a$ for all cases investigated. Particles are released at $\delta x_0 = 3 a$ and on different sides of the centreline, i.e. $z_0^{lag} = -z_0^{lead}$. Figure 16(a,b) shows the time evolution of the lateral position of both particles and their axial distance for a range of initial lateral positions, $[0.1 h, 0.45 h]$. Oscillations are visible in both the lateral position and the axial distance. The oscillation amplitude tends to be larger when the initial position $z_0$ is farther away from $z_{eq}$. As expected, for the special case $z_0 \approx z_{eq}$, there are virtually no oscillations since particles are initialised close to their equilibrium configuration. For the majority of cases, the initial lateral motion of both particles is towards the equilibrium position. However, lagging particles initially located between their equilibrium lateral position and the channel centreline behave differently: these particles first move away from their equilibrium position (towards the centreline) before changing direction and eventually converging at the equilibrium. Given that single particles in our simulations always move towards their equilibrium lateral position, this observation can be attributed to hydrodynamic particle–particle interactions.

Figure 16. Time evolution of (a) lateral particle position and (b) axial distance at ${La} = 90$, $\chi = 0.4$ and $\delta x_0 = 3a \approx \delta x_{eq}$. Initial lateral positions obey $z_0^{lag} = -z_0^{lead}$ and are varied in the range $[0.1h, 0.45h]$. (c) Lateral centre of mass position of the pair vs axial distance; the inset shows the zoomed area close to initial and equilibrium positions. (d) Time evolution of the half-lateral distance between particles compared with trajectories of single particles. Note that a reduced selection of initial positions are included in (c) to improve readability. The line colours in all panels correspond to the legend in (a). (a) Lateral positions vs time. (b) Axial distance vs time. (c) Lateral centre of mass position vs axial distance. (d) Half-lateral distance vs time.

The direction of the initial oscillation of the axial particle distance (figure 16b) differs either side of the lateral equilibrium position with growing amplitude as the difference between initial and equilibrium lateral position increases: when particles in a pair are initially closer to the centreline than their lateral equilibrium positions, particles first decrease their axial distance. Particles initially farther away from the centreline first increase their axial distance. The inset in figure 16(b) shows the two extreme cases with initial lateral positions of $0.1 h$ and $0.45 h$. Despite the amplitude of the first oscillation of both pairs being similar, the oscillation of the pair farthest away from the channel centreline damps more quickly than the pair closest to the centreline. We will explore the damping of the oscillations in more detail in § 4.5.

We now consider the centre-of-mass behaviour of the particle pair. Figure 16(c) shows the lateral centre of mass position vs the axial particle distance. Since particles are released with the same axial distance and at the same distance either side of the channel centre, the centre of mass always starts at the same point on the centreline (circle in figure 16c). The inset shows a zoomed area close to the initial and equilibrium positions. At early times, the centre of mass leaves the centreline, begins to oscillate and eventually spirals towards the equilibrium back on the centreline. The initial direction of the motion of the centre of mass depends on the initial lateral position of the particles. However, as the pairs approach the equilibrium positions at later times, all trajectories converge to a similar counterclockwise path. This convergence implies that the time history of the pair is eventually forgotten, which is in line with observations in figure 15 suggesting that the behaviour during the second phase is largely independent of that during the first phase of migration.

Finally, figure 16(d) shows the time evolution of the half-lateral distance between particles in a pair, compared with the lateral position of a single particle with the same initial position. Generally, both types of curves follow the same trend, i.e. to first order, particles in a pair and single particles behave similarly. Interestingly, particles in a pair reach their lateral equilibrium distance faster than the single particle converges to its lateral equilibrium position when particles are initially close to the centreline. This difference could be caused by the initial proximity of both particles which repel each other hydrodynamically, hence accelerating the initial lateral migration away from the centreline. Particles in a pair that are initially farther away from each other experience a weaker hydrodynamic repulsion and behave nearly like single particles. Another important observation is that the lateral oscillations seen in figure 16(a) are nearly completely absent in figure 16(d). On long time scales – of the order of hundreds of advection times – particles can change their lateral distance and migrate towards their equilibrium; but on the shorter time scale of the lateral oscillations – tens of advection times – particles do not oscillate relatively to each other. We hypothesise that a fast oscillation of the lateral particle distance is suppressed by the liquid between the particles: a fast lateral oscillation would require liquid being repeatedly squeezed in and out of the gap between the particles, which is energetically unfavourable. Therefore, the particles’ lateral positions are hydrodynamically coupled on a short time scale, which also explains why in some cases one particle in a pair initially moves away from its equilibrium position (figure 16a).

Taking all observations from figure 16 together, we can now explain the counterclockwise sense of the oscillation seen in figure 16(c). During periods when the lagging particle is closer to the centreline than the leading particle, the lateral centre of mass position is positive and the lagging particle is able to catch up with the leading particle ($\delta x$ decreases) since the lagging particle is exposed to a faster portion of the flow field. Once particles approach axially, the gap between both particles is decreasing, and liquid needs to be pushed out of the gap, causing a repulsion of the particles. This repulsion force affects both particles differently. The particle currently closer to the centreline experiences a smaller shear-gradient force and can move more easily away from the centreline than the particle currently in a higher shear-gradient region. As a consequence and supported by particle inertia, the lateral centre-of-mass location crosses the centreline and turns negative, and the leading particle becomes faster, therefore increasing the axial distance $\delta x$. Upon increasing the axial distance, the gap grows and liquid needs to move into the gap, therefore causing particle attraction. The process then continues with swapped roles of both particles and repeats itself. While inertia drives this oscillation, viscosity causes the damping and eventual convergence to the equilibrium state.

Having identified the importance of oscillations in the formation of a pair, we now focus our investigation on the dynamics of these oscillations and their effect on focusing time.

4.5. Capture oscillation dynamics and focusing time

Previous studies have identified the existence of a damped oscillation of the relative distance between captured particles (Schaaf et al. 2019; Hu et al. 2020; Udono 2020). To understand better how stable pairs form, we analyse the oscillations of the axial distance between particles when different parameters are varied: Laplace number, initial lateral position $z_0$ and initial axial distance $\delta x_0$. Note that the trajectories analysed in this section have already been presented in § 4.3 for variable Laplace number and § 4.4 for variable initial positions. We evaluated the damped frequencies by measuring the time periods between oscillation peaks, similar to the earlier analysis of rigid particle pairs (Schaaf et al. 2019). Figure 17(a) shows that the oscillation frequency increases with $La$ until approximately ${La} > 60$ beyond which the frequency becomes approximately constant. The oscillation frequency of the softest particle pair studied is approximately $50\,\%$ smaller than that of the most rigid particle pair studied. The inset in figure 17(a) reveals the results of a fast Fourier transform of the time data, confirming the findings in the main graph.

Figure 17. Oscillation frequency $\varOmega$ as a function of (a) Laplace number, (b) initial lateral position (${La} = 90$) and (c) initial axial distance (${La} = 36$). The inset in (a) shows results from a fast Fourier transform of the time data for a selection of Laplace numbers (colours correspond to those in figure 10a). (a) Frequency as function of softness. (b) Frequency as a function of initial lateral position (${La} = 90$). (c) Frequency as a function of initial axial distance (${La} = 36$).

We investigated the influence of initial configuration on the resulting oscillation frequency during pair capture. Figure 17(b,c), respectively, show the obtained frequencies as a function of initial lateral position for ${La} = 90$ and initial axial distance for ${La} = 36$. Figure 17(b) corresponds to the cases in figure 16, and figure 17(c) corresponds to the cases in figure 14. We see a slight increase in frequency with lateral position and a slight decrease in frequency with initial axial distance. However, the variation in both relationships is over an order of magnitude smaller than that of the variation with Laplace number. We conclude that particle softness is the key determining factor for the oscillation behaviour, and initial positions have a negligible effect on the late-stage capturing dynamics where oscillations are strong. These findings support our hypothesis that the particle dynamics during the second phase of migration are largely independent of the initial details and the first phase of migration.

Next we analyse the damping of the oscillations of the axial distance. To obtain the damping coefficient $\gamma$ for each case, we fitted decaying exponentials of the form $A \exp (-\gamma (t - t_0)) + \delta x_{eq}$ to the maxima and minima of the time data $\delta x(t)$. Since the initial axial distance $\delta x_0$ generally has a mismatch with the axial equilibrium distance $\delta x_{eq}$, the oscillations of $\delta x(t)$ show pronounced transients at early times. To improve results, we ignored the first two periods of each curve in the fitting process. The inset of figure 18(a) illustrates our procedure. Figure 18 shows the results for the damping ratio $\gamma$; panels (a–c) correspond to those in figure 17. Figure 18(a) reveals that the damping ratio increases with $La$. Unlike the frequency $\varOmega$ in figure 17(a), $\gamma$ keeps increasing beyond ${La} = 60$ and does not seem to converge to a ‘rigid limit’ in the range of Laplace numbers investigated. However, figure 18(b,c) shows that the damping ratio – as the frequency in figure 17(b,c) – is essentially independent of initial lateral position and initial axial distance. The behaviour of the damping coefficient provides further support to the importance of particle softness to the oscillation dynamics, while details of the initial conditions are less relevant.

Figure 18. Damping ratio $\gamma$ as a function of (a) Laplace number, (b) initial lateral position (${La} = 90$) and (c) initial axial distance (${La} = 36$). The inset in (a) depicts the process of obtaining the damping coefficient $\gamma$ (see main text for details). (a) Damping as a function of softness. (b) Damping as a function of initial lateral position (${La} = 90$). (c) Damping as a function of initial axial distance (${La} = 36$).

A key parameter in the formation of particle pairs is the focusing time for both the lateral position and the axial distance. In § 4.3 we presented the dependence of both focusing times on $La$. We showed that particles are focused laterally before the axial distance equilibrates. Figure 19 reveals the dependence of focusing time on particle softness, initial particle position and initial axial distance with moderate tolerance ($\delta x_{eq} \pm 0.05a$ and $z_{eq} \pm 0.02h$). For clarity, figure 19(a) contains the relevant subset of data from figure 12(c). As expected, figure 19(b,c) show that the initial conditions have a strong effect on the focusing time. While focusing times are minimal when particles are initialised close to their equilibrium lateral position ($z_{eq} \approx 0.3h$ in this case), focusing times increase with $\lvert z_0 - z_{eq} \rvert$ (figure 19b). Interestingly, focusing times become particularly large when particles are initially close to the centreline; probably because shear gradients are smaller in this region. This trend follows the change in oscillation amplitudes observed in figure 16. Likewise, focusing times are smallest when the initial axial distance matches the equilibrium distance ($\delta x_{eq} \approx 3a$ in this case), while the focusing time increases strongly with $\delta x_0$, in particular for $\delta x_0 > 9a$ (figure 19c). At lower values of $\delta x_0$, the increase in focusing time is much more gradual, perhaps offsetting the decrease in axial distance with the changing oscillation dynamics that are constant at larger initial axial distance (figure 15c).

Figure 19. Focusing time for lateral position of the leading particle and the axial distance as a function of (a) Laplace number, (b) initial lateral position (${La} = 90$) and (c) initial axial distance (${La} = 36$). (a) Focusing time as function of softness. (b) Focusing time as a function of initial lateral position. (c) Focusing time as a function of initial axial distance.

In our capsule model, we have chosen to constrain the bending modulus of the membrane by introducing a dependency on the elastic shear modulus so that it is no longer a free parameter. The relationship chosen between the strain and bending moduli has been derived from that of a red blood cell and, therefore, has physical significance. We suggest that in future work, the effect of bending rigidity should be investigated further.

Our analysis shows that both the initial positions and particle softness play an important role in the dynamics of stable pair formation. Particle softness dominates the characteristics of trajectory oscillations while initial positions play a larger role in determining focusing time. We conclude that designers of inertial microfluidic devices must also consider the particle positions prior to interaction, with device performance optimised through consideration of particle softness.

5. Conclusions

The formation of stable particle pairs is the foundation for many applications within the field of inertial particle microfluidics, including particle focusing, cytometry and separation. Identifying the conditions under which stable pairs form is crucial to improving the design process and performance of some inertial microfluidic devices. Despite most particles processed in microfluidic devices being soft biological cells, the majority of numerical investigations to date have concentrated on the modelling of rigid particles. It is known that particle deformability affects the migration behaviour under the influence of inertia. Here, we have investigated the effect of particle softness on the formation and stability of particle pairs in straight channels under moderate inertial flow conditions.

We used an in-house immersed-boundary-lattice-Boltzmann-finite-element solver to simulate single and pairs of soft capsules in channel flow under the influence of fluid inertia at channel Reynolds number $10$ and particle-to-channel confinement of $0.4$. The code has been benchmarked against previously published numerical results involving a pair of capsules interacting in inertial shear flow and a single soft particle migrating in a channel.

There are the following several important results.

1. (i) We first investigated the general effect of particle softness, characterised by the Laplace number, and initial particle position for a pair of equally soft particles. We observed two new trajectory types termed swap & capture and pass & capture that appear in some cases and have not been found for rigid particle pairs in the same geometry in earlier studies. Particle softness was found to increase the likelihood of a captured particle pair forming.

2. (ii) We observed that particle pairs that migrate to the channel centre are stable only in the lateral direction, but not in the axial direction (partially stable pairs). In contrast, particle pairs that migrate to an off-centre lateral equilibrium position form pairs that are stable in both the lateral and axial directions. We found that particle softness seems to affect the pair stability only through the resulting lateral equilibrium position. Generally, the lateral equilibrium positions of soft particles in a pair are nearly the same as those of single particles under the same conditions. Importantly, our simulations show that the stabilisation of the axial particle distance occurs only after particles have migrated laterally. Our findings suggest that the observations can be classified as zeroth-, first- and second-order effects which refer to effects caused by the background flow field, the behaviour of a single particle and the interaction of two particles, respectively.

3. (iii) The formation of a stable pair consists of two phases: an early and a late phase. During the early phase, particles migrate laterally and approach each other until the axial distance becomes small. The later phase involves a spiralling motion leading to a converged state. We found the late phase to be largely independent of the first phase, as long as the initial conditions lead to the eventual formation of a stable pair. During the spiralling motion towards the converged state, particles are tightly coupled through hydrodynamic interactions. The spiralling motion is driven by an interplay of flow-induced lift, particle inertia and viscous dissipation. Importantly, the pair formation time grows exponentially with the initial axial distance of the particles.

4. (iv) The particle oscillations during the late phase of pair formation are characterised by their frequency and damping coefficient. Both quantities increase with Laplace number until a rigid particle limit is reached, but they are largely independent of initial conditions. Finally, we show that the pair formation time is determined by the initial conditions while particle softness has only a mild influence.

Through our investigation, we have identified new physical effects in the formation of stable particle pairs that are softness-dependent. These effects should be considered in the design of inertial microfluidic devices, given that soft particles are used extensively in real-world applications. Since axial stabilisation comes after lateral migration, our observations put lower limits on the required time and distance necessary to generate stable pairs. Furthermore, devices relying on appropriate axial spacing require more consideration than devices that focus or separate particles by lateral position. Our findings could have significant ramifications for applications where a constant and reliable axial distance between particles is required, such as flow cytometry.

Future work could include a detailed analysis of the role of the flow field in the formation of pairs and trains of particles and the development of reduced-order models to predict pair formation without the need for resolved simulations. Since many inertial microfluidic applications involve heterogeneous particle mixtures, an investigation of heterogeneous pairs and trains over a larger range of Laplace numbers would provide further insight.