Cross-stream migration of a Brownian droplet in a polymer solution under Poiseuille flow

TL;DR

This study uses simulations to explore how polymer solutions influence the lateral migration of Brownian droplets in microchannels under flow, revealing polymer-induced focusing effects and shape dynamics.

Contribution

It demonstrates that adding polymers causes droplets to focus at the channel centerline, contrasting with inertial migration in Newtonian fluids, and explores the effects of polymer length and concentration.

Findings

Droplets migrate toward walls in Newtonian fluids due to inertia.

Polymers cause droplets to focus at the channel centerline.

Flow-rate and polymer properties influence droplet shape and position.

Abstract

The migration of a Brownian fluid droplet in a parallel-plate microchannel was investigated using dissipative particle dynamics computer simulations. In a Newtonian solvent, the droplet migrated toward the channel walls due to inertial effects at the studied flow conditions, in agreement with theoretical predictions and recent simulations. However, the droplet focused onto the channel centerline when polymer chains were added to the solvent. Focusing was typically enhanced for longer polymers and higher polymer concentrations with a nontrivial flow-rate dependence due to droplet and polymer deformability. Brownian motion caused the droplet position to fluctuate with a distribution that primarily depended on the balance between inertial lift forces pushing the droplet outward and elastic forces from the polymers driving it inward. The droplet shape was controlled by the local shear rate, and so its average shape depended on the droplet distribution.

Equations

