Flow-induced demixing of polymer-colloid mixtures in microfluidic channels Arash Nikoubashman, Nathan A. Mahynski, Amir H. Pirayandeh, and Athanassios Z. Panagiotopoulos Citation: The Journal of Chemical Physics 140, 094903 (2014); doi: 10.1063/1.4866762 View online: http://dx.doi.org/10.1063/1.4866762 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/9?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Erratum: Flow-induced demixing of polymer-colloid mixtures in microfluidic channels [J. Chem. Phys.140, 094903 (2014)] J. Chem. Phys. 141, 149906 (2014); 10.1063/1.4897157 Flow-induced migration of polymers in dilute solution Phys. Fluids 18, 031703 (2006); 10.1063/1.2186591 Elastic flow instability, curved streamlines, and mixing in microfluidic flows Phys. Fluids 16, 4028 (2004); 10.1063/1.1792011 Polymer chain dynamics in Newtonian and viscoelastic turbulent channel flows Phys. Fluids 16, 1546 (2004); 10.1063/1.1687415 Design of high-magnetic field gradient sources for controlling magnetically induced flow of ferrofluids in microfluidic systems J. Appl. Phys. 93, 7459 (2003); 10.1063/1.1557361

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

THE JOURNAL OF CHEMICAL PHYSICS 140, 094903 (2014)

Flow-induced demixing of polymer-colloid mixtures in microfluidic channels Arash Nikoubashman, Nathan A. Mahynski, Amir H. Pirayandeh, and Athanassios Z. Panagiotopoulosa) Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, USA

(Received 5 December 2013; accepted 6 February 2014; published online 6 March 2014) We employ extensive computer simulations to study the flow behavior of spherical, nanoscale colloids in a viscoelastic solvent under Poiseuille flow. The systems are confined in a slit-like microfluidic channel, and viscoelasticity is introduced explicitly through the inclusion of polymer chains on the same length scale as the dispersed solute particles. We systematically study the effects of flow strength and polymer concentration, and identify a regime in which the colloids migrate to the centerline of the microchannel, expelling the polymer chains to the sides. This behavior was recently identified in experiments, but a detailed understanding of the underlying physics was lacking. To this end, we provide a detailed analysis of this phenomenon and discuss ways to maximize its effectiveness. The focusing mechanism can be exploited to separate and capture particles at the sub-micrometer scale using simple microfluidic devices, which is a crucial task for many biomedical applications, such as cell counting and genomic mapping. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4866762] I. INTRODUCTION

Colloidal dispersions and polymers in solution are very interesting systems from both a thermodynamic and rheological point of view. In contrast to simple liquids, these fluids respond to external stresses, e.g., shear or pressure gradients, in a highly nonlinear fashion, which gives rise to phenomena including shear-thinning,1, 2 shear-thickening,1 and viscoelasticity.3 The nonlinear response to stress has its origins in the disparate length and time scales between the solvent and suspended particles. Although these non-Newtonian liquids are ubiquitous in nature and technology, e.g., blood4 and crude oil,5 the detailed, microscopic mechanisms responsible for their nonlinear behavior still remain incompletely understood. This is especially true for multicomponent systems such as colloid-polymer mixtures, where in addition to the solute-solvent interactions, the complex interactions between the different species play an important role. To better understand these effects, the dynamics of dilute colloid and polymer solutions have been studied extensively in the last decade with the advent of nano- and microfluidic devices6–8 and sophisticated computer simulations.9, 10 Early studies focused on infinitely dilute suspensions of rigid colloids and flexible polymers in a Newtonian solvent such as toluene or hexane, which is ideal for studying the (flow) properties of single solute particles.11–17 Among the major discoveries was the distinct cross-stream particle migration and resulting highly non-uniform solute distribution for both rigid11, 12, 16 and soft14, 15 particles when subjected to Poiseuille flow. This effect stems from the inhomogeneous a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2014/140(9)/094903/8/$30.00

solvent velocity, which induces a gradient in the stress field across the channel.18 This flow-induced migration has important consequences for many applications where the goal is to precisely control the lateral position of the solute, such as cell counting,19 filtration,20, 21 and genomic mapping.22 In principle, it is possible to overcome this lateral drifting and to spatially focus the solute particles on a streamline by imposing external forces, such as electric fields,23 or by introducing obstacles.24 In practice however, these so-called active methods often require significant modifications of the solute or the channel, and therefore add an unnecessary layer of complexity to the system. Alternatively, one can manipulate the lateral position of the solute particles by inducing nonlinear flows. Such conditions can, for instance, be achieved by tuning the solvent properties. Indeed, early experiments on microscopic particles dispersed in flowing viscoelastic media revealed a particle migration toward the region of minimal shear.25–31 Such a passive approach is a very auspicious candidate for spatial focusing, since it requires only minimal adjustments to the experimental setup. There is a tremendous interest in the biomedical industry of exploiting this mechanism for nanoparticles, due to the possibility of isolating and capturing specific cells and viruses.18, 32, 33 However, it is not immediately clear whether macroscopic viscoelastic effects translate to the nanoscopic colloidal scale, as the effects of Brownian motion become increasingly significant in the sub-micrometer regime. In addition, a continuum description of the surrounding viscoelastic medium becomes problematic at these length scales, since its constituent polymers are comparable in size to the solute particles. Therefore, it is more appropriate to regard these system as binary colloidpolymer mixtures instead. This in turn renders approaches based on continuum fluid mechanics only partially usable,

140, 094903-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

094903-2

Nikoubashman et al.

J. Chem. Phys. 140, 094903 (2014)

and thus one has to resort to alternative methods for modeling these systems, such as molecular dynamics (MD) simulations employed in this work. Indeed, only very recently the general feasibility of viscoelastic focusing has been experimentally demonstrated for nanometric rigid colloids and DNA molecules.34 However, experiments only provide limited insights into the system on a microscopic scale, and therefore the detailed mechanisms responsible for this flow-induced focusing remain largely unclear. In order to better understand this behavior, we studied the flow behavior of rigid spherical colloids dispersed in both Newtonian and viscoelastic solvents, using extensive computer simulations that closely replicate experimental conditions. Our simulations successfully reproduce the viscoelastic focusing onto a centerline seen experimentally at low Reynolds numbers.34 In addition, we found that the focusing cannot be maintained indefinitely, but fails when the flow rate exceeded some threshold value. We provide a detailed microscopic explanation for the pathways involved in these flow-induced transitions and discuss how the parameters can be optimized for viscoelastic focusing in microfluidic devices. The rest of this paper is organized as follows; in Sec. II, we introduce our simulation model and give a brief overview of the employed methodology. Then we present and discuss our results in Sec. III. Finally, we summarize our findings and provide a brief outlook for future problems in Sec. IV.

II. MODEL AND SIMULATION METHOD

Our model system consists of rigid spherical colloids with diameter σ c at a concentration of c = πρc σc3 /6 dispersed in a solvent, where ρ c is the number density of colloids. The solvent in which the colloids are dispersed is constructed on two length scales: an underlying Newtonian solvent described through point-like particles simulated via the Multi-Particle Collision Dynamics (MPCD) simulation scheme,10, 35, 36 and explicitly modeled polymer chains that provide viscoelasticity to the dispersion medium. The viscoelasticity of the solvent was adjusted through the polymer length Np and concentration m = πρp Np σm3 /6, where ρ p is the polymer number density and σ m the monomer diameter. This multiscale approach is motivated by the large time- and length-scale discrepancies between the (microscopic) solvent and the dispersed particles, which make fully atomistic simulations prohibitively time-consuming. Figure 1 shows a simulation snapshot which highlights the essential features of our model. In what follows, we provide a detailed description of the employed model and simulation method. We employed a bead-spring model for the polymers, where each chain is comprised of Np beads with diameter σ m = 1, which defines our unit of length. The bonds between two consecutive monomers are described by the finitely extensible nonlinear elastic (FENE) interaction potential37

UFENE (r) =

  − 12 kr02 ln 1 − ∞,

r2 r02



, r ≤ r0 r > r0

.

(1)

FIG. 1. Simulation snapshot of a colloid-polymer mixture (with concentration c = 0.025, m = 0.050) under Poiseuille flow at Re = 32. Colloids are colored in red, while individual polymers are color-coded. The solvent velocity profile along the channel is schematically shown on the right hand side.

Here, r0 is the maximum bond length and k is the spring constant, for which we have chosen the values r0 = 1.5 and k = 12 ε/σm2 to avoid any unphysical bond crossing even at high flow rates. Excluded volume between the colloids and monomers was modeled through the purely repulsive Weeks-ChandlerAnderson (WCA) potential,38 corresponding to good solvent conditions   σ 12  σ 6  4ε rij + ε, r ≤ 21/6 σij − rij , (2) Uij (r) = 0, r > 21/6 σij where r = |ri − rj | denotes the distance between particles i and j. Furthermore, the parameter ε = 5 kB T controls the strength of the mutual repulsion, while σ ij = (σ i + σ j )/2 is the mean diameter. Our simulations were conducted in a slit-like channel, where the confining walls were separated by a distance of Lx from each other with their normals parallel to the x-axis of the system. Each confining surface is considered to be an immobile monolayer of interaction sites (similar to the monomeric units), with homogeneous density ρσm3 = 1. Since such a detailed representation of the surfaces is computationally expensive and superfluous for the problem at hand, we employed an effective interaction potential obtained by integrating over all wall particles17, 39      √   σ 3 2 2 σij 9 π ε 15 − xij + 310 , x ≤ (2/5)1/6 σij 3 x Uw (x) = . 0, x > (2/5)1/6 σij (3) In Eq. (3) above, x denotes the distance from either wall, and the total external potential caused by both walls is given by Uext (x) = Uw (x) + Uw (Lx − x). We set Lx = 64 (which corresponds to Lx = 8 σ c for σ c = 8), and carefully checked for finite-size and incommensurability effects by conducting additional simulations with Lx = 96 (= 12 σc ) and Lx = 60 (= 7.5 σc ) at selected state points. Periodic boundary conditions were applied in both the y- and z-direction of the system, where we chose Ly = 32 and Lz = 120 to avoid unphysical self-interactions. For such a microfluidic channel, one often refers to the principal axes also as the gradient- (x), vorticity(y), and flow-directions (z). In our simulation scheme, the microscopic solvent particles are considered as ideal point particles with unit mass

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

094903-3

Nikoubashman et al.

J. Chem. Phys. 140, 094903 (2014)

