www.nature.com/scientificreports OPEN Particle Shape Influences Settling and Sorting Behavior in Microfluidic Domains Received: 18 October 2017 1 2 3 4 Hakan Başağaoğlu , Sauro Succi , Danielle Wyrick & Justin Blount Accepted: 18 May 2018 We present a new numerical model to simulate settling trajectories of discretized individual or a mixture Published: xx xx xxxx 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 . For example, ellipsoidal microparticles displayed longer blood circulation times than spherical particles due to less efficient phagocytosis by macrophages in the reticuloendothelial system . Hexagonal nanoparticles more effec- tively mitigated phagocytoses and remained in blood circulation longer than spherical particles . Unlike spherical particles, boomerang-shaped particles displayed a preferred direction of Brownian motion , 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 . 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 . In applications relevant to the design of biomedical devices, microuidic de fl vices with different geometric designs have been proposed to isolate circulating tumor cells (CTC) from healthy cells in blood samples through, for example, 7,8 vortex-aided particle separation , which could be useful for early cancer diagnosis and monitoring metastatic progres- sion or the efficiency of cancer treatments . Although the performance of the microuidic de fl vices in the segregation of CTC has been commonly tested with surrogate spherical particles, tumor cells oen exhi ft bit patient-specific arbitrary 10,11 shape proles, w fi hich do not conform to the spherical particle representation for tumor cells . The effect of non-spherical particle shapes on particle trajectories has been recently addressed in numerical 12,13 simulations. Settling dynamics and patterns of thin disks in an infinitely long viscous fluid domain and set- tling behaviors of individual spherical, cubical, or tetrahedral particles in an infinitely long fluidic domain with periodic lateral boundaries 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 mix - ture of DSP is very limited or perhaps non-existent in the literature. 1 2 Mechanical Engineering Division, Southwest Research Institute, San Antonio, TX, 78238, USA. Istituto Applicazioni del Calcolo, via dei taurini 19, 00185, Roma, Italy. Space Science Division, Southwest Research Institute, San Antonio, TX, 78238, USA. Defense Intelligence Solutions Division, Southwest Research Institute, San Antonio, TX, 78238, USA. Correspondence and requests for materials should be addressed to H.B. (email: email@example.com) ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 1 www.nature.com/scientificreports/ e ext Th ension of the lattice Boltzmann (LB) method for simulating flow of suspended bodies is a fast-growing 15 16,17 area of LB research , following the pioneering work of of Ladd . Considering broad uses of DSP in biomedi- cal 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 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. e p Th rimary purpose of this paper is to introduce the DSP-LBM and demonstrate its performance in simulat- ing the settling or flow of individual or a mixture of DSP under various combinations of properties associated with the particles, flow regimes, and the microuidic do fl main geometry. Using the DSP-LBM and a single chamber of the microuidic de fl vice geometry in ref. we numerically investigated the validity of recent findings and implica- tions in microuidic r fl esearch. These findings and implications involve: (i) when a large number of particles are released into a fluid in a microuidic de fl vice, 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 microuidic de fl vices with different geometric peculiarities. Methods 19–22 In the LB method , the mesodynamics of the Newtonian uid flo fl w can be described by a single relaxation time via the Bhatnagar-Gross-Krook (BKG) equation Δt eq ft (, re +Δ tt +Δ )( −= ft rr ,) [( ft ,) − ft (, r )], i ii i (1) where f (r, t) is the set of population densities of discrete velocities e at position r and discrete time t with a time i i eq 24 eq 2 increment of Δt, τ is the relaxation parameter, and is the local equilibrium , f fc =+ ωρ[1 () eu ⋅+ / ii s i i 24 2 () eu ⋅− )/2( cc uu ⋅ )/2], ω is the weight associated with e and c is the speed of sound, c =Δxt /( 3)Δ . The is s i i s s local fluid density, ρ, and velocity, u, at the lattice node are given by ρ = f and ρue = f + τρg, where g is ∑ ∑ i i i i 25 21 the strength of an external force . A D2Q9 (two-dimensional nine velocity vector) lattice 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 ∂ u + (u ⋅ ∇)u = −(∇P/ρ) + ν∇ u + g with the fluid kinematic viscosity, ν =Δ ct(0 τ −.5). Pressure, P, is computed via the ideal gas relation, . Pc = ρ e ext Th ension 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 sur- faces of a circular-cylindrical particle (hereaer ft , circular particle) by Ladd , we used discretized particle surfaces for 2D curved (e.g., circular)- and angular (e.g., hexagonal)-shaped particles in the DSP-LBM. A schematic rep- resentation 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 par- ticles 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 encloses the star-shape and the smaller circle with a radius of R that passes through the corners of the pentagon. es Th e two circles are related via R = ψR , in which ψ = cos(π/5) + [sin(π/5)]/[tan(π/10)]. The surface area of the S P star-shaped particle, A , is given by A = χ(R ) , in which χ =+ [( sint πψ /5)/ ][5/( an πψ /10] 4). The star-shaped S S S particle has five vertices located on the outermost tip of the triangles (v − v ), in addition to five vertices located S1 S5 1,5 on the corners of the inner pentagon (v − v ) (Fig. 1a). The coordinates (x , y ) of v , and v , where iε , are P1 P5 i i Si Pi computed by Eq. 2 and Eq. 3, respectively, x x cosi (( απ ˆ +− 21)/5) i c = + R , y y i c sini (( απ ˆ +− 21)/5) (2) x x cosi (2 απ ˆ +− (1)/5) i c = + R , y y sini (2 απ ˆ +− (1)/5) i c (3) α ˆ where is the initial tilt angle of the particle in the clockwise direction. x = (x , y ) is the center of mass of a par- c c c 1 1 N N ticle, and , where N is the number of vertices (N ) for DAsP or the number of bound- x = x yy = ∑ ∑ Ver c 1 i 1 i= c i= i N N ary nodes (N ) for DCsP. The mass of the star-shaped particle per unit particle thickness is given by Bnd m = χ ( R ) ρ . The moment of inertia, I , for the star-shaped particle was computed by p S p s 22 22 2 , in which a is the side length of the IA = (/ ρρ ac 24) 13 + ot ++ (5Ah /72)(4 3) aA ++ 2(λσ h/3) sP () T T p p pentagon, a = 2R sin(π/5), h is the height of an isosceles triangle, h = a/[2tan(π/10)], A is the area of the P P ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 2 www.nature.com/scientificreports/ Figure 1. A schematic representation of non-circular particle geometries in the DSP-LBM. α ˆ > 0 and α > 0 represent the initial tilt angle in the clockwise and counterclockwise directions. α ˆ =° 0 in (a). pentagon, Aa =+ (1/4)5 (5 + 25 ) , A is the area of the triangle, A = ah/2, λ = R cos(π/5), and σ = [cos( T T p 2 2 π/5) + cos(2π/5)] + [0.5 + sin(π/5) + sin(2π/5)] . Die ff rent 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 x x cosc () Φ− os() αα ˆˆ sins () Φ in() i c c/2 ii = + y y coss () ΦΦ in() αα ˆˆ sinc () os() d/2 i c ii (4) in which Φ = 2π(i − 1)/(N − 1). The mass of an elliptical particle per unit particle thickness is given by i Nbd m = A ρ , in which the surface area and its moment of inertia are computed by A = πcd/4 and , Ic =+ () d p E p E respectively. Intra-Particle Boundary Nodes (IPBN) and Extra-Particle Boundary Nodes (EPBN). e Th wind - ing number algorithm was implemented to determine whether a lattice node x = (x , y ) is enclosed by a pol- k k k ygon 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 , which is referred to as the winding number, m(x ). x is not enclosed by a polygon if m(x ) = 0. In the DSP-LBM, x and x + e form a intra-particle k k k k k i boundary nodes (IPBN) and extra-particle boundary nodes (EPBN) pair if m(x ) ≠ 0 and m(x + e ) = 0. The k k i 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 18,27 between the fluid and the mobile DSP, following the approach in refs in which the population densities near particle surfaces are modified to account for momentum-conserving particle-fluid collisions. Particle-fluid hydrodynamic forces, F , at the boundary nodes located halfway between the intra-particle lattice node, r , and 16,28,29 extra-particle lattice node, r + e , are computed by v i ρω ′ ⁎ i Fr =− 2( ft +Δ eu ,) t +⋅ () ee . rr vi ii i 2 bb (5) ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 3 www.nature.com/scientificreports/ Figure 2. (a) IPBNs and EPBNs of a discretized hexagonal particle geometry in the DPS-LBM. Blue lines are the hydrodynamic links along which the particle and fluid exchange momentum. (b) Momentum exchange between the particle and fluid at a boundary node marked by a square. The translational velocity, U , and the angular velocity of the particle, Ω , are advanced in time according to p p ΔΔ t t the discretized Newton’s equations of motion, and U () tt +Δ ≡+ UF () tt () +− () ρρ g pp T m ρ Δt , where m is the particle mass, I is the moment of inertia of the particle, and Ω () tt +Δ ≡Ω () tt + T () p p pp T u = U + Ω × ( r − r ). The new position of the center of mass of a particle is computed as b p p b c x (t + Δt) = x (t) + U (t)Δt. The population densities at r and r + e Δt are updated to account for particle-fluid c c p v v i hydrodynamics in accordance with 2ρω 2ρω ⁎ i ⁎ i ft ′+ (, rr Δ= tf )( ,) t −⋅ () ue ,( ft re +Δ ,) tt +Δ =′ ft (, re +Δ t ) +⋅ () ue . iv v r i vi iv i r i i 2 b i 2 b c c (6) s s New Locations of Vertices or Boundary Nodes. The locations of vertices or boundary nodes are th updated in each time step. The distance d = (d , d ) between the i vertex (or a boundary node) and the center i ix iy of mass of a particle, x is computed via d = x − x . After x (t + Δt) is computed, new positions of vertices (or c i i c c boundary nodes) are updated via dt cos(() Ω+ ϒΔ ) x x ix pi i c = + y y dt sin(() Ω+ ϒΔ ) i c iy pi (7) in which ϒ is the angle between (x − x ) and +x. For a hexagonal particle, for example, ϒ = α + (i − 1)π/3 for i i c i 1,6 iε . Model Validation. e DS Th P-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 . at two different Reynolds numbers, Re = 8.33 and Re = 1.03 (here, R = 2RU /ν, where R is the particle radius and U is the settling (terminal) velocity of the e s s particle). In ref. 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 o fl w domain was set to ∼30 W (adopted in all settling sim- ulations in this paper), where W is the channel width perpendicular to the main settling direction, and 2 2 30 R = 385 μm, ν = 0.01 cm , and |g| = 981 cm/s . ρ /ρ was adjusted to meet the reported Re values in ref. . For Re = 8.33, DSP-LBM (with ρ /ρ = 1.07) and FE solutions are in good agreement (Fig. 3b), although the DSP-LBM solution for Re = 6.65 (with ρ /ρ = 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 ρ /ρ = 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 of an elliptical particle in an initially quiescent fluid in a bounded domain (Fig. 3d) was com- pared against numerical solutions by Xia et al. . In these simulations, c/d = 2, W/c = 4 (Fig. 3d), density ratio of 2 2 α ˆ =° 45 ρ /ρ = 1.1, , ν = 0.01 cm /s, c = 0.1 cm, and |g| = 981 cm/s . 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. . Data availability. e d Th atasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 4 www.nature.com/scientificreports/ Figure 3. Numerical validations of the DSP-LBM with two benchmark problems, involving settling of a circular particle in (a–c) and an elliptical particle in (d–f). Figure 4. A schematic representation of initial orientations of non-circular particles in the settling simulation. e cen Th ter of mass of the particles was initially located on the mid-channel (y = W/2) near the inlet. 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. e s Th ame 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 , in which R is the equivalent radius of a circular particle that has the same surface area of a non-circular particle. −2 −3 2 2 In these simulations, R = 3.5 × 10 cm, the surface area of the particle is A = 3.9×10 cm , and g = 981 cm/s . e p The same A was specified for all DSP by setting R = 15.4, and R = 5.9 for the star particle; B =15, ζ = π/3, p s p φ = π/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 = 9.2 for the circular particle; and c = 26, d = 13 for the ellip- −3 tical particle (Fig. 1). Here, the length parameters are expressed in l.u. (1 l.u. = 3.846 × 10 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 con- fined channel for ρ /ρ = 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) aer ft 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 gravita- tional 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 ρ/ρ 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 ρ/ρ = 1.10, the small chip was responsible for the persistent periodicity in the boomerang particle’s settling trajectory, differ - ent 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. ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 5 www.nature.com/scientificreports/ Figure 5. Settling trajectories and cumulative angular rotations of DSP with ρ/ρ = 1.05 (a–b) and with ρ/ρ = 1.10 (c–d). W/R = 11.3 in these simulations. p e Figure 6. Settling trajectories of (a) boomerang particles with different surface concavity of its trailing edge −3 2 (simulations were performed with ρ/ρ = 1.10, W/R = 11.3, ζ = 60° and A = 3.9 × 10 cm ), and (b) rectangular p e particles of different aspect ratios (simulations were performed with ρ/ρ = 1.10 and W/R = 9.2). p e e o Th ther 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 ρ/ρ = 1.10, −3 2 W/R = 11.3, ζ = 60°, and A = 3.9×10 cm . DSP-LBM simulations show that the boomerang particle displayed e p gradually vanishing oscillations in its settling trajectory, similar to the triangular particle, if x was located inside the polygonal surface (for φ = 10° and 20°). The boomerang particle exhibited periodic oscillations in its settling trajectory if x was located on the polygonal surface (for φ = 30°) or outside the polygonal surface (for φ = 40°). −1 −1 The oscillation frequency,ϑ , dropped from 1.27 s to 1.13 s as x 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 ρ/ρ = 1.10 and W/R = 9.2. Although rectangular particles with different aspect ratios drifted toward the same p e ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 6 www.nature.com/scientificreports/ Figure 7. (a) Flow trajectories and (b) cumulative angular rotations of DSP at Re = 35. W/R = 11.3 and ρ /ρ = 1.0. equilibrium position at the centerline at x/6 W ∼ , 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 trajec- tory, 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 aer t ft he steady-flow field was established. A neutrally-buoyant spherical particle in a Poiseuille flow typically exhibits the Segre-Silberberg effect with an equilibrium settling position between the channel wall and centerline, in which the 33,34 equilibrium position varies with Re (here, Re = 2RU /ν, where U is the average steady fluid velocity prior to ss ss releases of the particles). For Re ∼. 01 and W/R = 6.6, the equilibrium position of the neutrally-buoyant spherical 33,35 particle was on the centerline in a tube . Consistent with these findings, different equilibrium positions of a circu- lar 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 parti- cle 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 . 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 x/2 W ∼ 5 their equilibrium position at . 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 boo- merang 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. e eff Th ect 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 par - ticles, regardless of their shapes, had the same surface area with R = 385 μm. The interparticle distance at the release location was 4R and the width and length of the domain was 40R × 80R . The fluidic domain was bounded in the e e e 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 , 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 veloci- ties 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 ∼09 .− 12 . . Accurate displacements and velocities are crit- ical in the design of engineered particles for targeted drug deliveries. Figure 8 demonstrated that non-circular shapes of particles have pronounced effects on the settling and flow behaviors of a mixture of DSP and their rep- resentation by circular shapes introduces errors in calculations of particles trajectories and travel times. ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 7 www.nature.com/scientificreports/ Figure 8. (a) Comparison of settling trajectories of seven DSP (in solid lines) to seven circular particles (in dashed lines) in an initially quiescent fluid in a confined domain, (b) the ratio of the settling velocity of a different-shape particle to its circular-shape counterpart released from the same point. (c) Comparison of flow trajectories of seven DSP (in solid lines) to seven circular particles (in dashed lines) in a Poiseuille flow with Re = 38, (d) the ratio of the translational velocity of a different-shaped particle to its circular-shape counterpart released from the same point. t is the total simulation time. 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 microuidic a fl nalyses: (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 uid fl 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 microuidics? T fl o answer these questions, DSP-LBM simulations were setup using a single chamber of the microuidic g fl eometry in ref. . 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 /s and c = 1,460 m/s, and the particles were neutrally buoyant. The average flow rate, u , of 52.14 m/s at the inlet in a single-chambered avg microuidic c fl hamber produced vortex structures, similar to the vortex structures in a multi-chambered micro- uidic de fl vice with m/s in Fig. 3 of ref. . u ∼ 1, 700 avg 7,8 Although steady vortex structures were previously envisioned to trap particles in microfluidic devices , Fig. 10 shows that vortices in a flowing fluid including mobile particles are indeed unsteady, even if the pres- sure differential at the inlet and outlet is held constant in time. Symmetry breaks in the flow domain with ini- tially 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 , 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. the flow field ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 8 www.nature.com/scientificreports/ Figure 9. Steady-flow field in a subsection of microuidic g fl eometry in ref. . All dimensions are scaled with respect to the large particle diameter, D. Small circles attached to particles are used to trace angular rotations of particles. Figure 10. Transient vortex structures. 7,8 involving multiple mobile particles was inherently transient, which contradicts the use of steady vortex regions in assessing particle entrapments in microuidics de fl vices. Moreover, the sorting mechanism related to correla- tions between the lateral displacements of particles to their sizes 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 oen t ft han small particles in Fig. 10. Thus, as compared to small particles, large particles had smaller resi- 7,8 dence times and were less-frequently trapped by transient vortices, different from earlier findings 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 microuidics in fl volving 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 microuidic do fl main 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 ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 9 www.nature.com/scientificreports/ Geometric Total Number of Average Number Total Number of Average Number Shape of LPs* Loops by LPs of Loops by a LP Loops by SPs* of Loops by a SP Circular 39 3.9 84 2.8 Elliptical 47 4.7 54 1.8 Hexagonal 41 4.1 90 3.0 Table 1. Number of trips (loops) the particles experience in a microuidic de fl vice. (*) LP stands for large particles of different geometric shapes. SP stands for small circular-cylindrical particles. Figure 11. Number of loops (trips) each particle experienced across the microuidic do fl main. of 0.11 μm , 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 par- ticle shape ae ff cted the residence time of all particles in the microu fl idic domain. For example, although Particle 38 was permanently trapped in the microuidic do fl main if the large particles were circular, it traveled through the microuidic do fl main 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 entrap- ments for smaller particles, while 21% less entrapments for larger particles. Although these findings require fur - ther 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 microuidics f fl or vortex-controlled, size-based separation of rigid particles. In brief, considering strong disparities between flow trajectories of the circular and non-circular particles in 10,11 microuidic do fl mains, the use of surrogate spherical particles to mimic tumor cells of abnormal shapes in microuidic fl experiments as in ref. could lead to misleading assessments on the performance of the microuidic fl designs proposed to isolate CTCs from healthy cells in biou fl ids. Here, we demonstrated that DSP-LBM could serve as a useful numerical tool for such analyses. References 1. Champion, J. A., Katare, Y. K. & Mitragotri, S. Particle shape: A new drug design parameter for micro- and nano-scale drug delivery. J. Control. Release 121, 3–9 (2007). 2. Sharma, G. et al. Polymer particle shape independently influences binding and internalization by macrophages. J. Control. Release 147, 408–412 (2010). 3. Lin, S. Y., Hsu, W. H., Lo, J. M., Tsai, H. C. & Hsiue, G. H. Novel geometry type of nanocarriers mitigated the phagocytosis for drug delivery. J. Control. Release 154, 84–92 (2011). 4. Chakrabarty, A. et al. Brownian motion of boomerang colloidal particles. Phys. Rev. Lett. 111, 160603 (2013). 5. Fan, Z., Senapati, D., Singh, A. K. & Ray, P. C. Theranostic magnetic coreâ€“plasmonic shell star shape nanoparticle for the isolation of targeted rare tumor cells from whole blood, fluorescence imaging, and photothermal destruction of cancer. Mol. Pharm. 10, 857–866 (2013). 6. Shah, S., Liu, Y., Hu, W. & Guo, J. Modeling particle shape-dependent dynamics in nanomedicine. J. Nanosci. Nanotechnol. 11, 919–928 (2011). 7. Paié, P., Che, J. & Carlo, D. D. Effect of reservoir geometry on vortex trapping of cancer cells. Microfluid. Nanofluid. 21, 104 (2017). 8. Zhou, J., Kasper, S. & Papautsky, I. Enhanced size-dependent trapping of particles using microvortices. Microfluid. Nanofluid. 15, 611–623 (2013). 9. Alixis-Panabières, C. & Pantel, K. Circulating tumor cells: liquid biopsy of cancer. Clin. Chem. 59, 110–118 (2013). 10. Park, S. et al. Morphological differences between circulating tumor cells from prostate cancer patients and cultured prostate cancer cells. PLoS One 9, e85264 (2014). ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 10 www.nature.com/scientificreports/ 11. Marrinucci, D. et al. Cytomorphology of Circulating Colorectal Tumor Cells:A Small Case Series. J. Oncol. 2010, 861341 (2010). 12. Auguste, F., Magnaudet, J. & Fabre, D. Falling styles of disks. J. Fluid Mech. 719, 388–405 (2013). 13. Chrust, M., Bouchet, G. & Dušek, G. Numerical simulation of the dynamics of freely falling discs. Phys. Fluids 25, 044102 (2013). 14. Rahmani, M. & Wachs, A. Free falling and rising of spherical and angular particles. Phys. Fluids 26, 083301 (2014). 15. Succi, S. Lattice Boltzmann 2038. Europhys. Lett. 109, 50001 (2015). 16. Ladd, A. J. C. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation. J. Fluid Mech. 271, 285–309 (1994). 17. Ladd, A. J. C. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results. J. Fluid Mech. 271, 311–339 (1994). 18. Nguyen, N.-Q. & Ladd, A. J. C. Lubrication corrections for lattice-Boltzmann simulations of particle suspensions. Phys. Rev. E 66, 046708 (2002). 19. Higuera, F. J. & Succi, S. Simulating the flow around a circular cylinder with a lattice Boltzmann equation. Europhys. Lett. 8, 517–521 (1989). 20. Benzi, R., Succi, S. & Vergassola, M. The lattice-Boltzmann equation: Theory and applications. Phys. Rep. 222, 145–197 (1992). 21. Succi, S. The lattice-Boltzmann Equation for Fluid Dynamics and Beyond, New York 2001. 22. Wolf-Gladrow, D. A. Lattice Gas Cellular Automata and Lattice Boltzmann Models, Springer-Verlag, Berlin Heidelberg 2000. 23. Bhatnagar, P. L., Gross, E. P. & Krook, M. A. A model for collision process in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94, 511–525 (1954). 24. Qian, Y. H., D’Humieres, D. & Lallemand, P. Lattice BGK models for Navier-Stokes equation. Europhys. Lett. 17, 479–484 (1992). 25. Buick, J. M. & Greated, C. A. Gravity in a lattice Boltzmann model. Phys. Rev. E 61, 5307–5320 (2000). 26. O’Rourke, J. Point in Polygon, in Computational Geometry, 2nd Edition 1998. 27. Başağaoğlu, H. & Succi, S. Lattice-Boltzmann simulations of repulsive particle-particle and particle-wall interactions: Coughing and choking. J. Chem. Phys. 132, 134111 (2010). 28. Aidun, C. K., Lu, Y. & Ding, E.-J. Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation. J. Fluid Mech. 373, 287–311 (1998). 29. Başağaoğlu, H., Carrola, J. T. Jr., Freitas, C. J., Başağaoğlu, B. & Succi, S. Lattice Boltzmann simulations of vortex entrapment of particles in a microchannel with curved and flat edges. Microfluid. Nanofluid. 18, 1165–1175 (2015). 30. Feng, J., Hu, H. H. & Joseph, D. D. Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid Part 1. Sedimentation. J. Fluid Mech. 261, 95–134 (1994). 31. Xia, Z. et al. Flow patterns in the sedimentation of an elliptical particle. J. Fluid Mech. 625, 249–272 (2009). 32. Segré, G. & Silberberg, A. Behavior of macroscopic rigid spheres in Poiseuille flow. J. Fluid Mech. 14, 136–157 (1962). 33. Karnis, A., Goldsmith, H. L. & Mason, S. G. The flow of suspensions through tubes: V. Inertial effects. Can. J. Chem. Engng. 44, 181–193 (1966). 34. Matas, J.-P., Morris, J. F. & Guazzelli, E. Inertial migration of rigid spherical particles in Poiseuille flow. J. Fluid Mech. 515, 171–195 (2004). 35. Yang, B. H. et al. Migration of a sphere in tube flow. J. Fluid Mech. 540, 109–131 (2005). 36. Asmolov, E. S. The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number. J. Fluid Mech. 381, 63–87 (1999). 37. Huhn, F. et al. Quantitative flow analysis of swimming dynamics with coherent Lagrangian vortices. Chaos 25, 087405 (2015). 38. Sajeesh, P. & Sen, A. K. Particle separation and sorting in microuidic de fl vices: a review. Microfluid. Nanofluid. 17, 1–52 (2014). Acknowledgements Funding for this research was provided by Southwest Research Institute’s Internal Research and Development Program, 18R-8602 and 15R-8651. S.S. wishes to acknowledge funding from the European Research Council under the European Union’s Horizon 2020 Framework Programme (No. FP/2014- 2020)/ERC Grant Agreement No. 739964 (COPMAT). The authors thank Miriam R. Juckett of Southwest Research Institute for reviewing the manuscript. Author Contributions H.B. and S.S. developed the numerical model. H.B. ran numerical simulations. D.W. helped with analyzing the results, J.B. involved in coding and simulation runs. All authors reviewed the manuscript. Additional Information Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-018-26786-7. Competing Interests: The authors declare no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre- ative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not per- mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2018 ScIenTIfIc REPO R ts | (2018) 8:8583 | DOI:10.1038/s41598-018-26786-7 11
Scientific Reports – Springer Journals
Published: Jun 5, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
All the latest content is available, no embargo periods.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud