Particle Shape Influences Settling and Sorting Behavior in Microfluidic Domains

We present a new numerical model to simulate settling trajectories of discretized individual or a mixture of particles of different geometrical shapes in a quiescent fluid and their flow trajectories in a flowing fluid. Simulations unveiled diverse particle settling trajectories as a function of their geometrical shape and density. The effects of the surface concavity of a boomerang particle and aspect ratio of a rectangular particle on the periodicity and amplitude of oscillations in their settling trajectories were numerically captured. Use of surrogate circular particles for settling or flowing of a mixture of non-circular particles were shown to miscalculate particle velocities by a factor of 0.9–2.2 and inaccurately determine the particles’ trajectories. In a microfluidic chamber with particles of different shapes and sizes, simulations showed that steady vortices do not necessarily always control particle entrapments, nor do larger particles get selectively and consistently entrapped in steady vortices. Strikingly, a change in the shape of large particles from circular to elliptical resulted in stronger entrapments of smaller circular particles, but enhanced outflows of larger particles, which could be an alternative microfluidics-based method for sorting and separation of particles of different sizes and shapes.

Flow and transport of engineered particles of different geometrical shapes are encountered in diverse biomedical applications. In targeted drug deliveries, the shape of engineered drug cargos has shown to have intriguing effects on their transport in blood vessels, adhesion onto channel walls, and targeting ability toward malignant cells 1 . For example, ellipsoidal microparticles displayed longer blood circulation times than spherical particles due to less efficient phagocytosis by macrophages in the reticuloendothelial system 2 . Hexagonal nanoparticles more effectively mitigated phagocytoses and remained in blood circulation longer than spherical particles 3 . Unlike spherical particles, boomerang-shaped particles displayed a preferred direction of Brownian motion 4 , which could have implications in design of new microscopic particles to deliver drugs or self-assemble into complex materials. A theranostic plasmonic shell-magnetic core star-shaped nanomaterial was used for targeted isolation and detection of rare tumor cells from a blood sample 5 . As for the adhesion kinetics of such engineered particles on channel walls, nanorod particles were numerically shown to adhere to channel walls easier than spherical particles due, in part, to larger surface area contacts with the channel walls as they tumble near the walls 6 .
In applications relevant to the design of biomedical devices, microfluidic devices with different geometric designs have been proposed to isolate circulating tumor cells (CTC) from healthy cells in blood samples through, for example, vortex-aided particle separation 7,8 , which could be useful for early cancer diagnosis and monitoring metastatic progression or the efficiency of cancer treatments 9 . Although the performance of the microfluidic devices in the segregation of CTC has been commonly tested with surrogate spherical particles, tumor cells often exhibit patient-specific arbitrary shape profiles, which do not conform to the spherical particle representation for tumor cells 10,11 .
The effect of non-spherical particle shapes on particle trajectories has been recently addressed in numerical simulations. Settling dynamics and patterns of thin disks 12,13 in an infinitely long viscous fluid domain and settling behaviors of individual spherical, cubical, or tetrahedral particles in an infinitely long fluidic domain with periodic lateral boundaries 14 were numerically investigated. However, to the best of our knowledge, numerical simulations of settling of a mixture of different-shaped particles (DSP), involving angular-and curved-shaped particles, in a bounded domain is unprecedented. Similarly, numerical simulations of flow trajectories of a mixture of DSP is very limited or perhaps non-existent in the literature.
The extension of the lattice Boltzmann (LB) method for simulating flow of suspended bodies is a fast-growing area of LB research 15 , following the pioneering work of of Ladd 16,17 . Considering broad uses of DSP in biomedical applications and the abundant experimental evidence for their shape-dependent distinct flow and transport behaviors, we extended the LB model (LBM) presented originally by Nguyen and Ladd 18 to simulate the settling and flow of DSP, including discretized angular-shaped particles (DAsP), involving star, boomerang, hexagonal, triangular, rectangular, and discretized curved-shaped particle (DCsP), involving circular and elliptical particles, consistent with the aforementioned shapes of engineered particles used in biomedical applications. The DSP-LBM is suitable for simulating settling and flow trajectories of any arbitrary-shaped particles, such as tumor cells.
The primary purpose of this paper is to introduce the DSP-LBM and demonstrate its performance in simulating the settling or flow of individual or a mixture of DSP under various combinations of properties associated with the particles, flow regimes, and the microfluidic domain geometry. Using the DSP-LBM and a single chamber of the microfluidic device geometry in ref. 7 we numerically investigated the validity of recent findings and implications in microfluidic research. These findings and implications involve: (i) when a large number of particles are released into a fluid in a microfluidic device, larger particles get selectively trapped by vortices, whereas smaller particles avoid entrapments; (ii) steady vortex structures can be used to quantify vortex-controlled, size-based separation of particles; and (iii) non-circular particles may be represented by circular particles in vortex-aided particle segregation via microfluidic devices with different geometric peculiarities.