m that do not interact through pairwise potentials with each other and the solute but rather through stochastic momentum exchanges. The MPCD algorithm consists of two alternating steps, a streaming and a collision step. During the former, all solvent particles propagate ballistically over √ a period of t, resulting in a mean free path of λ = t kB T /m. Then, the particles are sorted into cells with edge length a = σ m , and undergo collisions with solvent and solute particles within the same cell. Here, the detailed microscopic solvent collisions are mimicked through a collective rotation of the solvent velocities around a fixed angle α around a randomly chosen unit vector.35 Due to the cell-based nature of this algorithm, the spatial resolution of the hydrodynamic interactions is determined by the mesh size of the collision cells, a. It has been shown in Ref. 40 that Galilean invariance is violated for λ < a/2. Therefore, all lattice cells are shifted by a randomly chosen vector, drawn from a cube with edge length in the interval [−a/2, +a/2] before each collision step. While the aforementioned rules governing the solvent dynamics are general, the simulation of specific flow profiles requires special care and will be discussed briefly further below (see Refs. 10, 17, and 43 for a more detailed description of the simulation algorithm). For the conventional MD-timesteps, we used a Verlet integration schemewith tMD = 2 × 10−3 measured in the unit of time τ = mσm2 /(kB T ). For the MPCD algorithm, we chose a collision angle of α = 130◦ and a timestep of t = 4 × 10−2 , which leads to a mean free path of λ = 0.04. The density of the solvent was set to ρs = 10 σm−3 , resulting in a total number of Ns ≈ 2.5 × 106 solvent particles. These parameters lead to a dynamic viscosity of η = 19.7 τ kB T /σm3 , which was determined by fitting the velocity profile of the pure solvent under Poiseuille flow through Eq. (4). Furthermore, we density-matched the rigid colloids and monomers with the surrounding solvent to achieve neutrally buoyant conditions; each monomer has a mass equal to the average solvent mass per MPCD cell, Mm = 10, while we assigned Mc = ρs π σc3 /6 to each colloid. For experiments conducted at room temperature and channel widths of a few micrometers, these parameters describe a solvent with a viscosity of η ≈ 0.3 cP (comparable to hexane41 ). The systems were first equilibrated for 2.5 × 106 timesteps under quiescent conditions and then for 107 timesteps with the flow turned on. We found that this period was sufficient to reach a steady state for all investigated systems. Then we ran the simulations for additional 107 timesteps to perform the actual measurements. In the case of a pure Newtonian solvent confined in such a microfluidic channel, a parabolic velocity profile emerges between the two confining plates when a pressure drop is applied along the z-direction due to the no-slip boundary conditions at the confining interfaces (see Fig. 1). We emulated this behavior by applying a volume force (similar to gravity) to the solvent particles. The functional form of the velocity profile is given by vz (x) =

g (Lx − x)x, 2ν

(4)

where g is an acceleration constant that can be used to tune the strength of the flow, and ν = η/ρ s denotes the kine-

matic solvent viscosity. The velocity profile vanishes at the channel walls (x = 0 and x = Lx ), and attains its maximum value, vmax , at the center of the channel (x = Lx /2). The mean solvent velocity can be obtained by averaging Eq. (4), yielding v0 = 2vmax /3. The parabolic velocity profile further implies that a local shear rate γ˙ (x) = ∂vz (x)/∂x is established in the fluid, expressed as γ˙ (x) =

g (Lx − 2x). 2ν

(5)

Finally, we can then derive from Eq. (5) the average shear rate between the centerline and either channel wall as γ˙  = gLx /(4ν). In MPCD, no-slip boundary conditions can be achieved for planar walls coinciding with the collision cell boundaries by employing a bounce-back rule for the solvent particles. However, for a more general system geometry the walls may not coincide with, or even be parallel to the collision cell boundaries. Furthermore, partially occupied collision cells can also emerge from the cell-shifting, which is unavoidable for the small mean free paths λ employed in our simulations. Lamura et al.42 have demonstrated that in such a case the bounce-back rule has to be modified by refilling the boundary cells with virtual particles, and we employed this algorithm in our simulations. Additionally, thermostatting is required in any non-equilibrium MPCD simulation whenever either isothermal conditions are required or viscous heating can occur. Therefore, we rescaled the solvent velocities at the cellular level to maintain constant temperature.43 Finally, special care has to be taken when setting the flow strength to avoid any significant compressibility effects due to the ideal nature of the solvent particles. Hence, the maximum fluid velocity vmax should √ be well below the speed of sound in the medium cs = 5kB T /(3m). We found that the Mach number, Ma = vmax /cs , was less than 0.15 for the highest flow strengths employed in our simulations. This is well below the threshold Ma = 0.25 above which compressibility effects become significant,41 and therefore we can consider our solvent to be an incompressible liquid. In order to quantify the strength of the flow, we introduce the Reynolds number Re, which is a measure of the ratio of the inertial forces to the viscous forces, and define it as Re =

L3 Lx v0 = g x2 . ν 12ν

(6)

In the limit of Re 1, i.e., Stokes flow, neutrally buoyant solute particles are confined on a streamline and do not exhibit any cross-stream migration, because of the linearity and thus time reversibility of the system.44 However, for stronger flows (1 ≤ Re ≤ 100) the reciprocity of the system is broken, and lateral lift forces lead to a distinct cross-stream migration away from the channel’s centerline.11, 12, 15, 16 The strength of these lift forces is highly dependent on the size and the position of the solute particles in the microchannel, and expressions have been determined analytically by asymptotic expansion for small solute particles,45 by simulations,16, 46 and microfluidic experiments.46 These inertially driven lateral forces are however suppressed for solute particles dispersed in viscoelastic

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

094903-4

Nikoubashman et al.