$${\mathbf{F}}_{\mathrm{C}}=\{\begin{array}{ll}{a}_{ij}\left(1-r/r_{\mathrm{c}}\right)\hat{\mathbf{r}}& r\le r_{\mathrm{c}}\\ 0& r>r_{\mathrm{c}}\end{array},$$

$${\mathbf{F}}_{\mathrm{D}}$$

$${\mathbf{F}}_{\mathrm{R}}$$

$$w(r)=\{\begin{array}{ll}(1-r/r_{\mathrm{c}}{)}^s& r\le r_{\mathrm{c}}\\ 0& r>r_{\mathrm{c}}\end{array},$$

$$\sigma =\frac{L_z}{2}⟨p_{zz}-\frac{p_{xx}+p_{yy}}{2}⟩,$$

$$\sigma =0.75\rho k_{\mathrm{B}}T{r}_{\mathrm{c}}{\chi }^{0.26}(1-2.36/\chi {)}^{3/2},$$

$${\mu }_{\mathrm{s}}=\frac{315k_{\mathrm{B}}T}{128\pi {\gamma }_{\mathrm{s}}{r}_{\mathrm{c}}^3}+\frac{512\pi {\gamma }_{\mathrm{s}}{\rho }^2{r}_{\mathrm{c}}^5}{51975},$$

$${\mathbf{F}}_{\mathrm{B}}=-\kappa \left(r-b\right)\hat{\mathbf{r}},$$

$$u_x(z)=U\left[1-{\left(\frac{z}{H}\right)}^2\right],$$

$$u_x(z)=\frac{n}{n+1}{\left(\frac{\rho f_x{H}^{n+1}}{K}\right)}^{1/n}\left[1-{\left(\frac{|z|}{H}\right)}^{1+1/n}\right].$$

$$G_{\alpha \beta }=\frac{1}{N_{\mathrm{d}}}\sum _{i=1}^{N_{\mathrm{d}}}\Delta r_{i,\alpha }\Delta r_{i,\beta },$$

$$D\approx \frac{{\lambda }_1-({\lambda }_2+{\lambda }_3)/2}{{\lambda }_1+({\lambda }_2+{\lambda }_3)/2}.$$

$$\tan (2\theta )=\frac{2G_{xz}}{G_{xx}-G_{zz}}.$$

Full text

I Introduction

Particle migration in a microchannelStone and Kim (2001); Stone, Stroock, and Ajdari (2004) has important applications in separation technologies such as filtration Bhagat, Kuntaegowdanahalli, and Papautsky (2008), cell sorting Hur et al. (2011), and fractionation Giddings (1993). It also has implications for physical processes like the margination of cells in the blood streamKumar and Graham (2012a, b) and for multiphase flows in geological formations (enhanced oil recovery) Green and Willhite (1998); Wilson Jr. (1977); Tehrani (1996). Such cross-stream migration could be desirable if a separation is needed but undesirable if a homogeneous distribution is preferred, and it is important to understand and design the conditions under which migration occurs. Multiple mechanisms exist for cross-stream migration in microchannels Fu et al. (1999); Cohen (2005); Lee, Purdom, and Westervelt (2004); Das, Mandal, and Chakraborty (2018); Das and Chakraborty (2018), but in this article we will focus on particle migration that is passively controlled by a pressure- or gravity-driven flow Karimi, Yazdi, and Ardekani (2013); Amini, Lee, and Di Carlo (2014), which is attractive from an engineering perspective for its potential as a scalable, high-throughput technology.

Rigid particles in a Newtonian fluid are known to move across streamlines in parabolic (Poiseuille) flows due to lift forces at small but finite fluid inertia Leal (1980); Stone (2000). Inertial lift outward from the channel center is balanced by an inward force induced by hydrodynamic interactions with the walls, causing the particle to adopt an intermediate lateral position Matas, Morris, and Guazzelli (2004); Di Carlo et al. (2007); Amini, Lee, and Di Carlo (2014). This effect was first observed experimentally by Segré and Silberberg Segré and Silberberg (1961), who found that millimeter-sized spheres in pipe flow migrated to an annulus at roughly 60% of the pipe radius. The number and position of these “focusing” points depends on the channel geometry and flow, and has also been demonstrated for, e.g., parallel platesHo and (1974) and square ducts Di Carlo et al. (2009).

Deformable droplets in a Newtonian fluid exhibit an even richer set of behaviors than their rigid counterparts Leal (1980). Unlike rigid spheres, droplets can migrate across streamlines even in the Stokes flow (inertialess) limit due to their deformability. Chan and Leal showed that the direction of this migration depends on the viscosity ratio between the droplet and the fluid Chan and Leal (1979). Stan et al. found that chemical and surfactant-induced Marangoni effects also influenced droplet migration Stan et al. (2011, 2013). At finite fluid inertia, Legendre and Magnaudet demonstrated that there is lift on a droplet Legendre and Magnaudet (1997) analogous to the Saffman lift on a rigid particle Saffman (1965, 1968) but with a magnitude that depends on the viscosity ratio between the droplet and the fluid. Experiments Karnis, Goldsmith, and Mason (1966); Hur et al. (2011) and simulations Mortazavi and Tryggvason (2000); Chen et al. (2014); Pan et al. (2016); Marson et al. (2018) have shown that droplets undergo Segré–Silberberg-type migration in Poiseuille flow, and that the preferred lateral position depends on several dimensionless parameters, as recently discussed in detail by Marson et al. Marson et al. (2018).

High-throughput applications like filtration or sorting may require focusing particles onto the channel centerline Di Carlo et al. (2007); Karimi, Yazdi, and Ardekani (2013), which is not always achieved by inertial or deformation-induced migration in simple channel geometries. Considerable efforts have been dedicated to design various microfluidic device geometries that can manipulate particles in this way Di Carlo et al. (2007), but finding such geometries can be difficult and highly problem specific Amini, Lee, and Di Carlo (2014). Fortunately, it has been shown that the addition of polymers to the Newtonian solvent provides a simple mechanism, called viscoelastic focusing D’Avino, Greco, and Maffettone (2017), to drive particles toward regions of low shear.

Viscoelastic polymer solutions induce inward particle migration in Poiseuille flow due to a gradient in the first normal stress difference over the particle surface Leshansky et al. (2007). The elastic force exerted by the polymers competes directly with other forces acting on the droplet for the flow conditions, including inertial lift, deformation-induced forces, and wall forces, to set the lateral position of the particle. Such viscoelastic focusing of rigid particles has been demonstrated experimentally Karnis and Mason (1966); Gauthier, Goldsmith, and Mason (1971a, b); Tehrani (1996); Leshansky et al. (2007); Kim et al. (2012); Kim and Kim (2016) and using computer simulations D’Avino et al. (2012); De Santo et al. (2014); Nikoubashman et al. (2014a, b); Howard, Panagiotopoulos, and Nikoubashman (2015); Trofa et al. (2015). Interestingly, a neutral surface separating focusing points at the channel center and at the walls was discovered in simulations for certain classes of viscoelastic fluids D’Avino et al. (2012). Droplets under shear are also known to migrate in polymer solutions Hur et al. (2011); Gauthier, Goldsmith, and Mason (1971a, b); Chan and Leal (1979).

Most prior theoretical descriptions Chan and Leal (1979); Leal (1980); Leshansky et al. (2007); D’Avino, Greco, and Maffettone (2017) and simulations D’Avino et al. (2012); De Santo et al. (2014); Trofa et al. (2015) of viscoelastic focusing have adopted a continuum-level description. Such models neglect microscopic details and fluctuations of the macromolecular components of the viscoelastic medium and the particle or droplet. However, in microfluidic and nanofluidic devices, it can be necessary to consider such motion and interactions. For example, Brownian motion leads to considerable scattering in the position of a rigid sphere around the Segré–Silberberg annulus for Poiseuille flow in a pipe Prohm, Gierlak, and Stark (2012). Moreover, Brownian particles are often comparable in size to the macromolecular constituents of non-Newtonian fluids. At these length scales, Brownian spheres can exhibit anomolous motion in polymer solutions Mackay et al. (2003); Wong et al. (2004); Tuteja et al. (2007); Poling-Skutvik, Krishnamoorti, and Conrad (2015), which has been attributed to coupling between the motion of the sphere and the polymers Chen et al. (2018, 2019). It is then unclear whether well-established results for viscoelastic focusing of larger particles directly transfer to smaller particles in microchannels.

We previously demonstrated the applicability of viscoelastic focusing for Brownian rigid spheres with sizes comparable to the constituent polymer chains of a viscoelastic medium Nikoubashman et al. (2014a, b); Howard, Panagiotopoulos, and Nikoubashman (2015). However, we noted significant fluctuations of the particle around its focused position, in qualitative agreement with microfluidic experiments Kim et al. (2012). It is desirable to exploit the viscoelastic focusing mechanism to manipulate small Brownian droplets, which fluctuate in shape in addition to position, in microchannels. To our knowledge, this problem has gone relatively unexplored. We hypothesize that similar considerations may apply to the droplets as for the rigid spheres: namely, focusing onto the centerline should be improved by longer polymers and higher polymers concentrations. However, as for the rigid spheres, the distribution of the droplet position in the channel may be broad or narrow under certain flow conditions.

In this article, we test and confirm this hypothesis using particle-based computer simulations of the cross-stream migration of a Brownian fluid droplet in a polymer solution under Poiseuille flow. Although the droplet migrated outward in a Newtonian solvent (in agreement with prior simulations Marson et al. (2018)), we found that it focused onto the channel centerline in solutions of sufficiently long polymers at modest concentrations. The flow-rate dependence of the focusing was nontrivial due to a combination of effects from droplet deformation and the elastic force exerted by the polymers. We also varied the viscosity ratio between the droplet and the solvent, but did not observe any significant effect on the viscoelastic focusing in the flow regime considered. The droplet shape was controlled primarily by the local shear rate near the droplet, and so its average shape depended sensitively on the droplet distribution in the channel.

The rest of this article is organized as follows. We first describe the simulation model, including characterization of the fluid surface tension and viscosity. We then report our results, first analyzing the simulated flow fields and then systematically demonstrating the effects of polymer concentration, polymer chain length, and flow rate on the distribution of the droplet in the channel and its shape. We finally present our conclusions, suggesting avenues for future inquiry.

II Simulation model

A single fluid droplet was simulated in a Newtonian solvent and in a polymer solution using dissipative particle dynamics (DPD) simulations Hoogerbrugge and Koelman (1992); Español and Warren (1995); Groot and Warren (1997) DPD is a particle-based mesoscale simulation method that faithfully resolves hydrodynamic interactions, incorporates thermal fluctuations, and is well-suited for modeling multiphase fluids. In DPD, particles interact with each other through three pairwise forces: a conservative force $\mathbf{F}_{\rm C}$, a dissipative force $\mathbf{F}_{\rm D}$, and a random force $\mathbf{F}_{\rm R}$. As is typical, we modeled the conservative force acting on particle $i$ due to particle $j$ by a soft repulsion Groot and Warren (1997),

[TABLE]

where $a_{ij}$ sets the strength of the repulsion between particles $i$ and $j$, $r$ is the distance between the particle centers, $\mathbf{\hat{r}}$ is the unit vector to the center of particle $i$ from the center of particle $j$, and $r_{\rm c}$ is the cutoff radius for the interaction that sets the effective size of the particles.

The random and dissipative forces impart thermal fluctuations and drag while also acting as a thermostat on the DPD particles. These forces are applied in a pairwise manner that conserves momentum, with the forces on particle $i$ from particle $j$ given by

[TABLE]

where $\gamma_{ij}$ is the drag coefficient between particles $i$ and $j$, $w$ is a weight function, and $\Delta\mathbf{v}=\mathbf{v}_{i}-\mathbf{v}_{j}$ is the difference in the velocities of particles $i$ and $j$. To satisfy the fluctuation–dissipation theorem Español and Warren (1995), $\xi$ is an independent random variable for each pair of particles that has zero mean, $\langle\xi(t)\rangle=0$, and a variance $\langle\xi(t)\xi(t^{\prime})\rangle=2k_{\rm B}T\delta(t-t^{\prime})$ with $k_{\rm B}$ being Boltzmann’s constant and $T$ being the temperature. In this work, the drag coefficients were assigned per particle, $\gamma_{i}$, and the effective drag coefficient for a pair was determined by the mixing rule $\gamma_{ij}=2/(1/\gamma_{i}+1/\gamma_{j})$ Visser, Hoefsloot, and Iedema (2006).

The weight function $w$ modulates the dynamic properties, i.e., diffusivity and viscosity, of the fluid. We used the generalized weight function proposed by Fan et al. Fan et al. (2006),

[TABLE]

with $s=1/2$. This choice of $s$ increases the Schmidt number of the fluid compared to the standard DPD weight function Groot and Warren (1997) ($s=2$) to give a value closer to that of a real liquid. We also found that using $s=1/2$ gave better agreement with the no-slip boundary conditions at the microchannel walls than using $s=2$ (see below).

II.1 Fluid model

The polymer solution and droplet were modeled using three types of DPD particles: solvent (s) particles, polymer segment (p) particles, and droplet (d) particles. The model and results in this article will be reported in a fundamental system of units using $d$ as the unit of length, $m$ as the unit of mass, and $\varepsilon$ as the unit of energy, which gives $\tau=\sqrt{md^{2}/\varepsilon}$ as the unit of time. Throughout, the total density of DPD particles was $\rho=3.0/d^{3}$, all DPD particles had equal mass $1.0\,m$, the temperature was $T=1.0\,\varepsilon/k_{\rm B}$, and the cutoff radius was $r_{\rm c}=1.0\,d$. All simulations were performed using HOOMD-blue Anderson, Lorenz, and Travesset (2008); Glaser et al. (2015); Phillips, Anderson, and Glotzer (2011) (version 2.2.5) on multiple graphics processing units with a simulation time step of $0.01\,\tau$.

In order to choose the DPD repulsive parameters, we first computed the surface tension $\sigma$ between coexisting slabs of solvent and droplet particles. We fixed the repulsive parameter for particles of the same type to standard DPD values Groot and Warren (1997), $a_{\rm ss}=a_{\rm dd}=25\,\varepsilon/d$, but varied the cross-interaction strength, $a_{\rm sd}$. The drag coefficient should not affect the measured surface tension, which is a static property, and so was fixed to $\gamma_{\rm s}=\gamma_{\rm d}=1.0\,m/\tau$ to promote fast diffusion. The coexisting slabs were equilibrated by joining two cubic regions of edge length $30\,d$ to give an orthorhombic box centered around the origin with edge lengths $L_{x}=30\,d$, $L_{y}=30\,d$, and $L_{z}=60\,d$, where $x$, $y$, and $z$ denote the Cartesian coordinate axes. Particles were allowed to interdiffuse for $5\times 10^{4}\,\tau$ to equilibrate the joined slabs.

In this geometry, $\sigma$ can be computed from the pressure anisotropy Kirkwood and Buff (1949); Frenkel and Smit (2002),

[TABLE]

where $p_{\alpha\alpha}$ denotes the diagonal component of the stress tensor for index $\alpha$, and the prefactor of $1/2$ accounts for the presence of two interfaces due to the periodic boundary conditions. The cross-interaction strength was varied from $a_{\rm sd}=40\,\varepsilon/d$ to $100\,\varepsilon/d$, and the surface tension was measured using eq. (5) by sampling $p_{\alpha\alpha}$ every $0.05\,\tau$ during a $10^{5}\,\tau$ simulation. The measured surface tension (Fig. 1) is in good agreement with Groot and Warren’s empirical equation Groot and Warren (1997),

[TABLE]

with $\chi=0.286(a_{\rm sd}-a_{\rm ss})$ being their fit to the Flory–Huggins parameter for $\rho=3.0/d^{3}$.

As expected, the surface tension increased with increasing $a_{\rm sd}$ because the solvent and droplet particles became less miscible. The solvent and droplet particle density profiles near the interface (inset of Fig. 1) converged to similar values with increasing $a_{\rm sd}$. We desired a droplet that was sparingly soluble in the solvent but that could still deform under the flow rates accessible in the simulations. At $a_{\rm sd}=60\,\varepsilon/d$, the density of solvent particles dissolved in the droplet phase was already small ($4.66\times 10^{-4}/d^{3}$), and the surface tension $\sigma=2.65\,\varepsilon/d^{2}$ permitted modest deformation under viable flow rates. We accordingly selected $a_{\rm sd}=60\,\varepsilon/d$ for the cross-interaction strength.

We subsequently measured the shear viscosity of the solvent using reverse nonequilibrium simulations (RNES) Müller-Plathe (1999). Details of this method are well-described elsewhere Müller-Plathe (1999); Statt, Howard, and Panagiotopoulos (2018). We simulated a cubic box of edge length $20\,d$ containing only solvent particles with drag coefficients that varied from $\gamma_{\rm s}=1.0\,m/\tau$ to $50.0\,m/\tau$. Using RNES, we imposed a shear stress $\tau_{zx}$ on the solvent by periodically exchanging the $x$-momenta of one pair of particles from slabs of width $1.0\,d$ centered at $z=\pm 5\,d$. The swapped particles were the ones that most opposed the desired direction of flow ($x$) in each slab. We measured the velocity profile $u_{x}(z)$ between the exchange slabs ($|z|<3.5\,d$) every $10\,\tau$ over a $5\times 10^{5}\,\tau$ simulation, obtaining a Couette flow profile with a shear rate $\dot{\gamma}=\partial u_{x}/\partial z$ that decreased as the time between exchanges was increased from $0.05\,\tau$ to $0.5\,\tau$. The imposed shear stress was proportional to the measured shear rate, $\tau_{zx}=\mu_{\rm s}\dot{\gamma}$, as expected for a Newtonian fluid. The shear viscosity, $\mu_{\rm s}$, was then determined by a linear fit of $\tau_{zx}$ versus $\dot{\gamma}$.

As expected, the viscosity increased with increasing $\gamma_{\rm s}$ (Fig. 2). The simulated viscosity was generally lower than theoretically estimated Fan et al. (2006),

[TABLE]

particularly at high values of $\gamma_{\rm s}$. We selected $\gamma_{\rm s}=4.5\,m/\tau$ for the solvent particles, which gives a measured shear viscosity of $\mu_{\rm s}=1.73\varepsilon\tau/d^{3}$. In most cases, we used $\gamma_{\rm d}=\gamma_{\rm s}$, giving a droplet viscosity $\mu_{\rm d}=\mu_{\rm s}$, but we also varied $\gamma_{\rm d}$ to determine the effects of the viscosity ratio $\mu_{\rm d}/\mu_{\rm s}$ in select cases.

To model linear polymer chains of length $M$ that were fully soluble in the solvent but insoluble in the droplet, the polymer segment (p) particles were treated as if they were solvent (s) particles in the DPD interactions, i.e., $a_{\rm pp}=a_{\rm ps}=25\,\varepsilon/d$, $a_{\rm pd}=60\,\varepsilon/d$, and $\gamma_{\rm p}=4.5\,m/\tau$. Bonds within a chain were modeled by adding a harmonic spring force $\mathbf{F}_{\rm B}$ to $\mathbf{F}_{\rm C}$ for connected pairs of particles. The force on particle $i$ bonded to particle $j$ was

[TABLE]

with spring constant $\kappa=100\,\varepsilon/d^{2}$ and $b=0.7\,d$ Kranenburg, Nicolas, and Smit (2004).

II.2 Flow in microchannel

We simulated gravity-driven flow of the droplet and polymer solution in a parallel plate microchannel. The full system was initialized as follows. We first dispersed solvent particles with a total density of $\rho=3.0/d^{3}$ into a three-dimensional, periodic simulation box of dimensions $L_{x}=80\,d$, $L_{y}=40\,d$, and $L_{z}=42\,d$ and equilibrated the solvent for $1000\,\tau$. We chose $x$ as the direction of flow in the microchannel, and the parallel plates had normals along $z$. We constructed the microchannel walls by freezing solvent particles having $|z|\geq H=20\,d$, zeroing their velocities, and switching their types to be wall (w) particles Pivkin and Karniadakis (2005); Fedosov, Pivkin, and Karniadakis (2008). (The total channel width was $2H$.) The wall particles interacted with the fluid as if they were solvent particles, i.e., $a_{\rm sw}=a_{\rm pw}=25\,\varepsilon/d$, $a_{\rm dw}=60\,\varepsilon/d$, and $\gamma_{\rm w}=4.5\,m/\tau$. Mutual DPD interactions between wall particles were excluded. To help enforce no-slip and no-penetration boundary conditions at the walls, solvent, polymer, and droplet particles were additionally reflected from the planes at $z=\pm H$ using bounce-back rules Fedosov, Pivkin, and Karniadakis (2008); Revenga, Zúñiga, and Español (1999).

We selected particles near the origin of the channel to form a droplet of radius $R=4.0\,d$, giving a droplet blockage ratio of $R/H\approx 0.2$. Due to the sparing solubility of the droplet particles in the solvent, we first estimated the number of particles required to form such a droplet volume using the lever rule with the coexistence densities shown in Fig. 1. This procedure gave a droplet with a radius initially larger than $R$, but some particles later dissolved into the solvent so that the droplet reached its target radius. We then randomly created linear polymers of length $M$ from the remaining solvent particles. To build each chain, we first randomly removed $M$ solvent particles. They were reinserted as polymer segment (p) particles between the channel walls in a randomly generated chain conformation having a bond length of $0.7\,d$ between connected particles. The number of polymer chains $N_{\rm p}$ was chosen to give the desired polymer weight fraction, $\phi_{\rm p}=N_{\rm p}M/\rho V$, where $V=2L_{x}L_{y}H$ is the volume of the microchannel. In most simulations, we used $\phi_{\rm p}=5.0\%$ or $10.0\%$, but also tested $\phi_{\rm p}=0.0\%$ (no polymer), $2.5\%$ and $7.5\%$ for selected conditions. The complete configuration, including the solvent, droplet, and polymers, was equilibrated for $5000\,\tau$.

Flow was generated by applying a constant body force, $f_{\rm x}$, in the $x$-direction for all solvent, polymer, and droplet particles. For the pure solvent, applying such a force in conjunction with no-slip boundary conditions at the channel walls gives the standard parabolic (Poiseuille) velocity field,

[TABLE]

where $U=\rho f_{x}H^{2}/2\mu_{\rm s}$ is the maximum velocity at the channel centerline for this flow field. To help enforce the wall boundary conditions in the simulations, the frozen wall particles were assigned velocities $v_{x}(z)=-u_{x}(2H-|z|)$ based on their positions in the wall Fedosov, Pivkin, and Karniadakis (2008). Additionally, $u_{x}(z)$ was initially superimposed onto the polymer solution and droplet to accelerate the approach to a steady flow profile during a $5000\,\tau$ simulation. Fig. 3 shows an example configuration for the polymer solution under flow at steady state.

We repeated this procedure 5 times for each combination of chain length $M$, polymer concentration $\phi_{\rm p}$, body force $f_{\rm x}$, and droplet viscosity $\mu_{\rm d}$ studied to generate independent starting configurations. Production simulations of $10^{5}\tau$ were performed for each configuration. The droplet properties were sampled every $50\,\tau$, while the properties of the entire solution were recorded every $2500\,\tau$. The computational workflow and data were managed using the signac framework Adorf et al. (2018).

III Results and discussion

III.1 Flow field

We first measured the average flow field in the microchannel, including the solvent, polymers, and droplet. The flow was unidirectional, $u_{x}(z)$, and is shown for various polymer chain lengths at the largest polymer concentration simulated ($\phi_{\rm p}=10.0\%$) in Fig. 4a and for various concentrations of the longest polymers simulated ($M=80$) in Fig. 4b. The body force in Fig. 4 was $f_{x}=0.005\,\varepsilon/d$, which was the largest value we simulated and where any wall slip or non-Newtonian flow effects should be most pronounced. This upper bound for $f_{x}$ in our simulations was determined by trial and error so that no droplet breakup occurred.

The velocity profile in the absence of polymer was parabolic, as expected, and was also in quantitative agreement with eq. (9) using the measured $\mu_{\rm s}$ (Fig. 2). This indicates that the no-slip boundary conditions are well-enforced and also validates the RNES measurement of $\mu_{\rm s}$. The addition of polymers with $M=10$ resulted in a lower maximum velocity $U$ at the centerline, consistent with the expected higher viscosity of a polymer solution (Fig. 4a) Rubinstein and Colby (2003). Increasing the length of the polymers from $M=10$ to $M=80$ further lowered $U$. Additionally, the velocity profiles became less parabolic and developed a flattened region near $z=0$, consistent with an increasingly non-Newtonian character of the fluid. Similar trends were observed when varying the concentration of the $M=80$ polymers from $\phi_{\rm p}=0.0\%$ to $10.0\%$ (Fig. 4b), with higher polymer concentrations giving less parabolic flow profiles.

To quantify this non-Newtonian behavior, we modeled the shear stress in the polymer solutions using a power law, $\tau_{zx}=K\dot{\gamma}^{n}$, where $n$ is the flow behavior index and $K$ is a prefactor giving correct dimensions to $\tau_{zx}$. In a parallel plate channel, the flow field for a power-law fluid is

[TABLE]

A Newtonian solvent has $n=1$ and $K=\mu_{\rm s}$, and eq. (10) reduces to eq. (9), whereas shear-thinning fluids have $n<1$. We determined $K$ and $n$ by fitting the flow fields in Fig. 4 through eq. (10), recovering exponents ranging from $n=0.92$ for $M=10$ to $n=0.66$ for $M=80$ when $\phi_{\rm p}=10.0\%$ (Fig. 4a). Likewise, $n$ decreased from $0.81$ at $\phi_{\rm p}=2.5\%$ to $0.69$ at $\phi_{\rm p}=7.5\%$ for the $M=80$ polymers (Fig. 4b). There is a small but noticeable deviation of the measured velocity from the fit using eq. (10) for $M=80$ at $\phi_{\rm p}=10.0\%$, suggesting that the shear stress may have a more sophisticated functional form than the power-law model. Nonetheless, the fitted exponents give us a useful qualitative characterization of the polymer solutions.

Longer polymer chains shear thin more readily under flow than shorter chains because they have longer relaxation times, $\tau_{\rm p}$, that cause them to deform and align with the flow at smaller $\dot{\gamma}$ Doi and Edwards (1986). The dimensionless Weissenberg number, ${\rm Wi}=\dot{\gamma}\tau_{\rm p}$, characterizes this relationship. When ${\rm Wi}\ll 1$, the rate of deformation is slow compared to the polymer relaxation and primarily coil conformations are expected, whereas for ${\rm Wi}\gg 1$, the polymers are expected to be highly deformed. We approximate the shear rate by $\dot{\gamma}\approx U/H$, and estimate the polymer relaxation time from the Zimm model for a Gaussian polymer chain in dilute solution Doi and Edwards (1986), $\tau_{\rm p}\approx\mu_{\rm s}b^{3}M^{3/2}/k_{\rm B}T$. We find that ${\rm Wi}\approx 1.3$ for the $M=10$ polymers and ${\rm Wi}\approx 24$ for the $M=80$ polymers for the conditions in Fig. 4a ($f_{x}=0.005\,\varepsilon/d$). Hence, more significant shear-alignment is expected for the longer chains, which should result in more shear thinning (smaller values of $n$). This expectation is consistent with the shape of the flow fields in Fig. 4 and the fitted values for $n$.

III.2 Droplet distribution

Having characterized the flow in the microchannel, we measured the center-of-mass position of the droplet between the channel walls, $z_{\rm c}$. The droplet was identified for each configuration using a clustering procedure Pedregosa et al. (2011); Ester et al. (1999) in order to exclude droplet particles dissolved in the solvent from subsequent analysis. We analyzed the absolute value $|z_{\rm c}|$ based on the symmetry of the microchannel and to improve sampling. Previous studies Marson et al. (2018); Pan et al. (2016) have reported the average center-of-mass position, $\langle|z_{\rm c}|\rangle$, which is the first moment of the distribution of $|z_{\rm c}|$. However, a Brownian droplet can adopt a variety of distributions in the channel depending on the conditions, and we found that $\langle|z_{\rm c}|\rangle$ was not sufficiently discriminating between these. For example, a uniformly distributed droplet has $\langle|z_{\rm c}|\rangle\approx(H-R)/2$, which is indistinguishable from a droplet which is strongly focused at this position throughout the entire simulation. We accordingly computed the distribution of $|z_{\rm c}|$ using a bin size of $1.0\,d$, and will focus most of our discussion around such distributions.

We first considered the distribution of the droplet in solutions of $M=80$ polymer chains of increasing concentration $\phi_{\rm p}$ (Fig. 5) for the flow conditions shown in Fig. 4b. In the neat solvent, the droplet migrated outward from the channel centerline, showing a strongly preferred position of $|z_{\rm c}|\approx 7.5\,d$. Such outward migration is consistent with prior theoretical and simulation work for droplets Chen et al. (2014); Marson et al. (2018); Pan et al. (2016). The droplet can migrate by two mechanisms: (1) deformation due to the flow, even in the creeping flow limit, and (2) lift forces at finite inertia. We define a channel Reynolds number, ${\rm Re}=2\rho UH/\mu_{\rm s}$, and a droplet Reynolds number ${\rm Re}_{\rm d}={\rm Re}(R/H)^{2}$ Di Carlo et al. (2007); Marson et al. (2018). When the channel Reynolds number is sufficiently small, the flow is expected to be laminar. When ${\rm Re}_{\rm d}$ is small, inertial forces on the droplet are not significant and results from the Stokes flow limit are expected to apply. As an upper bound, we find ${\rm Re}\lesssim 100$ and ${\rm Re}_{\rm d}\lesssim 4$ for the investigated flow rates, consistent with the laminar flow of Fig. 4 but suggesting that inertial lift on the particle may be significant. This estimate is in accord with the observed migration of the droplet away from the centerline.

The addition of polymers to the channel at increasing polymer concentration $\phi_{\rm p}$ dramatically altered the preferred position of the droplet. The droplet distribution significantly broadened at the lowest concentration ($\phi_{\rm p}=2.5\%$). Interestingly, this included an increased probability of finding the droplet near the wall, beyond the preferred peak in the neat solvent, which we speculate may be partially due to polymer-mediated depletion interactions Asakura and Oosawa (1954). Depletion, often discussed in the context of rigid spherical colloids in solution with smaller polymer chains, induces an effective attraction between otherwise hard particles (the colloids) due to volume exclusion of a secondary species (the polymers) Russel, Saville, and Schowalter (1989). In this case, the effective attraction is between the droplet and the wall because the polymers are insoluble in the droplet and cannot penetrate the channel boundaries. Such an attraction near the wall was more pronounced for simulations at smaller $f_{x}$ (not shown here).

Continuing to add polymer increasingly focused the droplet onto the channel centerline with a narrowing distribution of $|z_{\rm c}|$. The increased polymer concentration had competing effects on the droplet migration. On the one hand, the depletion force scales with $\phi_{\rm p}$ Asakura and Oosawa (1954), which increases the outward force on the droplet near the wall. On the other hand, the increased $\phi_{\rm p}$ lowered the maximum velocity in the channel, decreasing the outward inertial lift on the droplet Di Carlo et al. (2007). Concurrently, the increased polymer concentration also increased the inward elastic force on the droplet Howard, Panagiotopoulos, and Nikoubashman (2015). The net result of these interactions is an increased inward force for larger $\phi_{\rm p}$, which improves the droplet focusing onto the centerline, consistent with with our previous work on viscoelastic focusing of rigid particle Nikoubashman et al. (2014a, b); Howard, Panagiotopoulos, and Nikoubashman (2015).

We next considered the impact of chain length $M$ on droplet focusing at two polymer concentrations, $\phi_{\rm p}=5.0\%$ and $\phi_{\rm p}=10.0\%$. Longer chains are expected to have better droplet focusing for three reasons: (1) the elastic force should scale with $M$ Howard, Panagiotopoulos, and Nikoubashman (2015), (2) longer chains deform at lower shear rates and so act more non-Newtonian, and (3) the maximum velocity was found to be lower for longer chains, reducing the outward inertial lift. The measured distributions of $|z_{\rm c}|$ (Fig. 6) are clearly consistent with this hypothesis. The addition of polymers with $M=10$ did not have a significant impact on the droplet distribution compared to the neat solvent. This may not be surprising given that ${\rm Wi}=1.3$ for the $M=10$ polymers, and the solution is nearly Newtonian. However, adding polymers of increasing length improved the focusing onto the centerline in a monotonic fashion for a given concentration.

We note, however, that there are additional concentration effects that influenced when polymers of a given size became effective focusers. This is most apparent for the $M=40$ chains. At $\phi_{\rm p}=5.0\%$, the droplet had a broad distribution of $|z_{\rm c}|$ and a most probable position of $8.5\,d$. However, at $\phi_{\rm p}=10.0\%$, the droplet was strongly focused onto the centerline. We speculate that this difference in behavior is due to an increase in elastic force with concentration, which was sufficient to overcome the inertial lift at $\phi_{\rm p}=10.0\%$ but too weak at $\phi_{\rm p}=5.0\%$.

To understand this flow rate and concentration dependence in more detail, we computed the droplet distribution for the $M=40$ chains at varying $f_{x}$ for $\phi_{\rm p}=5.0\%$ and $10.0\%$, which we compare to the distributions without any polymer (Fig. 7). Without polymer, the droplet initially migrated outward as $f_{x}$ increased, but the peak of this distribution moved inward with additional increases in $f_{x}$ (Fig. 7a), consistent with the simulations of Marson et al.Marson et al. (2018) At $\phi_{\rm p}=5.0\%$ (Fig. 7b), there was an initial trend to focus when $f_{x}\lesssim 0.003\,\varepsilon/d$. However, at larger $f_{x}$, the droplet began to migrate outward, suggesting that inertial lift dominated over the available elastic force. In contrast, the droplet distribution sharpened around the channel centerline at $\phi_{\rm p}=10.0\%$ for all $f_{x}$ considered here (Fig. 7c). It is possible that there is a sufficiently large $f_{x}$ that could exceed the inward elastic force at this concentration. However, the droplet may breakup under shear before such a force can be applied.

We finally tested the sensitivity of the viscoelastic focusing to the viscosity ratio between the droplet and the solvent. In the Stokes flow limit, Chan and Leal showed that inward or outward droplet migration can be obtained in a Newtonian solvent based on this ratio Chan and Leal (1979). However, recent simulations by Marson et al. suggest that such differences may not be as significant in the inertial regime Marson et al. (2018). Our primary interest is in how the droplet focusing may change in a non-Newtonian polymer solution. We varied the droplet viscosity ratio, $\mu_{\rm d}/\mu_{\rm s}$, from 0.54 ($\gamma_{\rm d}=1.0\,m/\tau$) to 5.3 ($\gamma_{\rm d}=40.0\,m/\tau$) for the $M=80$ polymer solutions at $\phi_{\rm p}=10.0\%$, which focused the droplet when $\mu_{\rm d}/\mu_{\rm s}=1.0$. Since the effective viscosity of the polymer solution is higher than that of the pure solvent, $\mu_{\rm d}/\mu_{\rm s}$ should be considered an upper bound on the viscosity ratio between the droplet and the polymer solution.

We found no significant differences between the droplet distributions under these conditions, and so we omit the data here for brevity. This result may not be unexpected given the qualitative picture of the focusing mechanism. The polymers are insoluble in the droplet, and so they primarily influence the fluid around it. (The flow inside the droplet is affected by $\mu_{\rm d}$, but such effects may be secondary.) Given that the viscosity ratio did not significantly alter the droplet distribution in the inertial regime of the pure solvent for Marson et al. Marson et al. (2018), it is then not surprising that the viscosity ratio also does not significantly change the droplet distribution in the polymer solution. Indeed, letting $\mu_{\rm d}/\mu_{\rm s}\to\infty$ should recover the rigid particle limit of viscoelastic focusing to which we have already drawn analogy. However, we do anticipate that the viscosity ratio may still influence the droplet distribution more significantly in other flow regimes (e.g., Stokes flow limit) that were not accessible to us in our simulations. In these cases, the migration forces controlled by the viscosity ratio, with inward or outward direction Chan and Leal (1979), would either work cooperatively or antagonistically with the elastic force of a sufficiently deformed polymer.

III.3 Droplet shape

We have concentrated our discussion thus far on how polymers influence the droplet distribution in the microchannel, but have not yet considered how the polymers influence the droplet shape and orientation in the flow. To characterize the droplet shape, we first computed its gyration tensor $\mathbf{G}$,

[TABLE]

where $\Delta\mathbf{r}_{i}$ is the vector from the droplet center of mass to particle $i$, $\alpha$ and $\beta$ are indices in the usual tensor notation, and $N_{\rm d}$ is the number of particles in the droplet. We then computed the eigenvalues $\lambda_{\alpha}^{2}$ of $\mathbf{G}$, whose corresponding eigenvectors give the principle moments of the droplet, and sorted them in descending order, $\lambda_{1}^{2}\geq\lambda_{2}^{2}\geq\lambda_{3}^{2}$. We determined the Taylor deformation parameter Taylor (1932, 1934), a dimensionless measure of the asphericity of the droplet, as

[TABLE]

For a sphere, $\mathbf{G}$ is diagonal ($\lambda_{1}=\lambda_{2}=\lambda_{3}$) and $D=0$, while a prolate spheroid has $\lambda_{1}>\lambda_{2}=\lambda_{3}$ and $D\to 1$ when the aspect ratio between the major and minor axes of the spheroid increases. Hence, larger values of $D$ correspond to droplets that have more significant deformation. We additionally determined the inclination angle of the droplet relative to the flow direction, $\theta$, using $\mathbf{G}$ Ripoll, Winkler, and Gompper (2006):

[TABLE]

$\theta\approx 0^{\circ}$ for a sphere (no preferred orientation) or for an object completely aligned with the flow, but $\theta\neq 0^{\circ}$ for particles that align with a relative tilt.

It is well-established that $D$ increases with $\dot{\gamma}$ for a droplet in an unbounded shear flow Taylor (1932, 1934); Rallison and Acrivos (1978). In Poiseuille flows, the shear rate varies across the channel, and accordingly, the droplet may experience a different deformation based on its lateral position. Fig. 8a shows the average deformation $\langle D\rangle$ versus the average center-of-mass position of the droplet $\langle|z_{\rm c}|\rangle$, clearly indicating that the droplet is (on average) more deformed when it is (on average) farther from the centerline, where the shear rate is higher. The droplet is additionally (on average) more inclined relative to the flow when it is more deformed (Fig. 8b).

It is tempting to find a parameter to collapse the data in Fig. 8 onto a single curve, e.g., using the capillary number Rallison and Acrivos (1978); Pan, Phan-Thien, and Khoo (2014); Marson et al. (2018). Unfortunately, such an analysis is again considerably complicated using only average quantities because of the droplet distribution. Indeed, the average properties computed in Fig. 8 and in prior studiesMarson et al. (2018) are intimately connected to the droplet distribution, which sets the preferred droplet location and as a consequence, the shear rates it experiences.

To deconvolve the droplet shape from the droplet distribution, we averaged $D$ and $\theta$ as functions of the shear rate $|\dot{\gamma}|$ at the droplet center-of-mass using the flow fields measured in the simulations (Fig. 4). Figure 9 shows the results of this procedure for the pure solvent and the $M=80$ polymers with $\phi_{\rm p}=10.0\%$. For a given polymer solution, the data were well-collapsed across all body forces when plotted against $|\dot{\gamma}|$. Moreover, all data nearly collapsed onto a single curve for small $\dot{\gamma}$, i.e., near the channel centerline. However, there were some noticeable differences between different polymer solutions at larger $|\dot{\gamma}|$, which corresponded to larger $f_{x}$ and positions closer to the channel walls. This discrepancy is not surprising since in that regime wall effects on the droplet may be significant.

In general, the droplet deformation $D$ increased monotonically with shear rate (Fig. 9a), but never reached zero even as $\dot{\gamma}\to 0$ due to the finite size of the droplet. In contrast, the inclination angle $\theta$ showed a maximum at intermediate shear rates (Fig. 9b). On visual inspection of the trajectories, this change in orientation with the flow appeared to be due to alignment of the droplet with the walls. Other than these boundary effects, we found that the droplet shape was not strongly influenced by the presence of the polymers in solution. Instead, the deformation and orientation correlated strongly with the shear rate due to the imposed flow. The primary roles of the polymers in setting the average deformation (Fig. 8) were then as viscosity modifiers, altering the flow field (Fig. 4) for a given body force, and as focusers for the droplet, which caused the droplet to experience a given shear rate with higher probability.

IV Conclusions

We used dissipative particle dynamics computer simulations to show that a Brownian droplet in a dilute polymer solution migrates toward the center of a parallel-plate microchannel under gravity-driven flow. The droplet had a distribution of positions in the channel that sharpened near the center for longer polymer chains at higher concentrations, but had a nontrivial dependence on the flow rate due to droplet and polymer deformation. The average droplet shape depended on the droplet distribution because its local deformation was controlled by the shear rate. Our simulations demonstrate the applicability of the viscoelastic focusing mechanism for Brownian droplets that are comparable in size to the polymer chains in the viscoelastic fluid.

Polymer-induced migration may play an important role in droplet migration and mobility in small channels flooded with complex fluids, such as those encountered in oil recovery from geological formations or in membrane filtration. In this work, we have neglected polymer solubility and adsorption with the droplet, the presence of surfactants, complex microchannel boundaries, and collective interactions between droplets at finite concentration. Such effects are surely present in many applications, and an intriguing avenue of future research is to determine how they may assist or hinder polymer-induced droplet migration in microchannels. Viscoelastic focusing will likely also influence the migration of other rigid and deformable objects in these complex fluid mixtures, including colloids, star polymers Srivastva and Nikoubashman (2018), dendrimers, cells, and micelles. Controlling the distribution of these objects in a mixture through general inertial and viscoelastic focusing mechanisms presents an opportunity to effect a separation.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements.

We happily thank Florian Müller-Plathe and Athanassios Panagiotopoulos for discussions inspiring this research. Work by M.P.H. was supported as part of the Center for Materials for Water and Energy Systems, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0019272. T.M.T. acknowledges financial support from the Welch Foundation (Grant No. F-1696), and A.N. acknowledges financial support from the German Research Foundation (DFG) under Project No. NI 1487-2/1. The simulations are part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.