Methods
In the LB method [19][20][21][22] , the mesodynamics of the Newtonian fluid flow can be described by a single relaxation time via the Bhatnagar-Gross-Krook (BKG) equation 23 is the set of population densities of discrete velocities e i at position r and discrete time t with a time increment of Δt, τ is the relaxation parameter, and f i eq is the local equilibrium 24 , The local fluid density, ρ, and velocity, u, at the lattice node are given by ρ = ∑ f i i and ρ τρ , where g is the strength of an external force 25 . A D2Q9 (two-dimensional nine velocity vector) lattice 21 was adopted in numerical simulations. Through the Chapman-Enskog approach, the LB method for a single-phase flow recovers the Navier-Stokes equation in the limit of small Knudsen number for weakly compressible fluids, in which ∇ ⋅ ∼ u 0 and ∂ t u + (u ⋅ ∇)u = −(∇P/ρ) + ν∇ 2 u + g with the fluid kinematic viscosity, ν τ = Δ − . c t( 05) s 2 . Pressure, P, is computed via the ideal gas relation, ρ = P c s 2 . The extension of the LBM to the DSP-LBM involves (i) geometric description of DSP to locate the vertices for DAsP or boundary nodes for DCsP, (ii) calculations of the position of intra-and extra-particle boundary nodes in the vicinity of arbitrary-shaped particle surfaces across which the particle and fluid exchange momentum, and (iii) calculations of new positions of the center of mass of a particle and its vertices based on particle-fluid hydrodynamics.
Geometric Description of 2D Different-shaped Particles. Similar to geometric construction of surfaces of a circular-cylindrical particle (hereafter, circular particle) by Ladd 18 , we used discretized particle surfaces for 2D curved (e.g., circular)-and angular (e.g., hexagonal)-shaped particles in the DSP-LBM. A schematic representation of non-circular particle geometries are shown in Fig. 1, which are subsequently used to locate vertices of DAsP and boundary nodes of DCsP. We provide geometric descriptions for the star-shaped and elliptical particles next, but geometric descriptions of the remaining particles are provided in Supplementary Information-1.
The star-shaped particle geometry is represented by five isosceles triangles connected to a pentagon at the center, as shown in Fig. 1a. The geometry is constructed by two circles; the bigger circle with a radius of R S encloses the star-shape and the smaller circle with a radius of R P that passes through the corners of the pentagon. These two circles are related via R S = ψR P , in which ψ = cos(π/5) + [sin(π/5)]/[tan(π/10)]. The surface area of the star-shaped particle, A S , is given by . The star-shaped particle has five vertices located on the outermost tip of the triangles (v S1 − v S5 ), in addition to five vertices located on the corners of the inner pentagon (v P1 − v P5 ) (Fig. 1a). The coordinates (x i , y i ) of v Si , and v Pi , where iε 1,5 , are computed by Eq. 2 and Eq. 3, respectively,  ) for DAsP or the number of boundary nodes (N Bnd ) for DCsP. The mass of the star-shaped particle per unit particle thickness is given by m p = χ(R S ) 2 ρ p . The moment of inertia, I s , for the star-shaped particle was computed by pentagon, = + + A a (1/4) 5 (5 2 5) P 2 , A T is the area of the triangle, A T = ah/2, λ = R p cos(π/5), and σ = [cos( π/5) + cos(2π/5)] 2 + [0.5 + sin(π/5) + sin(2π/5)] 2 .
Different from a star-shaped particle, the elliptical particle geometry is described by boundary nodes, N bnd , along the discretized curved surfaces, the length of its long-and short-axes (c and d), and the initial tilt angle, α (Fig. 1b). The coordinates of its boundary nodes are computed by . The mass of an elliptical particle per unit particle thickness is given by