media;25–31 these systems exhibit an inhomogeneous stress under flow due to the polymeric nature of the solvent, which results in a solute migration toward the region of lowest γ˙ . For steady Poiseuille flow, this region corresponds to the centerline, where the local shear-rate vanishes completely. Expressions for the viscoelastic lift force have been derived for dispersed particles with σ c 1 μm,28–30 but these findings are not readily applicable to nanoparticles, since Brownian motion plays a significant role at these length scales.30 Another complication arises from the fact that nanometric solute particles have a comparable size as the polymers constituting the viscoelastic solvent. Hence, it is more appropriate to regard the system as a colloid-polymer mixture than as colloids embedded in a continuous viscoelastic medium. An important measure for polymer solutions is the Weissenberg number Wi = τp γ˙ , which is the product of polymer relaxation time τ p and (local) shear-rate. It has been shown in Ref. 47 that the MPCD algorithm correctly reproduces Zimm dynamics for sufficiently small solvent mean free paths λ  0.1. Since we have λ = 0.04 in our simulations, we can safely assume that our polymers conform with the Zimm model, and therefore use the following expression for τ p :48 τp,n =

3 ηRg,0

. (7) kB T n3γ Here, τ p,n is the relaxation time of the nth mode, and γ = 0.588 is the Flory exponent for good solvent conditions. For steady shear, the dynamics of the polymer are dictated by its longest relaxation time, and therefore we can approximate τ p ≈ τ p,1 . III. RESULTS

In order to elucidate the origin and properties of nanoscale viscoelastic focusing for colloid-polymer mixtures in microfluidic channels under Poiseuille flow, we systematically varied the colloid and polymer concentrations (0 ≤ c ≤ 0.075 and 0 ≤ m ≤ 0.075) and the flow strength (0 ≤ Re ≤ 128). In addition, we also studied the influence of the colloid diameter (σ c = 6, 8, 10) and polymer length (Np = 5, 10, 20, and 40). If not otherwise stated, we employed σ c = 8 and Np = 40, where the polymers had an equilibrium radius of gyration 2Rg,0 ≈ σ c . For the highest investigated concentrations with colloids of diameter σ c = 8 and chains of length Np = 40, these values correspond to 69 colloids and 880 polymers. We first studied the behavior of pure systems of colloids or polymers, i.e., m = 0 or c = 0, in order to verify our algorithm and establish a reference for the mixed systems. Figure 2 schematically shows the particle distributions along the gradient direction of the flow as a function of Re and concentrations. In the pure colloid case, the solute particles were evenly distributed in the channel at zero flow [see Fig. 2(a)]. When the flow was turned on, colloids migrated toward the channel walls and formed hexagonally packed layers in the flow-vorticity plane close to the surfaces due to the shearinduced inertial lift forces. As Re was increased, the colloids were pushed further to the walls and formed a more tightly packed hexagonal structure at the surfaces. In the case of a

J. Chem. Phys. 140, 094903 (2014)

FIG. 2. Schematic representation of the density distributions ρ(x) of colloids (continuous, red) and polymers (dashed, blue) along the gradient direction (x) of the flow. (a), (d), and (g) correspond to the pure colloid case (m = 0), (b), (e), and (h) correspond to the composite system, and (c), (f), and (i) correspond to the pure polymer case (c = 0). (a)–(c) show ρ(x) at quiescent conditions, while (d)–(f) and (g)–(i) show the system at Reynolds numbers below and above the critical value Rec , respectively.

pure polymer solution, we observed a homogeneous density distribution ρ p (x) along the x-direction under quiescent conditions [see Fig. 2(c)]. As we slowly increased the flow strength, polymers migrated away from the centerline to the walls as shown in panels (f) and (i) of Fig. 2. This cross-stream migration is caused by the inhomogeneity of γ˙ (x), which leads to a position-dependent deformation and shear-alignment of the polymers. This in turn leads to a gradient in the chain mobility which manifests itself in the partial evacuation of the centerline as Re is increased. These results confirm previous studies of colloidal dispersions11, 12, 16 and polymer solutions13–15 under flow in the limit of infinite dilution, and we can therefore conclude that the characteristic features at  → 0 are indeed qualitatively preserved for the significantly higher  investigated here. However, upon mixing the colloids and polymers, the resulting density profiles are not superpositions of the individual ρ(x), indicating a strong interaction between the two species. Indeed, it is already well established through experiments49 and theory50, 51 that, under quiescent conditions, athermal colloid-polymer mixtures phase separate due to depletion effects at sufficiently high concentrations [see Fig. 2(b)]. This purely entropic phenomenon arises as colloids approach one another, excluding polymers from existing in the gap between them, resulting in an osmotic pressure difference between the bulk and locally depleted region which drives the colloids to aggregate. The phase behavior of this binary mixture strongly depends on the size ratio between the two components, generally classified in terms of q = 2Rg,0 /σ c . The so-called “colloid” limit, where q  1, has been more extensively studied than its counterpart the “protein” limit (q 1), since the polymer’s internal degrees of freedom play a less consequential role in the overall free energy landscape at densities where

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

Nikoubashman et al.

0 0

0.25

0.5 x / Lx

0.75

1

FIG. 3. Colloid density ρ c along the gradient direction for a colloid-polymer mixture with c = 0.025 and m = 0.050.

