Combined effects of fluid type and particle shape on particles flow in microfluidic platforms
Hakan Başağaoğlu, Justin Blount, Sauro Succi, Christopher J. Freitas
NNoname manuscript No. (will be inserted by the editor)
Combined Effects of Fluid Type and Particle Shape on Particles Flow inMicrofluidic Platforms
Hakan Ba¸sa˘gao˘glu · Justin Blount · Sauro Succi · Christopher J. Freitas
Received: date / Accepted: date
Abstract
Recent numerical analyses to optimize the design of microfluidic devices for more effective en-trapment or segregation of surrogate circulating tumor cells (CTCs) from healthy cells have been reportedin the literature without concurrently accommodating the non-Newtonian nature of the body fluid and thenon-uniform geometric shapes of the CTCs. Through a series of two-dimensional proof-of-concept simulationswith increased levels of complexity (e.g., number of particles, inline obstacles), we investigated the validity ofthe assumptions of the Newtonian fluid behavior for pseudoplastic fluids and the circular particle shape fordifferent-shaped particles (DSP) in the context of microfluidics-facilitated shape-based segregation of parti-cles. Simulations with a single DSP revealed that even in the absence of internal geometric complexities ofa microfluidics channel, the aforementioned assumptions led to 0.11-0.21 W ( W is the channel length) errorsin lateral displacements of DSPs, up to 3-20% errors in their velocities, and 3-5% errors in their travel times.When these assumptions were applied in simulations involving multiple DSPs in inertial microfluidics withinline obstacles, errors in the lateral displacements of DSPs were as high as 0.78 W and in their travel timesup to 23%, which led to different (un)symmetric flow and segregation patterns of DSPs. Thus, the fluid typeand particle shape should be included in numerical models and experiments to assess the performance ofmicrofluidics for targeted cell (e.g., CTCs) harvesting. Keywords
Computational methods in fluid dynamics · Hydrodynamics, hydraulics, hydrostatics
PACS
PACS 47.11.-j · Microfluidic devices with distinct geometric peculiarities have been proposed and tested for size-based and/orshape-based segregation of targeted cells in diverse applications. A microfluidic device with a narrow channelconnected to an expanded region with multiple outlets was used to sort out
Euglena gracilis , microalgeaexplored for biodiesel and biomass production, based on their geometric shapes with different cell aspectratios [1]. Similarly, various microfluidics methods and geometric designs [2,3] have been developed in cancerresearch to segregate rare circulating tumor cells (CTCs), typically occuring 0–10 CTCs/mL of blood[4,
H. Ba¸sa˘gao˘gluMechanical Engineering Division, Southwest Research Institute, San Antonio, TX 78238 USA; Currently at Edwards AquiferAuthority, San Antonio, TX 78215Tel. : +1-210-4775105, email: [email protected]. BlountDefense Intelligence Solutions Division, Southwest Research Institute, San Antonio, TX 78238 USAS. SucciFondazione Istituto Italiano di Tecnologia, Center for Life Nanoscience at la Sapienza, Rome, vle Regina Margherita, 00165-Italy; Istituto Applicazioni del Calcolo, Via dei Taurini 19, 00185, Roma, ItalyC. J. FreitasMechanical Engineering Division, Southwest Research Institute, San Antonio, TX 78238 USA a r X i v : . [ phy s i c s . c o m p - ph ] J un Hakan Ba¸sa˘gao˘glu et al. simultaneously accommodate the non-Newtonian behavior ofthe fluid flow and non-uniform shapes of particles. Although the fluid was assumed to be Newtonian inassessing the performance of microfluidic cell capture or segregate devices in [13,21,23,24,25,26], potentialerrors associated with this assumption have not been reported to date.Thus, the main motivation of this paper is to assess potential errors in trajectories and velocities ofa mixture of non-uniform shaped particles in a pseudoplastic fluid in microchannels, if the pseudoplasticfluid is approximated by a Newtonian fluid and the geometric shape of the particles is assumed to becircular. This assessment is crucial as the performance of microfluidic cell capture or segregation deviceshas been commonly tested using a Newtonian fluid as in [13,20,21,23,24,25,26]. The secondary motivationis to present a new numerical model for simulating flow of a mixture of DSPs in a non-Newtonian fluidin microfluidics with complex geometrical features, which is suitable for simulating fate and transport ofCTCs in body fluid. Unlike the particles flow model that simulates particles as thin solid shells filled with aviscous fluid [23], the DSP-LBM simulates particles with intra-particle non-viscous ‘ghost’ fluid that does notcontribute to particle-fluid hydrodynamics. Therefore, the DSP-LBM is readily suitable for simulating high-frequency rotations of settling or flowing multiple non-uniform shaped particles, including both discretizedcurve-shaped and angular-shaped particles [14].In this paper, we report for the first time the DSP-LBM simulations of a mixture of DSPs in non-Newtonian fluid flow in a microflow channel with an array of inline obstacles. Because CTCs are less de-formable than white blood cells [25], DSP-LBM simulations focused on the behavior of mixture of rigid,non-uniform shaped particles in non-Newtonian fluid flow in microchannels. 2D DSP-LBM simulation withNewtonian and non-Newtonian fluids were performed and compared to investigate the effect of (i) non-Newtonian fluid behavior (pseudoplastic or dilatant) on the lateral displacements of an individual DSP in luid Type and Particle Shape Affect Trajectories of Particles in Microfluidic Platforms 3 a microchannel; (ii) flow strength, inertial focusing, fluid type (pseudoplastic vs. Newtonian), and particleshape on the travel times and flow trajectories of a mixture of DSPs in a microchannel without inline ob-stacles; and (iii) the order of equally-spaced DSPs released from multiple-ports near the inlet on the flowtrajectories of DSPs in a Newtonian or pseudoplastic fluid in inertial microfluidics with I-shaped inline ob-stacles. DSP-LBM simulations demonstrated that different geometric shapes of the particles (i.e., surrogatecells) and the non-Newtonian behavior of the fluid should be accommodated in microfluidic experimentsand numerical simulations to assess or optimize the performance microfluidic device geometries for enhancedsize-based and shape-based cell enrichment from body fluids.
In the LBM [27,28,29,30], the mesodynamics of the incompressible, non-Newtonian fluid flow [31,32,33,34]can be described by a single relaxation time (SRT) [35] f i ( r + e i (cid:52) t, t + (cid:52) t ) − f i ( r , t ) = (cid:52) tτ ∗ [ f eqi ( r , t ) − f i ( r , t )] , (1)where f i ( r , t ) is the complete set of population densities of discrete velocities e i at position r and discretetime t with a time increment of (cid:52) t , τ ∗ is the relaxation parameter associated with non-Newtonian fluid flow,and f eqi is the local equilibrium [36], described by f eqi = ω i ρ (cid:18) e i · u c s + ( e i · u ) c s − u · u c s (cid:19) . (2)where ω i is the weight associated with e i and c s is the speed of sound, c s = (cid:52) x/ (cid:112) (3) (cid:52) t , and the local fluiddensity, ρ , and velocity, u , at the lattice node are given by ρ = (cid:80) i f i and ρ u = (cid:80) i f i e i + τ ∗ ρ g , where g is thestrength of an external force [37] and τ ∗ = 0 . ν ∗ (cid:0) (cid:52) t/ (cid:52) x (cid:1) . The kinematic viscosity of the non-Newtonianfluid is described as ν ∗ = (cid:104) n − | Π D | n − (cid:105) ξ [34,38], in which ξ is the consistency, Π D is the second invariantof the rate of strain tensor, n is the fluid-type identifier, n < n = 1, and n > Π D is computed as Π D = 12 (cid:16) [ tr ( D )] − tr ( D ) (cid:17) , (3)where D = (cid:0) ∇ u ( x ) + ( ∇ u ( x )) T (cid:1) , u = ( u, v ), x = ( x, y ), T is the transpose, and tr is the trace of thematrix D . A D2Q9 (two-dimensional nine velocity vector) lattice [29] was used in numerical simulations. ξ = ν for the Newtonian fluid, in which ν is the kinematic viscosity of the Newtonian fluid. τ ∗ is relatedto ν via τ ∗ = 0 . τ − . ν ∗ /ν , in which τ is the relaxation parameter associated with the Newtonianfluid. Through the Chapman-Enskog approach, the LB method for a single-phase non-Newtonian fluid flowrecovers the Navier-Stokes equation in the limit of small Knudsen number for weakly compressible fluids( (cid:52) ρ/ρ ∼ M ∼ × − , where M is the Mach number), in which ∇ · u ∼ ∂ t u + ( u · ∇ ) u = − ∇ Pρ + ν ∗ ∇ u + g . Pressure differential, P , is computed via the ideal gas relation, P = c s ρ . Different-shaped particles . 2D simulations of flow of particles in Newtonian or non-Newtonian fluids wereperformed in this paper using the DSP-LBM that accommodates particle-fluid hydrodynamics of DSPs [14].Simulations were conducted using discretized angular shaped particles (DAsPs), encompassing rectangularand hexagonal particles, and discretized curved-shaped particles (DCsPs), encompassing circular-cylindricaland elliptical particles. The DSP-LBM calculates first the coordinates of the boundary nodes, ( x i , y i ), onthe circumference of DAsPs, using the information on the center of the mass of a particle x c = ( x c , y c ) andother geometric shape-specific parameters. ( x i , y i ) for circular and elliptical particles are computed by Eqs.4 and 5, respectively, (cid:20) x i y i (cid:21) = (cid:20) x c y c (cid:21) + R p (cid:20) cos (2 π ( i − / ( N bnd − sin (2 π ( i − / ( N bnd − (cid:21) , (4) Hakan Ba¸sa˘gao˘glu et al. (cid:20) x i y i (cid:21) = (cid:20) x c y c (cid:21) + (cid:20) cos ( Φ i ) cos (ˆ α ) − sin ( Φ i ) sin (ˆ α ) cos ( Φ i ) sin (ˆ α ) sin ( Φ i ) cos (ˆ α ) (cid:21) (cid:20) c/ d/ (cid:21) , (5)where R p is the radius of a circular particle, c and d are the length of the major and minor axes ofan elliptical particle, ˆ α is the initial tilt angle of an elliptical particle in the clockwise direction, Φ i =2 π ( i − / ( N Nbd − N Bnd is the number of boundary nodes. Different from DCsPs, the DSP-LBMcalculates first the coordinates of vertices, ( x v , y v ), for DAsPs, based on the information on ( x c , y c ) and othergeometric shape-specific parameters. ( x v , y v ) for hexagonal and rectangular particles are computed by Eqs.6 and 7, respectively, (cid:20) x v Hi y v Hi (cid:21) = (cid:20) x c y c (cid:21) + L (cid:20) cos ( α + ( i − π/ sin ( α + ( i − π/ (cid:21) , (6) (cid:20) x v Rj y v Rj (cid:21) = (cid:20) x c y c (cid:21) + √ l + w (cid:20) ϕ xj ϕ yj (cid:21) , (7)where i(cid:15) [1 , L is the side length of a hexagonal particle, j(cid:15) [1 , l and w are the long and short sidelengths of a rectangular particle, ϕ xj = [ cos ( α + θ ) , − cos ( α − θ ) , − cos ( α + θ ) , cos ( α − θ )] and ϕ yj =[ sin ( α + θ ) , sin ( α − θ ) , sin ( α − θ ) , − sin ( α + θ ) , − sin ( α − θ )], and α is the initial tilt angle of a particle(rectangular or hexagonal) in the counterclockwise direction. x c = N (cid:80) Ni =1 x i and y c = N (cid:80) Ni =1 y i , where N is the number of vertices for DAsPs or the number of boundary nodes for DCsPs. Intra-Particle Boundary Nodes (IPBNs) and Extra-Particle Boundary Nodes (EPBNs) . EachDAsP or DCsP is represented by a polygon in the DSP-LBM. The winding number algorithm[39] is used todetermine the location of lattice nodes inside and closest to the particle surface. These nodes are labeled asIPBNs and denoted by r v . A lattice node outside the particle surface and separated from the closest IPBNby e i is labeled as EPBN and represented by ( r v + e i ). Each ( r v , r v + e i ) pair forms a hydrodynamic linkacross the particle surface along which the mobile DAsP or DCsP exchanges momentum with the bulk fluid.Particle-fluid momentum exchanges occur at all boundary nodes located at r b = 0 . r v + ( r v + e i )]. Particlemotion is determined by local velocities, u r b , computed from particle-fluid hydrodynamics at each r b . Particle-fluid Hydrodynamics . Following the approach in [40,41,42], population densities near particlesurfaces are modified to account for momentum-conserving particle-fluid collisions. Particle-fluid hydrody-namic forces, F r b , at each u r b is computed by F r b = − (cid:20) f (cid:48) i ( r v + e i (cid:52) t, t ∗ ) + ρω i c s ( u r b · e i ) (cid:21) e i . (8)Eq.8 was derived based on the premise that the fluid occupies the entire flow domain to ensure thecontinuity in the flow field to avoid large artificial pressure gradients that may arise from the compres-sion and expansion of the fluid near particle surface; however, the fluid inside the particle does not con-tribute to F r b . The translational velocity, U p , and the angular velocity of the particle, Ω p , are computed by U p ( t + (cid:52) t ) ≡ U p ( t ) + (cid:52) t (cid:104) F T ( t ) m p + ( ρ p − ρ ) ρ p g (cid:105) and Ω p ( t + (cid:52) t ) ≡ Ω p ( t ) + (cid:52) tI p T T ( t ), where F T and T T arethe total hydrodynamic force and torque on the particle exerted by the surrounding fluid, respectively. m p is the particle mass, I p is the moment of inertia of the particle (Table 1), and u b = U p + Ω p × ( r b − r c ).Table 1: Mass and moment of inertia of DSPs. Particle Shape Particle Mass per Unit Particle Thickness, m p Moment of Inertia, I p Circular ∗ πR ρ p (1 / m p R Elliptical π ( cd/ ρ p m p (cid:0) c + d (cid:1) Hexagonal (3 / √ L ρ p ( m p / L (cid:2) cot ( π/ (cid:3) Rectangular lwρ p m p (cid:0) l + h (cid:1) (*) the circular particle is treated as a thin solid disk.luid Type and Particle Shape Affect Trajectories of Particles in Microfluidic Platforms 5 Local forces and torques associated with particle-fluid hydrodynamics at r b , covered/uncovered nodesdue to the particle motion, and steric (repulsive) interaction between particles in close proximity or particlesin close contact with stationary solid objects contribute to F T and T T , F T = (cid:88) r b F r b + (cid:88) r c,ub F r c,ub + (cid:88) | r pw |≤| r it | F r pw + (cid:88) | r pp (cid:48) |≤| r it | F r pp (cid:48) , (9) T T = (cid:88) r b ( r b − r c ) × F r b + (cid:88) r c,ub ( r c,ub − r c ) × F r c,ub + (cid:88) | r pw |≤| r it | ( r w − r c ) × F r pw + (cid:88) | r pp (cid:48) |≤| r it | ( r p (cid:48) − r c ) × F r pp (cid:48) , (10)where F r c,ub = ± ρ (cid:16) u r c,ub − U p (cid:17) / (cid:52) t is the force induced by covered, r cb , and uncovered, r ub lattice nodesdue to due to particle motion [43,44,45]. Steric interaction forces, F r i , between the particles and betweenthe particles and stationary solid zones, including channel walls and inline obstacles, are expressed in termsof two-body Lennard-Jones potentials [14,42] such that F r i = − ψ (cid:16) | r i || r it | (cid:17) − n , where | r i | is the distancebetween a particle surface node and the neighboring particle surface node ( r i = r pp (cid:48) ) or between a particlesurface node and the stationary solid node on channel walls or inline obstacles ( r i = r pw ); p is the particleindex; | r it | is the repulsive threshold distance; n is the unit vector along r i ; and ψ is the repulsive strengthbetween the particles and between the particles and stationary solid nodes.The new position of the center of mass of a particle is computed as x c ( t + (cid:52) t ) = x c ( t ) + U p ( t ) (cid:52) t .The population densities at r v and r v + e i (cid:52) t are updated to account for particle-fluid hydrodynamics inaccordance with [40] f (cid:48) i ( r v , t + (cid:52) t ) = f i ( r v , t ∗ ) − ρω i c s ( u r b · e i ) , (11) f i ( r v + e i (cid:52) t, t + (cid:52) t ) = f (cid:48) i ( r v + e i (cid:52) t, t ∗ ) + 2 ρω i c s ( u r b · e i ) . (12)where f (cid:48) i corresponds to population densities that propagate in − e i after collision and t ∗ represents thepost-collision time. In the end of each time-step, the location of vertices on angular-shaped surfaces orboundary nodes on curved surfaces are updated using the geometrical relations in Eqs. 4-7. The distance d i = ( d ix , d iy ) between the i th vertex (or a boundary node) and x c is computed via d i = x i − x c . After x c ( t + (cid:52) t ) is computed, new positions of vertices (or boundary nodes) are updated by (cid:20) x i y i (cid:21) = (cid:20) x c y c (cid:21) + (cid:20) d ix cos (( Ω p + Υ i ) (cid:52) t ) d iy sin (( Ω p + Υ i ) (cid:52) t ) (cid:21) (13)where Υ i is the angle between x i − x c and + x . . The earlier version of the 2D LBM accommodating circular particles only successfully capturedthe 3D particle flow dynamics (e.g., particle velocities and trajectories) in microfluidic experiments [46],by implementing Reynolds number ( Re )-based dimensional scaling [47], given that a circular-cylindricalparticle would have the same wake structure as the spherical particle, but at a lower Re . The DSP-LBMwas validated in [14] against two benchmark problems, involving previously reported numerically-computedsettling trajectories of a circular particle [47], and the trajectories and angular rotations of a settling ellipticalparticle in an initially quiescent Newtonian fluid [48]. Ref. [14] reported also successful comparisons ofexperimentally-determined [49] and DSP-LBM simulated terminal velocities of a gravity-driven settling ofspherical particles 5% or 10% denser than the bulk fluid. DSP-LBM simulations in [14] were practicallyinsensitive to the grid resolution. The same lattice spacing and particles’ dimensions in [14] were adopted Hakan Ba¸sa˘gao˘glu et al. in DSP-LBM simulations in this paper, in which R p =10 lattice units, N Bnd =100 for DCsPs, and the aspectratio of the elliptical and rectangular particles is 2. N Bnd =100 led to a lattice spacing of ∼ . − . < Non-Newtonian Fluid Flow Module . The generalized analytic solution for the steady velocity profile ofnon-Newtonian ( n (cid:54) = 1) or Newtonian ( n = 1) Poiseuille flow is given by [32,50]: u ( y ) = (cid:18) ξ | g | (cid:19) n (cid:18) W (cid:19) n (cid:18) nn + 1 (cid:19) (cid:34) − (cid:18) | y | W (cid:19) n (cid:35) , (14)where W is the channel width and y is vertical distance, orthogonal to the main flow direction, from one ofthe channel walls. DSP-LBM simulations of flow of a Newtonian fluid with ν =1 cm / s were performed ina channel with W =0.05 cm. Local kinematic viscosities for non-Newtonian fluid flow were computed from ν ∗ = (cid:104) n − | Π D | n − (cid:105) ξ . In these simulations, the maximum steady fluid velocity, u max , was set to 9 cm/s bysetting ξ = 0 . ν for n = 0 . ξ = 1 . ν for n = 1 (Newtonian), and ξ = 1 . ν for n = 1 . ξ is often determined experimentally, in DSP-LBM simulations ξ for eachnon-Newtonian fluid was determined by calibrating the magnitude of the maximum fluid velocity againstthe analytic solution in Eq.14, in which calibration does not affect the shape of the velocity profile, as shownin [38]. DSP-LBM simulations with n = 0 . n = 1 . ν ∗ were set to 10 − and 0 . ν ∗ [31,38].Fig. 1 shows that the steady velocity profiles computed by DSP-LBM closely matched the analyticsolutions given by Eq. 14. The successful model validations, involving particle settling in a Newtonian fluidin [38] and the steady non-Newtonian and Newtonian velocity profiles in Fig. 1, suggest that the SRT isappropriate for DSP-LBM simulations; therefore, the computationally demanding multi-relaxation time wasnot adopted in simulations in this paper. This is consistent with the discussion in [51].Fig. 1: Normalized steady velocity profiles of pseudoplastic ( n = 0 . n = 1), and dilatant( n = 1 .
2) Poiseuille flows. Fluid velocities are normalized with respect to the maximum fluid velocity, u max ,at the midchannel. luid Type and Particle Shape Affect Trajectories of Particles in Microfluidic Platforms 7 The DLD has been implemented to segregate particles based on their shapes or sizes in a Newtonian fluidin microfluidic devices with specifically designed internal geometric features [8,9,10]. The following researchinquiries were investigated in this section: (i) how differently a circular particle would undergo lateral dis-placements after being released into Newtonian or non-Newtonian fluid flow in a microchannel even in theabsence of geometrically complex internal geometric features? and (ii) how would particle trajectories andvelocities vary with the fluid types and particle shapes? The answers to these questions would provide insightsinto whether the assumption of a Newtonian fluid for non-Newtonian fluids and/or the circular particle shapefor non-circular particles would be reliable in optimizing microfluidic designs for entrapment or segregationof arbitrary-shaped CTCs in body fluid.DSP-LBM simulations with a single DSP in Newtonian or non-Newtonian Poiseuille flow were performedto address these inquiries. The channel length, L , (in the direction of the main flow) and width, W , were set to100 R p and 10 R p . A periodic boundary condition was imposed at the inlet and outlet, and a no-slip boundarycondition was imposed along the channel walls. A circular, elliptical, hexagonal, or rectangular particle wasreleased at a point 20% off the midchannel into a pseudoplastic ( n = 0 . n = 1 . n = 1 .
2) fluid with an average steady velocity, u avg , 6.0 cm/s prior to particle release. R p = 500 µ m, ν = 0 . /s (for a Newtonian fluid), and the resultant flow Reynolds number, Re = u avg R p /ν = 30. All DSPs hadthe same surface area of 7 . × − cm . This simulation was set-up such that the circular particle (thereference particle) in a Newtonian fluid drifted to an equilibrium position ∼
8% off the midchannel near theoutlet, exhibiting the S´egre-Silberberg effect that a rigid spherical and neutrally buoyant particles typicallydisplays in slow flow [52,53].In Figs. 2a-d, DSPs displayed steady equilibrium with monotonic drift to the equilibrium position awayfrom the midchannel in pseudoplastic fluid flow, steady equilibrium with transient overshoot in Newtonianfluid flow, and weak to strong oscillatory motion in dilatant fluid flow in a 2D microchannel when u avg wasthe same in all simulations. The circular particle exhibited larger oscillations and transient overshoot duein part to its smaller I p than the other particles, which led to the least resistance to angular rotations. Incontrast, the rectangular and elliptical particles have higher I p ’s, and hence, they exhibit relatively strongerresistance to rotations. Although the hexagonal particle has slightly higher I p than the circular particle, itexhibits persistent rotations as it travels along its equilibrium position in a channel [14] due to asymmetricposition of its vertices about its equilibrium position.Among these simulations, the DSP experienced the highest inertial effects in dilatant fluid flow due tothe largest fluid velocity differentials on its opposite surfaces orthogonal to the main flow direction as itapproached and migrated around the midchannel (Fig. 1). Such transitions in migration modes of DSPs aresimilar to changes in trajectories of gravity-driven settling of a circular particle when the inertial effects wereelevated by increasing the particle density [47,14] or changes in trajectories of neutrally buoyant circular orDSPs flowing in a Newtonian fluid with higher flow strengths [43,46]. Thus, elevated inertial effects on thetrajectories of a mobile particle in a dilatant fluid caused by relatively sharper gradients of fluid velocitiesthan in Newtonian or pseudoplatic fluids are comparable to the elevated inertial effects on the settlingtrajectories of a denser particle in an initially quiescent fluid or flow trajectories of a particle in higher Re -flow. Overall transitions in flow trajectories of a single DSP in pseudoplastic fluid flow to flow trajectories indilatant fluid flow are consistent with settling trajectories of a circular particle [47], or an elliptical particle[48], or a particle of other geometric shapes as the particle density is increased [14].In Figs. 2a-d, after arriving at x/R p ∼
60, the circular particle traveled in the left side of the channel( y/R p >
5) in a dilatant fluid; whereas, it remained in the right side of the channel in Newtonian or pseu-doplastic fluids. If a dilatant fluid was approximated in numerical simulations by or replaced in microfluidicexperiments with a Newtonian fluid, the lateral displacement of the circular particle would be underesti-mated, for example, by a distance of 2 . R p (=0.21 W ) at x/R p ∼ . R p (=0.11 W ) at x/R p ∼ y/R p = 5 at x/R p = 250) in a dilatantfluid, only non-circular particles were drifted to the midchannel in a Newtonian fluid at x/R p = 250. In pseu-doplastic and Newtonian fluid flows, the circular and non-circular particles drifted to different quasi-steady Hakan Ba¸sa˘gao˘glu et al. equilibrium positions at x/R p = 250, revealing the sensitivity of the equilibrium position of a particle to itsgeometric shape.Figs. 2a-d also revealed that the S´egre-Silberberg effect is not only related to the particle size-based Re ,but also to the fluid type and particle shape. Unlike in dilatant fluid flow, all DSPs exhibited the S´egre-Silberberg effect in pseudoplastic fluid flow. However, only the circular particle displayed the S´egre-Silberbergeffect in the Newtonian fluid flow. (a) Particle trajectories (b) Particle trajectories(c) Particle velocities (d) Particle velocities Fig. 2: Trajectories of DSPs in (a) dilatant or Newtonian and (b) Newtonian or pseudoplastic fluid flow. Theratio of velocities of DSPs in (c) dilatant and (d) pseudoplastic fluid flow to their velocities in Newtonianfluid flow. The first subscript of U denotes the particle shape and the second subscript denotes the fluid type(D, N, and Ps corresponds to dilatant, Newtonian, and pseudoplastic fluids). U CD , for example, representsthe velocity of a circular particle in a dilatant fluid.The ratio of the velocity of a particular DSP in a dilatant fluid to its velocity in a Newtonian fluid inFigs. 2c-d shows that approximating a dilatant fluid with a Newtonian fluid would result in errors up to of3-11% (the largest for the elliptical particle and the smallest for the rectangular particle) in local particlevelocities. Similarly, if pseudoplastic fluid flow is approximated by non-Newtonian fluid flow, errors in localparticle velocities would be up to 9-20% (the largest for the rectangular particle and the smallest for thecircular particle). The travel times of the circular, hexagonal, and elliptical particles were 2.8%, 1.1%, and0.8% longer in dilatant fluid flow than in Newtonian fluid flow. In contrast, the rectangular particle traveled1.5% faster in a Newtonian fluid than in a dilatant fluid. On the other hand, the travel times of DSPs were luid Type and Particle Shape Affect Trajectories of Particles in Microfluidic Platforms 9 R p errors in lateral displacements of particles, up to 3-20% errors in particle velocities,and 3-5% errors in travel times of DSPs. DSP-LBM simulations were performed to investigate the combined effects of the fluid type and flow strengthon the flow trajectories and velocities of a mixture of DSPs in a microchannel (Fig. 3) by accommodat-ing interparticle interactions in Eq. 10. Because the blood is considered as a pseudoplastic fluid, only apseudoplastic fluid is used for non-Newtonian fluid flow simulations in the following sections.Fig. 3: A schematic representation of the microfluidic channel geometry and release locations of particles inDSP-LBM simulations. All particles have the same surface area of 7 . × − cm . The initial tilt angle ofnon-circular particles is 0 ◦ . ψ =1 and r it =2.5 lattice spacing in simulating particle-particle and particle-wallinteractions [14].The properties of the fluid and particles, and the specifications of the flow domain boundaries in simu-lations in Section 4 were adopted for simulations in this section, except for the channel width. The channelwidth was widened from 10 R p to 15 R p to initially place 4 particles with the center-to-center separationdistance of 3 R p (Fig. 3). The initial separation distance of 3 R p was chosen to be slightly longer than thelength of the long axis of the elliptical particle (2 . R p ) and the long side of the rectangular particle (2 . R p ).Different from simulations in Section 4, circular, rectangular, elliptical, and hexagonal particles (CREH con-figuration) were simultaneously released from a multiple-port near the inlet into the fluid after the steadyflow field was established. DSP-LBM simulations were performed with a pseudoplastic or Newtonian fluidflowing in a microchannel with an average steady velocities of 6.0 cm/s, corresponding to Re = 30, (slowflow) and 12.0 cm/s, corresponding to Re = 60, (fast flow) prior to releases of DSPs into the fluid. Thissimulations were set up such that the circular particle (the reference particle) drifted monotonically to themidchannel at low inertial effects at Re = 30, but its trajectory displayed transient overshoot about themidchannel at higher inertial effects at Re = 60 [43,47,46].At Re = 30, DSPs drifted toward the midchannel in Fig. 4a due to the combined inertial and wall effects.Particles’ drifts were more pronounced in a Newtonian fluid flow as the velocity gradients of the fluid aroundthe midchannel were sharper in a Newtonian fluid than in a pseudoplastic fluid (Fig. 1). As a result, DSPsdrifted to semi-equilibrium positions off the midchannel in a pseudoplastic fluid. However, at Re = 60 withhigher inertial effects, DSPs displayed oscillatory trajectories with transient overshoot about the midchannel(Fig. 4b) in a Newtonian fluid, in an agreement with the behavior of circular particle in a Newtonian fluid Fig. 4: Flow trajectories of DSPs in pseudoplastic and Newtonian fluids with an average steady velocity of(a) slow flow ( Re = 30) and (b) fast flow ( Re = 60). The ratio of velocities of DSPs in a pseudoplastic fluidto their velocities in a non-Newtonian fluid in (c) slow flow ( Re = 30) and (d) fast flow ( Re = 60).with high inertial effects in [47]. On the contrary, the flow trajectories of DSPs in a pseudoplastic fluid wereless sensitive to the flow strengths considered here. Relatively low velocity gradients orthogonal to the mainflow direction in a pseudoplastic fluid was not sufficient to impose large uneven hydrodynamic forces on theopposite sides of the particle for the particle to gain sufficient angular momentum to drift to the midchannel.Thus, when DSPs were away from walls, they displayed small lateral displacements in a pseudoplastic fluid,regardless of the flow strengths considered.Additional observations from Figs. 4a-d are noteworthy. At Re = 30, the trajectories of the ellipticalparticle in pseudoplastic and Newtonian fluids were similar with the resultant spatial discrepancies in itslateral displacements within R p . Conversely, the trajectories of the remaining particles in pseudoplastic andNewtonian fluids exhibited larger disparities as high as 1.6 R p for the rectangle, 1.8 R p for the hexagonal, and3.3 R p for the circular particles. Therefore, if pseudoplastic fluid flow is approximated by Newtonian fluidflow in these simulations, the maximum error in lateral displacements of the particles would be in the rangeof 0.1 W to 0.2 W . Although the trajectories of the elliptical particle in pseudoplastic and Newtonian fluidswere practically identical at Re = 30, its lateral displacements differed as high as 3.4 R p at Re = 60 due tolarger angular momentum the elliptical particle exhibited at Re = 60. In contrary to the elliptical particle,the trajectories of the rectangular particle in pseudoplastic and Newtonian fluids at Re = 30 and at Re = 60were similar. Because the elliptical and rectangular particles have the same aspect ratio and their initialrelease locations were symmetric about the midchannel, the particles’ shape and their release positions (Fig. luid Type and Particle Shape Affect Trajectories of Particles in Microfluidic Platforms 11
3) were critical for their subsequent migration pathways. Fig. 4c-d shows also that shortly after the particleswere released, only the circular particle traveled slower in a pseudoplastic fluid than in a Newtonian fluid atboth Re = 30 and Re = 60. DSP-LBM simulations were performed to investigate the combined effects of the fluid type, inline obstacles,and the order of particles’ position at the release location on the flow trajectories of a mixture of DSPsin microfluidics with I-shaped inline obstacles. The geometric peculiarities of the microfluidics and inlineobstacles are shown shown in Fig. 5. The properties of the particles and fluids, flow domain boundaries,and the slow flow field condition in simulations in Section 5 were adopted for DSP-LBM simulations inthis section. An array of I-shaped inline obstacles was originally used in [10] for shape-based segregationof particles through the DLD method. Horizontal and vertical separation distances between inline obstaclesin Fig. 5 were chosen to be larger than the long axis of the elliptical particle (2 . R p ) and the long side ofthe rectangular particle (2 . R p ) to avoid particle filtration. The particles in the second multiple port wereinitially rotated by 45 ◦ to reflect potential uncertainties associated with the initial orientations of DSPs. (a) (b) Fig. 5: (a) A schematic representation of the inertial microfluidics geometry and the release locations ofparticles in DSP-LBM simulations, (b) geometric specifics of an I-shaped inline obstacle. R p is the radius ofthe circular particle. All particles have the same surface area of 7 . × − cm .After the steady flow field was established, DSPs were released simultaneously into the fluid near theinlet from two multiple-ports (Fig. 5). We considered two scenarios to investigate the effects of the orderof DSPs at the release location on the particles’ trajectories in pseudoplastic and Newtonian fluid flows. Inthe first scenario, the first multiple-port closer to the inlet involved rectangular, circular, hexagonal, andelliptical particles (the RCHE configuration) with the center-to-center separation distance of 3 R p from theright wall to the left wall. In the second multiple-port, the order of release positions of DSPs was flipped andall DSPs were tilted by 45 ◦ in the clockwise direction. This scenario will be referred to as ‘RCHE’ hereafter,in reference to the particles arrangement in the first multiple-port. In the second scenario, the first multiple-port involved circular, rectangular, elliptical, and hexagonal (the CREH configuration) particles with thecenter-to-center separation distance of 3 R p from the right wall to the left wall. In the second multiple port,the order of release positions of DSPs were flipped and all DSPs were tilted by 45 ◦ in the clockwise direction.This scenario will be referred to as ‘CREH’ hereafter. DSPs in the following discussion are labelled withthe first letter of their geometric shape followed by a number indicating from which multiple-port they arereleased.DSP-LBM simulations in Fig. 6 demonstrate geometric shape-based separation of particles in both pseu-doplastic and Newtonian fluids. In both the CREH and RCHE configurations, most of the particles segregatedtoward the left wall were angular-shaped; whereas, most of the particles segregated toward the right wallwere curve-shaped. However, in pseudoplastic fluid flow, geometric shape-based particle segregation was sen-sitive to the initial configuration of the particles. For example, most of the particles segregated toward the left wall were angular-shaped for the RCHE configuration; whereas, the opposite was true for the CREHconfiguration.Fig. 6 also show that the particles released from the ports closer to the walls kept moving closer tothe walls and avoided flow pathways in between inline obstacles, regardless of the fluid type and geometricshape of the particles. Conversely, the particles released from the ports closer to the midchannel exhibiteddifferent migration patterns in Newtonian and pseudoplastic fluids, which varied also with the geometricshape of the particles. The RCHE configuration in a pseudoplastic fluid (Fig. 6a) resulted in most symmetricparticles’ trajectories, in which three of the particles (E1, R2, C2) were separated toward the left wall and theother three (E2, R1, C1) were separated toward the right wall as they passed through the zone of obstacleswhile the remaining two particles (C2,H2) exhibited nearly symmetric flow trajectories between the inlineobstacles. (a) RCHE in pseudoplastic fluid flow(b) RCHE in Newtonian fluid flow(c) CREH in pseudoplastic fluid flow(d) CREH in Newtonian fluid flow Fig. 6: Flow trajectories of DSPs in pseudoplastic or Newtonian fluids in inertial microfluidics with inlineobstacles until one of the DSPs reached the exit-end. DSPs in two different configurations (RCHE or CREH)were released after the steady flow field was established. The average steady velocity was the same inpseudoplastic and Newtonian fluid flows, corresponding to the slow flow field Re = 30 in Fig. 4. luid Type and Particle Shape Affect Trajectories of Particles in Microfluidic Platforms 13 Unlike in the pseudoplastic fluid flow, the flow trajectories of C2 and H2 in a Newtonian fluid crossedover near the first two obstacles as C2 and H2 separated to the opposite half of the channel with respect totheir release locations. Subsequently, C2 gradually separated toward the right wall as it passed through thezone of obstacles while H2 migrated in between the obstacles. If pseudoplastic fluid flow was approximatedby Newtonian fluid flow, the lateral displacements of C2 and H2 would be off by 6.4 R p (=0.43 W ) and 5.1 R p (=0.34 W ) at x = 75 R p , outside the zone of obstacles. Such deviations would be non-negligible errors inlateral displacements of the particles in inertial microfluidics. Moreover, across the zone of the obstacles,more particles were separated toward the right wall in a Newtonian fluid than in a pseudoplastic fluid. Thecircular particles (C1 in the pseudoplastic fluid flow, and C1 and C2 in Newtonian fluid flow) were the fastestmoving particles. The travel time of the particles in the pseudoplastic and Newtonian fluid flow differed bya factor of 0.91-1.11 for the RCHE and 0.77-1.15 for the CREH configurations.The effects of the order of the particles at the release location on the particles’ trajectories in a pseudo-plastic fluid is evident from Fig. 6a-c. Unlike C2 and H2 in the RCHE configuration, R2 and E2 released fromthe second multiple-port near the midchannel in the CREH configuration did not display symmetric flowtrajectories. Nor did the particles evenly separated to the walls in the RCHE configuration. Two rectangularparticles, R1 and R2, were the fastest moving particles in a pseudoplastic fluid flow in the CREH configu-ration that had the same release positions with the fastest moving C1 and C2 in the RCHE configuration.Thus, the release position of the particles appear to be more critical than their shape in determining themaximum travel time of the particles in pseudoplastic fluid flow in these simulations.The order of the release positions of the particle was also important in Newtonian fluid flow (6b-d), wheremore particles were segregated toward the left wall with negligible subsequent drifts towards the midchannelfor the CREH configuration than for the RCEH configuration. Although the assumption of Newtonian fluidfor the pseudoplastic fluid would lead to small error of ∼ R p (0.1W) in the lateral displacement of E2 at x = 75 R p , the error in the lateral displacements of the R2, R1, and E1 outside the zone of obstacles would beas high as 4.1 R p (=0.27W), 9.9 R p (=0.66W), and 11.6 R p (=0.78W). Hence, both the geometric shape andthe fluid type could have significant effects on the flow trajectories, lateral displacements, and travel times ofDSPs in inertial microfluidics with inline obstacles. Therefore, negligence of the pseudoplastic nature of thefluid and/or the exact geometric shape of the cells (or surrogate particles) would not be practical in assessingor optimizing geometric design of microfluidics proposed or designed for cell segregation. The DSP-LBM withnon-Newtonian fluid flow simulation capabilities circumvents such problems in practice. Recent numerical analysescitePCD17,HC17,KMFD16,MSAM12,JSH15,DTF15,SAM18 to assess the performance of particular microflu-idic device designs in separating CTCs from healthy cells have been reported without concurrently accommo-dating non-circular geometric shapes of CTCs and the non-Newtonian behavior of body fluids. We upgradedthe DSP-LBM [14] for non-Newtonian fluid flow simulations and used it to numerically investigate the reli-ability of the assumptions of a Newtonian fluid for pseudoplastic fluids and circular-shape for non-circularparticles in assessing the performance of microfluidic devices with simplified geometries for shape-basedsegregation of particles.Numerical results demonstrated that the S´egre-Silberberg that the neutrally buoyant particles were pre-viously shown to display in slow flow is not only associated with the flow strength ( Re ), but also with thecombination of particle shape and fluid type. The simulations demonstrated that if a smooth-walled channelfilled with a dilatant, Newtonian, or pseudoplastic fluid flowing with the same average fluid velocity, theDSP would experience higher inertial effects in a dilatant fluid, as opposed to lower inertial effects in apseudoplastic fluid.The results also revealed that the lateral displacements, velocities, and travel times of individual particlesdiffered due to the geometric shape of the particle and the fluid type (Newtonian vs. non-Newtonian). Theaforementioned assumptions resulted in errors as high as 0.21 W in lateral displacements of the particlesin a microchannel. Simulations with a mixture of different-shaped particles in a microchannel showed thatthe lateral displacements of some of the particles in a mixture were practically insensitive to the fluid type.Moreover, unlike a mixture of particles in a Newtonian fluid, the lateral displacements of particles in a pseudoplastic fluid were nearly insensitive to an increased inertial effect. Yet, when these assumptions wereapplied, errors in the lateral displacements of the particles varied in the range of 10-20% of the channelwidth. Although the trajectories of an elliptical particle in a mixture of DSPs was insensitive to the fluidtypes in slow flow, at higher inertial effect its lateral displacements in the Newtonian and pseudoplastic fluidsdiffered by 0.23W.The discrepancy in segregation patterns and lateral displacements of the particles were more pronouncedinertial microfluidics with an array of inline obstacles. Numerical simulations revealed that not only theparticle shape and fluid type, but also the order of the particles at the release location had significanteffect on the particles’ trajectories and their separation patterns around the zone of obstacles. Although theorder and orientations of non-circular particles at the release ports are difficult, if at all possible, to controlin microfluidic experiments, numerical simulations revealed their non-negligible effects on segregation ofparticles across the inertial microfluidics. Errors in lateral displacements of particles would be as high as0.78 W in these simulations, if the aforementioned assumptions are made.In brief, DSP-LBM simulation results demonstrated that the assumption of circular particle shape fornon-circular particles and the Newtonian fluid type for non-Newtonian body fluids are not necessarily reliableand practical in assessing or optimizing microfluidic device designs for segregation of CTCs from healthycells. Microfluidic experiments with a mixture of arbitrary-shaped particles in non-Newtonian fluids underdifferent flow conditions can be used to test our numerical findings. Acknowledgements
Funding for this research was provided by Southwest Research Institute’s Internal Research and Devel-opment Program, 18R-8602 and 15R-8651. SS kindly acknowledges funding from the European Research Council under theEuropean Union’s Horizon 2020 Framework Programme (No. FP/2014- 2020)/ERC Grant Agreement No. 739964 (COPMAT).