Intra-Particle Boundary Nodes (IPBN) and Extra-Particle Boundary Nodes (EPBN). The wind-
ing number algorithm 26 was implemented to determine whether a lattice node x k = (x k , y k ) is enclosed by a polygon in Fig. 1, described by a series of boundary nodes for DCaP or vertices for DAsP along the particle surface. The algorithm computes the number of times the polygon winds around x k , which is referred to as the winding number, m(x k ). x k is not enclosed by a polygon if m(x k ) = 0. In the DSP-LBM, x k and x k + e i form a intra-particle boundary nodes (IPBN) and extra-particle boundary nodes (EPBN) pair if m(x k ) ≠ 0 and m(x k + e i ) = 0. The IPBNs and EPBNs for a discretized hexagonal particle and the momentum exchanges between the particle and the fluid at the mid-point of hydrodynamic links connecting an IPBN and an EPBN are shown in Fig. 2.
Particle-fluid Hydrodynamics. Particle-fluid hydrodynamic calculations rely on momentum exchanges between the fluid and the mobile DSP, following the approach in refs 18,27 in which the population densities near particle surfaces are modified to account for momentum-conserving particle-fluid collisions. Particle-fluid hydrodynamic forces, F r b , at the boundary nodes located halfway between the intra-particle lattice node, r v , and extra-particle lattice node, r v + e i , are computed by 16,28,29  The translational velocity, U p , and the angular velocity of the particle, Ω p , are advanced in time according to the discretized Newton's equations of motion, where m p is the particle mass, I p is the moment of inertia of the particle, and The new position of the center of mass of a particle is computed as The population densities at r v and r v + e i Δt are updated to account for particle-fluid hydrodynamics in accordance with 16 New Locations of Vertices or Boundary Nodes. The locations of vertices or boundary nodes are updated in each time step. The distance d i = (d ix , d iy ) between the i th vertex (or a boundary node) and the center of mass of a particle, x c is computed via in which ϒ i is the angle between (x i − x c ) and +x. For a hexagonal particle, for example, ϒ i = α + (i − 1)π/3 for iε 1,6 .
Model Validation. The DSP-LBM was validated with two benchmark problems. First, the settling trajectory of a circular particle in an initially quiescent fluid in a bounded domain (Fig. 3a) computed by the DSP-LBM was compared against the finite-element (FE) solutions by Feng et al. 30 at two different Reynolds numbers, Re = 8. 33 and Re = 1.03 (here, R e = 2RU s /ν, where R is the particle radius and U s is the settling (terminal) velocity of the particle). In ref. 30 the values of R and ν in FE simulations were not provided, but only Re values were reported. In the DSP-LBM simulations, the length of the bounded flow domain was set to ∼ W 30 (adopted in all settling simulations in this paper), where W is the channel width perpendicular to the main settling direction, and R = 385 μm, ν = 0.01 cm 2 , and |g| = 981 cm/s 2 . ρ p /ρ was adjusted to meet the reported Re values in ref. 30 . For Re = 8.33, DSP-LBM (with ρ p /ρ = 1.07) and FE solutions are in good agreement (Fig. 3b), although the DSP-LBM solution for Re = 6.65 (with ρ p /ρ = 1.05) matched the FE solution for Re = 8.33 better. The FE solution for Re = 1.03 was in a good agreement with the DSP-LBM solution (with ρ p /ρ = 1.01) for Re = 1.68 (Fig. 3c).
In the second validation test, the DSP-LBM simulation of the settling trajectory and angular rotations θ α = + Ω Δ t ( ) p of an elliptical particle in an initially quiescent fluid in a bounded domain (Fig. 3d) was compared against numerical solutions by Xia et al. 31 . In these simulations, c/d = 2, W/c = 4 (Fig. 3d), density ratio of ρ p /ρ = 1.1, α =°45 , ν = 0.01 cm 2 /s, c = 0.1 cm, and |g| = 981 cm/s 2 . Figure 3e,f show that settling trajectory and angular rotations of an elliptical particle computed by DSP-LBM are in good agreement with the simulation results by Xia et al. 31 .
Data availability. The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Results
Settling Trajectories of Different-Shaped Particles. The DSP-LBM was used to simulate the settling trajectories and velocities of DSP as a function of particle density. The same problem set-up in Fig. 3d was used, but the elliptical particle was replaced by particles of different shapes. The blockage ratio is defined as W/R e , in which R e is the equivalent radius of a circular particle that has the same surface area of a non-circular particle. In these simulations, R e = 3.5 × 10 −2 cm, the surface area of the particle is A p = 3.9×10 −3 cm 2 , and g = 981 cm/s 2 .
The same A p was specified for all DSP by setting R s = 15.4, and R p = 5.9 for the star particle; B =15, ζ = π/3, φ = π/6 for the boomerang particle; L = 10.1 for the hexagonal particle; a = 24.8 for the triangular particle; l = 23.1, w = 11.5 for the rectangular particle; R = R e = 9.2 for the circular particle; and c = 26, d = 13 for the elliptical particle (Fig. 1). Here, the length parameters are expressed in l.u. (1 l.u. = 3.846 × 10 −3 cm) and angles are described in radians. The initial orientation of the particles are shown in Fig. 4.
DSP-LBM simulation results in Fig. 5 unveiled three distinct shape-dependent-particle behaviors in a confined channel for ρ/ρ p = 1.05: (i) the boomerang and triangular particles exhibited an initial large displacement from the centerline toward the channel wall at y = W, followed by oscillatory trajectories about the centerline while displaying the largest cumulative angular rotations; (ii) after a large displacement toward the wall at y = 0, the elliptical and rectangular particles with the same aspect ratio [(c/b) = (l/w) = 2] drifted toward the centerline and displayed nearly zero angular rotations as they gradually oriented their principal axis normal to the gravitational field; and (iii) the hexagonal and star particles settled near the centerline similar to a circular particle, but they displayed non-zero cumulative angular rotations, unlike the circular particle.
As ρ/ρ p increased from 1.05 to 1.10 (i.e., higher inertial effect), the particles exhibited more oscillations in their settling trajectories as they gradually drifted to the mid-channel. The most striking finding was the effect of the small triangular chip (BDC in Fig. 1c) on the settling trajectory of the boomerang particle. For ρ/ρ p = 1.10, the small chip was responsible for the persistent periodicity in the boomerang particle's settling trajectory, different from slowly decaying oscillations in the triangular particle's settling trajectory. Thus, DSP-LBM simulations revealed that a small chip in the boomerang geometry is a key design criteria, controlling the amplitude and frequency of the oscillations in settling trajectories of the boomerang particle.  The other design criteria for engineered DSP may include the (linearized) surface concavity of the boomerang particles and the aspect ratio of rectangular particles. The effect of the surface concavity of the trailing edge of the boomerang particle, controlled by its inner angle (φ) on its settling trajectory, is shown in Fig. 6a for ρ/ρ p = 1.10, W/R e = 11.3, ζ = 60°, and A p = 3.9×10 −3 cm 2 . DSP-LBM simulations show that the boomerang particle displayed gradually vanishing oscillations in its settling trajectory, similar to the triangular particle, if x c was located inside the polygonal surface (for φ = 10° and 20°). The boomerang particle exhibited periodic oscillations in its settling trajectory if x c was located on the polygonal surface (for φ = 30°) or outside the polygonal surface (for φ = 40°). The oscillation frequency,ϑ, dropped from 1.27 s −1 to 1.13 s −1 as x c moved from the polygonal surface (D in Fig. 1c) to an exterior point outside the polygonal surface.
The effects of the aspect ratio of a rectangular particle on its settling trajectories are shown in Fig. 6b for ρ/ρ p = 1.10 and W/R e = 9.2. Although rectangular particles with different aspect ratios drifted toward the same  equilibrium position at the centerline at ∼ x W / 6, the rectangular particle with the largest aspect ratio exhibited the largest initial displacement from the centerline and more frequent and largest oscillations in its settling trajectory, which could be critical in multi-particle flows.
The effects of particle shape on the settling (terminal) velocity of an individual particle are provided in Supplementary Information-2. As shown in Supplementary Information-3, the overall settling trajectories of DSP in these simulations are deemed to be independent of grid resolution for all practical purposes.