phase separation typically occurs.51, 52 We chose to focus our efforts in this regime as well, since it is experimentally the more common of the two. As we turned on the flow, a significant fraction of colloids migrated to the centerline of the microchannel, thereby expelling the polymers toward the confining walls. This configuration is schematically shown in Fig. 2(e), while Fig. 3 shows the density profiles from simulations. Here, we can clearly identify a colloid- and polymer-rich region with a very sharp interface between them. This flow-induced demixing occurred for all investigated polymer lengths and particle concentrations, with some quantitative differences which will be discussed further below. As we increased the flow strength above some critical Rec , we measured a steady evacuation of the centerline with a concomitant increase of colloids in the vicinity of the channel walls, until the channel center was entirely devoid of colloids [cf. Fig. 2(h)]. In order to understand the origin of the flow-induced colloid focusing, let us consider a system where the colloids are initially placed in the centerline of the channel [cf. Fig. 2(e)]; in the case of a Newtonian solvent, inertial lift forces would immediately push the solute particles toward the channel walls upon turning on the flow.11, 12, 15, 16 However, the addition of polymers can reinforce this otherwise unstable stream of colloids propagating in the centerline, since they act as a buffer impeding the cross-stream migration that would otherwise draw the colloids to the edges of the channel. The rigid colloids do not suffer from any intrinsic penalty by existing in the centerline owing to their lack of internal degrees of freedom. Instead, they are subjected solely to viscous drag due to the no-slip boundary conditions at their surfaces. This friction can be reduced for colloids that propagate in the wake of others, promoting the emergence of well-ordered layers. Concomitantly, the polymers are exposed to higher local shear rates by migrating closer to the walls, resulting in a stronger alignment in flow direction. This in turn leads to an overall reduction of the intermolecular friction, since these off-center planes of aligned polymers simply slip past one another as opposed to tumbling in the channel center (see Fig. 8). However, this balance cannot be sustained for arbitrarily high flow rates as viscous forces become progressively less significant as the dynamics are dominated by inertial effects instead. Hence, the inertial lift forces experienced by the colloids become stronger with increasing flow strength, eventually allowing them to break through the polymer layers and migrate to the channel edges. In order to quantify this be-

1

Np = 40 Np = 20 Np = 10 Np = 5

0.8 / c

2

0.6

ctr c

4

havior, we measured the fraction of colloids in the centerline, ctr c , as a function of Re. Here, we assign a colloid i to the centerline when c (xi ) > 0.1 c (L/2) at its lateral position xi . Although we do not a priori impose that this be symmetric, we find that this occurs naturally to within a few per cent. Figure 4(a) shows ctr c for all investigated values of Np , and it is clearly visible that the range of Re over which colloids can be trapped in the centerline increases with Np . The reason for this is twofold; first of all, the shear-alignment of long polymers requires considerably lower γ˙ compared to short ones, and therefore the stabilizing polymer buffer emerges at lower Re. Second, it is more difficult for a colloid to migrate through a layer of long polymers as opposed to short ones, since a layer of short shear-aligned polymers contains less bonds on average than a layer of long polymers with the equivalent monomer density. Thus, for short chains there exist more gaps between them through which a colloid can migrate while causing less significant perturbations to the

0.4 0.2 (a) 0 16

32

48

64

80

96

112

128

Re 1

c= c= c= c=

0.8 / c

Re = 16 Re = 48 Re = 112

0.075, 0.075, 0.075, 0.025,

m= m= m= m=

0.025 0.050 0.075 0.050

0.6

ctr c

c(x) c

3

6

J. Chem. Phys. 140, 094903 (2014)

0.4 0.2 (b) 0 16

32

48

64

80

96

112

128

Re 60 Rec

094903-5

40 20

Simulation data Logarithmic curve fit

(c) 0 0

200

400

600

800

1000

1200

1400

p

FIG. 4. Panel (a): Fraction of colloids located in the channel center ctr c vs. Re for all investigated values of Np at fixed c = 0.075 and m = 0.050. The symbols show our simulation data, while the continuous lines are fits to a modified error function. Panel (b): Same as (a) but for all investigated values of c and m at fixed Np = 40. Panel (c): Logarithmic curve fit of Rec as a function of τ p for c = 0.075 and m = 0.050.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

Nikoubashman et al.

J. Chem. Phys. 140, 094903 (2014)

vL ∝

τp σc2 ∇ γ˙ (x)2 .

(9)

Equation (9) neglects the effect of the colloid diameter σ c , which could influence the efficiency of focusing as well. On the one hand, Eq. (8) suggests viscous lift forces ∝ σc3 toward the channel center in the regime Re 1. On the other hand, it is well known from experiments46 and simulations16 of colloids in Newtonian solvents that the inertial lift forces which push the colloids away from the centerline are proportional to σc3 in the regime Re 1. Hence, these two effects are competing against each other at intermediate Re, and need to be balanced in order to achieve an optimum focusing efficiency. To this end, we conducted additional simulations for σ c = 6 and σ c = 10. Figure 5 shows ctr c vs. Re, and it is clear

c= 6 c= 8 c = 10

/ c

0.8 0.6 0.4

(8)

This lift velocity vL has been derived in the limit of Re 1 and Pe 1, where inertial forces and thermal fluctuations can be neglected entirely. Equation (8) suggests that all colloids will eventually migrate to the center and remain there. This is obviously not the case for nanoscale particles, since Brownian motion can override vL and therefore drive particles out of the centerline. Figure 4(b) shows ctr c for the investigated concentrations c and m , and it is readily apparent that increasing the polymer concentration leads to a significant increase of both ctr c and Rec . Clearly, the overall reduction of viscous drag achieved by replacing polymers with colloids in the centerline increases with m , until the centerline is entirely packed. Moreover, it also becomes harder for colloids to drift away from the centerline due to the denser polymer buffer layer. This interpretation of the data is corroborated by the simulation results with lower c = 0.025, where ctr c decays slightly faster compared to the runs at c = 0.075, which might seem counterintuitive at first glance. In this regard, colloids provide a considerable amount of very mobile excluded volume which allows the polymer to evacuate unfavorable channel regions, whereupon taking the polymer’s place a colloid incurs no additional hydrodynamic penalties owing to its lack of internal degrees of freedom. From these measurements, we can suggest an empirical expression for Rec depending on the polymer relaxation time τ p [via Np and Eq. (7)] and polymer concentration m . To this end, we define Rec as the position of the inflection 55 point in ctr c (Re), and determine its value by fitting our data through a modified error function. Following our previous arguments, we can identify the two limiting cases as Rec → 0 for Np → 0 and Rec → ∞ for Np → ∞. The expression ln (1 + τ p ) satisfies these two conditions, and fits our data well [see Fig. 4(c)]. Furthermore, we found that Rec ∝ m for the m values investigated here. Hence, we posit the following empirical relationship for Rec : Rec ∝ m ln(1 + τp ).

1

ctr c

otherwise flow-aligned polymers. Such perturbations would incur a significant drag penalty, and impede such motion. These findings are consistent with results from microfluidic experiments on micrometer-sized colloids30, 31 and continuum fluid dynamics calculations,53, 54 where an accelerated particle motion toward the centerline was observed with increasing τ p and γ˙

0.2 0 16

32

48

64

80

96

112

128

Re FIG. 5. Fraction of colloids located in the channel center ctr c vs. Re for all investigated values of σ c at c = 0.075, m = 0.050, and Np = 40.

that ctr c decays slower for larger solute particles. We also observe a slight linear increase of Rec . This discrepancy stems from the fact that the predictions for Re 1 are neglecting viscoelastic effects, which still play a significant role even at high flow rates. The microscopic nature of our simulations allows us to study the colloid trajectories and polymer conformations in detail, which provide valuable additional insights into the origin of the flow-induced demixing. In order to better understand the breakup mechanism at Re > Rec , we conducted additional simulations in which we initialized the system from a focused state [cf. Fig. 2(e)], and tracked the lateral position of individual colloids as a function of time. Figure 6 shows the trajectory of one such representative particle at intermediate and high Reynolds numbers. For Re = 48 < Rec , we can see that the colloid remains on its initial central position over the whole course of the simulation. As we increased Re beyond its critical value, the particle abruptly left its centerline position and migrated toward the channel walls until it reached a distance of Lx /4, and remained there for ≈4000 τ . Subsequently, the colloid continued its outward motion until it eventually reached its steady state position close to the wall. A similar transition can be observed for the highest flow 1

Re=48 Re=96 Re=128

0.75 x / Lx

094903-6

0.5 0.25 0 0

5000

10000

15000

20000

25000

t/ FIG. 6. Lateral colloid position x as a function of time t for c = 0.075 and m = 0.050. The system has been initialized from a focused configuration. Gray solid lines indicate layers where colloids aggregate and spend most of their time. Separating these are layers of high polymer density as drawn schematically in Fig. 2, which act as barriers impeding the movement of the colloids toward the channel walls. These are progressively less significant as Re increases.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

Nikoubashman et al.

094903-7

J. Chem. Phys. 140, 094903 (2014)

strengths investigated here (Re = 128). These data clearly show the colloid propagating through well-defined layers of polymers. At low Re, the colloids remain focused in the centerline, however, upon increasing the flow strength they migrate through the layers with increasing ease. Next, we studied the polymer conformation at different colloid concentrations and flow strengths. Figure 7(a) shows the spatial distribution of the diagonal components of the gyαα , ration tensor normalized by its value at rest, Rgαα (x)/Rg,0 between the confining walls at Re = 16 for a polymer with Np = 40 beads. We also measured the bond stretching due to the applied flow, and found that the stretching was below 1% for the highest flow strengths applied. It is clear that the polymer elongates in the flow direction and accordingly contracts in the remaining ones, though the reduction is always more severe in the gradient direction owing to confinement. However, this trend is significantly reduced very near the center of the channel due to the vanishing shear rate there; instead of forming a shear-aligned filament, the polymer follows the parabolic velocity profile of the solvent and resembles a “hairpin.” We observed the same qualitative behavior for all investigated values of Np . Furthermore, it is clear that these trends are independent of colloid concentration, which conforms with previous simulation results for single polymers in solutions of colloids under shear flow.56 Figure 7(b) shows the normalized Rgαα as a function of Re and Wi, where each data point has been obtained by 4

c =0.0, c =0.025, c =0.075,

(a)

3.5

2.5 2 1.5 1 0.5 0 0

0.25

0 3.5

35

0.5 x / Lx

0.75

105 140

70

1

175

2Rgxz mg = tan(2χ ) = zz . Wi Rg − Rgxx

210 0.25

(b)

10

0.2 c=0.0, c=0.025, c=0.075,

2.5 2

m =0.050 m =0.050 m =0.050

0.15

1

0.1

1.5

/

(Re) / Rg,0

3

Rg

(10)

Figure 8 shows χ for the cases c = 0 and c = 0.075 at Re = 96 as a function of x, and it is clear that, except near the centerline, the polymers are nearly perfectly aligned with the flow. Far from the centerline, complete shear-alignment already took place at very low Re ≈ 16 and no significant improvement was achieved upon increasing the flow strength (not shown here). However, in the center of the channel, where γ˙ = 0, polymer alignment is reduced dramatically and tumbling leads to large variances in χ . This effect is amplified as Re → 0 whereupon lower shear rates weaken the coupling between the flow and polymer orientation. Tumbling polymers collide with their shear-aligned neighbors and thus increase their intermolecular friction; while this is unfavorable, it is inevitable given the polymer density in the channel. This trend is quantified in the inset of Fig. 8, where we plotted mg as a function of Wi. It is well known from theory48 and experiments57 that close to equilibrium, i.e., Wi → 0, the deformation of a sheared linear polymer is given by Rgxz ∝ γ˙ and (Rgzz − Rgxx ) ∝ γ˙ 2 . Hence, mG is independent of Wi in this regime. For larger shear rates, the orientational resistance