Flow Trajectories of Individual Particles of Different Shapes.
In the DSP settling problems discussed above, the fluid was initially quiescent. To simulate shape-dependent flow trajectories of DSP, the particles whose initial orientations shown in Fig. 4 were released into a Poiseuille flow from a point 20% off the centerline after the steady-flow field was established. A neutrally-buoyant spherical particle in a Poiseuille flow typically exhibits the Segre-Silberberg effect 32 with an equilibrium settling position between the channel wall and centerline, in which the equilibrium position varies with Re 33,34 (here, Re = 2RU ss /ν, where U ss is the average steady fluid velocity prior to releases of the particles). For ∼ .
Re 0 1 and W/R e = 6.6, the equilibrium position of the neutrally-buoyant spherical particle was on the centerline in a tube 33,35 . Consistent with these findings, different equilibrium positions of a circular particle in a Poiseuille flow computed by DSP-LBM as a function of Re are shown in Supplementary Information-4. Among them, Re = 35.2, corresponding to the average steady fluid velocity of 4.97 cm/s prior to the particle release, was chosen and the flow trajectories of DSP were simulated (Fig. 7). At Re = 35.2, the circular particle exhibited slowly diminishing overshots about the centerline in its flow trajectory due to combined effects of inertial and wall effects. At much higher Re, however, the wall effect may be confined to near-wall layers only 36 .
When compared to the settling trajectories of DSP (Fig. 5), the flow trajectories of the DSP are more sensitive to the particle shape in a flowing fluid. DSP followed distinct flow trajectories at Re = 35.2 before they drifted to their equilibrium position at ∼ x W / 25. Only DCsPs exhibited overshots in their flow trajectories. Although the settling trajectories of the elliptical and rectangular particles with the aspect ratio of 2 were similar, their flow trajectories were different, revealing the significant effect of the (discretized) curved particle surface on particle trajectories in a shear flow. Similarly, the settling trajectories of the star and hexagonal particles were similar, unlike their trajectories in a shear flow. Uniform and repetitive oscillations in the settling trajectories of boomerang and triangular particles were replaced by non-uniform oscillations in their trajectories in a shear flow. Circular, star, hexagonal, and boomerang particles displayed the largest cumulative angular rotations at Re = 35.2 while the boomerang and triangular particles exhibited the largest cumulative angular rotations as they settled.
Settling and Flow of a Mixture of DSP. The effect of particle shapes on the settling and flow behavior of a mixture of DSP was numerically demonstrated here for the first time. Four simulations were setup, through which trajectories and velocities of seven settling or flowing DSP were compared to those of seven circular particles. All particles, regardless of their shapes, had the same surface area with R e = 385 μm. The interparticle distance at the release location was 4R e and the width and length of the domain was 40R e × 80R e . The fluidic domain was bounded in the settling simulation. A periodic boundary condition was implemented at the inlet and outlet for the flow simulation for which Re = 38. Steric interaction forces, based on two-body Lennard-Jones potentials 27 , were used to avoid unphysical overlapping of particles when they are in near contact, as described in Supplementary Information-5. Figure 8a,b show that use of multiple surrogate circular particles in place of a mixture of non-circular particles led to not only misrepresentation of settling trajectories of DSP, but also underestimation of their settling velocities by a factor of up to 2.2 (large velocity ratios near the bottom boundary can be ignored as some particles rested on the bottom while the others continued to roll, which resulted in large velocity ratios). Similarly, Fig. 8c,d show that if non-circular particle shapes are overlooked, lateral displacements in computed trajectories significantly differed and particle velocities deviated by a factor of ∼ . − .