mG

Rg

(x) / Rg,0

3

m=0.050 m=0.050 m=0.050

averaging over the channel width Lx . It is clearly visible that the polymer chains are already almost fully aligned in the flow direction at very low Re, and that the alignment quality does not improve significantly for higher Re. This trend becomes clear, considering that the corresponding Wi 1. Furthermore, Fig. 7(b) reveals that the average polymer extension in the flow direction slightly increases for higher colloid concentrations c at intermediate Reynolds numbers Re < Rec . This is due to the fact that more colloids occupy the centerline and thereby push more polymers toward the channel edges, where the local shear rate is higher. This effect vanishes for Re > Rec , and Rgαα becomes independent of m . We then quantified the average alignment of the polymers with respect to the flow by measuring orientational resistance mg , based on the angle χ , between the eigenvector of the gyration tensor with the largest eigenvalue and the flow direction.57, 58 For the system at hand, geometric considerations yield

0.05

0.1

1

10 100

c=0.0, c=0.075,

m =0.050 m =0.050

0

1

-0.05

0.5

-0.1

0 0

16

32

48

64 Re

80

96

112

128

FIG. 7. Panel (a): Diagonal components of the polymer gyration tensor Rgαα at Re = 16, normalized by their values at rest, along the gradient direction for a pure polymer system (lines) and mixtures (symbols). Panel (b): Rgαα averaged over the channel width Lx as a function of Re and Wi. In both panels, from top to bottom: flow (z), vorticity (y), and gradient direction (x).

-0.15 0

0.25

0.5 x / Lx

0.75

1

FIG. 8. Main plot: orientational resistance, mg , of polymers with Np in the absence (red) and presence (blue) of colloids at Re = 96 (main). Inset: mG from simulations as a function of Wi (symbols) together with its power law fit (dashed line).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

094903-8

Nikoubashman et al.

can be described by a power law mG (Wi) ∝ Wiμ , where a value of μ = 0.54 ± 0.03 has been reported in Ref. 59 for ultradilute polymer solutions in Couette flow. We have measured a somewhat higher exponent μ = 0.76 ± 0.02, which is likely due to the higher polymer concentration and inhomogeneous shear field in our simulations. From these data, we can safely conclude that the presence of the spherical colloids does not have any significant impact on the shape of the polymers for the parameters investigated. IV. CONCLUSIONS

We performed extensive molecular dynamics simulations, coupled with Multi-Particle Collision Dynamics to take into account solvent-mediated hydrodynamics, to study the viscoelastic focusing of nanoparticles in microfluidic devices. We explicitly modeled the constituent polymers of the viscoelastic medium and provided a detailed microscopic explanation of the phenomenon. We found that the emergence of the steady stream of colloids in the centerline at low and intermediate Reynolds numbers Re is due to the cooperative minimization of viscous drag. Commensurately, as Re exceeded some critical value Rec , inertial influences overrode this effect leading to an abrupt decay in the centerline density peak. The specific value of Rec strongly depends on the properties of the viscoelastic medium, and we found that the colloid focusing could be maintained over a broader Re range for longer polymers and higher polymer concentrations. Although our investigation have provided a clearer picture of the mechanisms behind viscoelastic focusing of nanoparticles, there still remain some open questions, such as the influence of enthalpic polymer-colloid interactions. Furthermore, it is interesting to investigate the viscoelastic focusing of flexible and rigid polymers instead of rigid colloids, since such scenarios are highly relevant for many biomedical applications. ACKNOWLEDGMENTS

A.N. acknowledges financial support from the Princeton Center for Complex Materials (PCCM), a U.S. National Science Foundation (NSF) Materials Research Science and Engineering Center (Grant No. DMR-0819860). In addition, N.A.M. received support from the U.S. National Science Foundation (Grant No. CBET-1033155). 1 N.

J. Wagner and J. F. Brady, Phys. Today 62(10), 27 (2009). M. Erwin, D. Vlassopoulos, and M. Clotre, J. Rheol. 54, 915 (2010). 3 R. B. Bird, R. C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids (Wiley Interscience, New York, 1987). 4 B. M. Johnston, P. R. Johnston, S. Corney, and D. Kilpatrick, J. Biomech. 39, 1116 (2006). 5 E. S. Boek, H. K. Ladva, J. P. Crawshaw, and J. T. Padding, Energy Fuels 22, 805 (2008). 6 G. M. Whitesides and A. D. Stroock, Phys. Today 54(6), 42 (2001). 7 T. M. Squires and S. R. Quake, Rev. Mod. Phys. 77, 977 (2005). 8 S. Buzzaccaro, E. Secchi, and R. Piazza, Phys. Rev. Lett. 111, 048101 (2013). 9 S. Chen and G. D. Doolen, Annu. Rev. Fluid Mech. 30, 329 (1998). 10 G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, Adv. Polym. Sci. 221, 1 (2009). 11 G. Segre and A. Silberberg, Nature 189, 209 (1961). 2 B.

J. Chem. Phys. 140, 094903 (2014) 12 J.-P.

Matas, J. Morris, and É. Guazzelli, J. Fluid. Mech. 515, 171 (2004). Khare, M. D. Graham, and J. J. de Pablo, Phys. Rev. Lett. 96, 224505 (2006). 14 R. Chelakkot, R. G. Winkler, and G. Gompper, Europhys. Lett. 91, 14001 (2010). 15 S. Reddig and H. Stark, J. Chem. Phys. 135, 165101 (2011). 16 C. Prohm, M. Gierlak, and H. Stark, Eur. Phys. J. E 35, 80 (2012). 17 A. Nikoubashman, C. N. Likos, and G. Kahl, Soft Matter 9, 2603 (2013). 18 A. Karimi, S. Yazdi, and A. M. Ardekani, Biomicrofluidics 7, 021501 (2013). 19 S. C. Hur, H. T. K. Tse, and D. D. Carlo, Lab Chip 10, 274 (2010). 20 A. A. S. Bhagat, S. S. Kuntaegowdanahalli, and I. Papautsky, Phys. Fluids 20, 101702 (2008). 21 Y. S. Choi, K. W. Seo, and S. J. Lee, Lab Chip 11, 460 (2011). 22 K. Jo, Y.-L. Chen, J. J. de Pablo, and D. C. Schwartz, Lab Chip 9, 2348 (2009). 23 N. Narayanan, A. Saldanha, and B. K. Gale, Lab Chip 6, 105 (2006). 24 L. R. Huang, E. C. Cox, R. H. Austin, and J. C. Sturm, Science 304, 987 (2004). 25 K. Karnis and S. G. Mason, Trans. Soc. Rheol. 10, 571 (1966). 26 F. Gauthier, H. L. Goldsmith, and S. G. Mason, Rheol. Acta 10, 344 (1971). 27 F. Gauthier, H. L. Goldsmith, and S. G. Mason, Trans. Soc. Rheol. 15, 297 (1971). 28 B. P. Ho and L. G. Leal, J. Fluid. Mech. 76, 783 (1976). 29 M. A. Tehrani, J. Rheol. 40, 1057 (1996). 30 A. M. Leshansky, A. Bransky, N. Korin, and U. Dinnar, Phys. Rev. Lett. 98, 234501 (2007). 31 S. Yang, S. S. Lee, S. W. Ahn, K. Kang, W. Shim, G. Lee, K. Hyun, and J. M. Kim, Soft Matter 8, 5011 (2012). 32 D. C. Colter, I. Sekiya, and D. J. Prockop, Proc. Natl. Acad. Sci. U.S.A. 98, 7841 (2001). 33 X. Cheng, D. Irimia, M. Dixon, K. Sekine, U. Demirci, K. Zamir, R. G. Tompkins, W. Rodriguez, and M. Toner, Lab Chip 7, 170 (2007). 34 J. Y. Kim, S. W. Ahn, S. S. Lee, and J. M. Kim, Lab Chip 12, 2807 (2012). 35 A. Malevanets and R. Kapral, J. Chem. Phys. 110, 8605 (1999). 36 A. Malevanets and R. Kapral, J. Chem. Phys. 112, 7260 (2000). 37 M. Bishop, M. H. Kalos, and H. L. Frisch, J. Chem. Phys. 70, 1299 (1979). 38 J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971). 39 D. A. Lenz, R. Blaak, and C. N. Likos, Soft Matter 5, 2905 (2009). 40 T. Ihle and D. M. Kroll, Phys. Rev. E 63, 020201 (2001). 41 E. F. Cooper and A. F. A. Asfour, J. Chem. Eng. Data 36, 285 (1991). 42 A. Lamura, G. Gompper, T. Ihle, and D. M. Kroll, Europhys. Lett. 56, 319 (2001). 43 A. Nikoubashman and C. N. Likos, J. Chem. Phys. 133, 074901 (2010). 44 L. G. Leal, Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes (Cambridge University Press, Cambridge, 2007). 45 E. S. Asmolov, J. Fluid. Mech. 381, 63 (1999). 46 D. Di Carlo, J. F. Edd, K. J. Humphry, H. A. Stone, and M. Toner, Phys. Rev. Lett. 102, 094503 (2009). 47 M. Ripoll, K. Mussawisade, and G. Gompper, Europhys. Lett. 68, 106 (2004). 48 M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, 1999). 49 R. Tuinier, P. A. Smith, W. C. K. Poon, S. U. Egelhaaf, D. G. A. L. Aarts, H. N. W. Lekkerkerker, and G. J. Fleer, Europhys. Lett. 82, 68002 (2008). 50 G. J. Fleer and R. Tuinier, Phys. Rev. E 76, 041802 (2007). 51 N. A. Mahynski, T. Lafitte, and A. Z. Panagiotopoulos, Phys. Rev. E 85, 051402 (2012). 52 V. J. Anderson and H. N. W. Lekkerkerker, Nature (London) 416, 811 (2002). 53 P. Y. Huang, J. Feng, H. H. Hu, and D. D. Joseph, J. Fluid Mech. 343, 73 (1997). 54 P. Y. Huang and D. D. Joseph, J. Non-Newtonian Fluid Mech. 90, 159 (2000). 55 A. Nikoubashman, G. Kahl, and C. N. Likos, Phys. Rev. Lett. 107, 068302 (2011). 56 H. Chen and A. Alexander-Katz, Phys. Rev. Lett. 107, 128301 (2011). 57 A. Link and J. Springen, Macromolecules 26, 464 (1993). 58 M. Ripoll, R. G. Winkler, and G. Gompper, Phys. Rev. Lett. 96, 188302 (2006). 59 R. G. Teixeira, H. P. Babcock, E. S. G. Shaqfeh, and S. Chu, Macromolecules 38, 581 (2005). 13 R.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.70.241.163 On: Sat, 20 Dec 2014 07:40:03

Flow-induced demixing of polymer-colloid mixtures in microfluidic channels.

We employ extensive computer simulations to study the flow behavior of spherical, nanoscale colloids in a viscoelastic solvent under Poiseuille flow. ...
935KB Sizes 0 Downloads 3 Views