Discussion
In the preceding sections, DSP-LBM simulations demonstrated significant effects of particle shapes on the settling or flow trajectories of an individual particle or a mixture of DSP. Using the DSP-LBM, we investigated here the validity of recent findings and implications in microfluidic analyses: (i) would steady vortex structures alone be used to quantify vortex-controlled size-based sorting of particles? (ii) would larger particles be selectively entrapped in steady vortex regions despite the cumulative effects of particle-fluid hydrodynamics on the fluid velocity in relatively dense suspensions? and (iii) would the findings from vortex-controlled size-based separation of circular particles be extensible to non-circular particles in microfluidics? To answer these questions, DSP-LBM simulations were setup using a single chamber of the microfluidic geometry in ref. 7 . After the steady-flow field was established, 10 large particles of 0.38 μm in diameter and 30 small particles of 0.19 μm in diameter were released into a microfluidic chamber from random locations at the inlet. The dimensions of the microfluidic domain and the steady flow field are shown in Fig. 9. The fluid was water with ν = 0.01 cm 2 /s and c s = 1,460 m/s, and the particles were neutrally buoyant. The average flow rate, u avg , of 52.14 m/s at the inlet in a single-chambered microfluidic chamber produced vortex structures, similar to the vortex structures in a multi-chambered microfluidic device with ∼ u 1, 700 avg m/s in Fig. 3 of ref. 7 . Although steady vortex structures were previously envisioned to trap particles in microfluidic devices 7,8 , Fig. 10 shows that vortices in a flowing fluid including mobile particles are indeed unsteady, even if the pressure differential at the inlet and outlet is held constant in time. Symmetry breaks in the flow domain with initially symmetric vortex structures, disappearance or changes in the location of vortices, and formation (birth) of new vortices as a result of cumulative effects of interparticle and particle-fluid hydrodynamics are evident from Fig. 10. Particle motion in this case is largely determined by momentum exchanges between the particles and unsteady discrete vortices, similar to the underlying reasoning of a steadily swimming fish in a water with discrete vortices 37 , for which Lagrangian coherent structures are typically used to decompose unsteady fluid flows into dynamically different regions. In brief, for the initial flow condition given in Fig. 9 as in ref. 7   involving multiple mobile particles was inherently transient, which contradicts the use of steady vortex regions 7,8 in assessing particle entrapments in microfluidics devices. Moreover, the sorting mechanism related to correlations between the lateral displacements of particles to their sizes 38 is not applicable for multi-particle simulations in a fluidic domain in Fig. 9.
Next, DSP-LBM simulations were used to investigate if the larger circular particles are selectively trapped in unsteady vortex regions (Fig. 10). Particles leaving the flow domain were allowed to re-enter from the inlet. Simulations continued up to 7.4 μs, which was long enough for some particles to travel through the entire domain 8 times, referred to as 8 loops here. Flow trajectories of some of the large and small particles are shown in Supplementary Information-6. Table 1 reports that large circular particles left the flow domain on average 39% more often than small particles in Fig. 10. Thus, as compared to small particles, large particles had smaller residence times and were less-frequently trapped by transient vortices, different from earlier findings 7,8 that relied on the assumption of particle entrapments by steady vortices. However, the use of steady vortex structures to assess particle entrapments may still be valid for microfluidics involving dilute suspensions in lower Re flows.
Finally, the effect of the geometric shape of the large particles on the vortex entrapments of small and large particles were investigated for the microfluidic domain shown in Fig. 9. In DSP-LBM simulations, the shape of the large particles was either circular, elliptical (with an aspect ration of 1.2), or hexagonal with the surface area  of 0.11 μm 2 , while the small particles were circular. This simulation was setup to mimic a small number of large, non-circular tumor cells dispersed in a large number of small, circular healthy cells. Figure 11 shows that the particle shape affected the residence time of all particles in the microfluidic domain. For example, although Particle 38 was permanently trapped in the microfluidic domain if the large particles were circular, it traveled through the microfluidic domain 8 times if the large particles were hexagonal. Moreover, Table 1 shows that large hexagonal particles, when compared to the simulation with large circular particles, resulted in shorter average residence times for all particles with 5% and 7% increases in the number of loops for small and large particles, respectively. Strikingly, the use of large elliptical particles, instead of large circular particles, resulted in 36% enhanced entrapments for smaller particles, while 21% less entrapments for larger particles. Although these findings require further systematic experimental and numerical analyses to confirm, DSP-LBM simulations showed for the first time that by changing the shape of large particles from circular to elliptical, the smaller particles could be selectively entrapped by transient vortices while the larger particles could be effectively flushed out, which is in contrast to current and proposed uses of microfluidics for vortex-controlled, size-based separation of rigid particles.
In brief, considering strong disparities between flow trajectories of the circular and non-circular particles in microfluidic domains, the use of surrogate spherical particles to mimic tumor cells of abnormal shapes 10,11 in microfluidic experiments as in ref. 7 could lead to misleading assessments on the performance of the microfluidic designs proposed to isolate CTCs from healthy cells in biofluids. Here, we demonstrated that DSP-LBM could serve as a useful numerical tool for such analyses.  Figure 11. Number of loops (trips) each particle experienced across the microfluidic domain.