- Email: [email protected]
Slip-shear and inertial migration of finite-size spheres in plane Poiseuille flow
Computational Materials Science 176 (2020) 109542
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Slip-shear and inertial migration of finite-size spheres in plane Poiseuille flow
T
Yuanding Huang, Xuezeng Zhao, Yunlu Pan
⁎
Key laboratory of Micro-Systems and Micro-Structures Manufacturing Ministry of Education and School of Mechatronics Engineering, Harbin Institute of Technology, Harbin, Heilongjiang Province 150001, PR China
ARTICLE INFO
ABSTRACT
Keywords: Inertial migration Sphere Slip-shear Rotation Plane Poiseuille flow
The effect of the slip-shear as well as the rotation on the inertial migration of finite-size spheres with blockage ratio (ratio of the sphere diameter to the channel gap) 0.2 are investigated in a plane Poiseuille flow at Re = 100–500 using dissipative particle dynamics simulations, by means of changing the driving force, the sphere density, and the fluid-sphere boundary slip condition, respectively. Neutrally buoyant spheres with equal density with the fluid and with no-slip boundary condition are observed to equilibrated at lateral positions around midway between the channel centerline and the wall, and the positions move closer to the channel centerline with increasing Re. The spheres always lag behind the surrounding fluid. A correction to the slip is proposed by considering the effect of the distance between the sphere and the wall. A forward driving force or a larger density facilitates the forward motion of the spheres and leads to their leading relative to the surrounding fluid, which drives the spheres to lateral equilibrium positions close to the wall; and vice versa. All spheres rotate with angular velocities approximate to but a bit less than the half of the local fluid shear rate at the sphere centroid, no mater the spheres are neutrally or non-neutrally buoyant, expect the spheres much close to the wall. The lift induced by rotation is found one to two orders of magnitude less than the lift induced by slip-shear due to the small angular slip velocity between the sphere and the surrounding fluid; while the slip-shear-induced lift is found critical to determine the lateral equilibrium positions of finite-size spheres cooperating with the wallinduced lift and the shear-gradient-induced lift as Re varies. An empirical theory considering the last three lifts above mentioned is built, which interprets the relation between the slip-shear and the equilibrium positions of the spheres. Above relation holds for spheres with no-slip boundary condition, whilst for neutrally buoyant spheres with full-slip boundary condition, though they always lag behind the surrounding fluid, they migrate towards the wall as Re increases, as their hydrodynamic interactions with the fluid are changed.
1. Introduction
throughput cytometry, enrichment of target cells, and detection of infected or diseased cells (e.g., mutant red blood cells, cancer cells, etc.) [1,7,8]. For inertial migration, neutrally buoyant solid spheres flowing in a straight channel or tube Poiseuille flow, is the most common and to some extent, the simplest case, as it avoids influences such as the noncentrosymmetric particle shape, deformation of deformable particles, viscoelasticity of vesicles, cells, or non-Newtonian medium fluid, and second flow induced by channel curvature or disturbance structure, etc. [5]. However, even for this simple case, complicated phenomena have been observed by substantial subsequent studies after the pioneering work of Segré and Silberberg [6]. It was found that rigid spherical particles are effectively focused when fluid inertia is in the intermediate 1) and the range (~1 < Re < ~100 ) between the Stokes regime (Re
Neutrally buoyant (external force- and torque-free) solid particles, droplets, vesicles, and cells, can migrate laterally across streamlines and even be focused at well-defined off-center positions when flowing with the medium fluid in a channel or a tube [1–5]. This counterintuitive phenomenon was first reported by Segré and Silberberg [6] in 1961, who observed that initially uniformly dispersed rigid spheres with millimeter size in a tube Poiseuille flow concentrated into an annulus at a position around 60% of the tube radius away from the tube centerline. They attributed it to the inertia of the moving fluid, thereafter this phenomenon is referred to as inertial migration or inertial focusing. Microfluidic devices utilizing inertial migration have been used to focus, separate, and sort particles, enabling applications such as high-
⁎
Corresponding author. E-mail address: [email protected] (Y. Pan).
https://doi.org/10.1016/j.commatsci.2020.109542 Received 18 November 2019; Received in revised form 13 January 2020; Accepted 15 January 2020 0927-0256/ © 2020 Elsevier B.V. All rights reserved.
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
turbulent regime (Re~2000 ) [5], where Re = f U m Lc / µ f is the channel Reynolds number, with f being the density of the fluid, Um being the maximum fluid velocity, Lc being the channel characteristic dimension, or rather, the inner diameter of a tube or the hydraulic diameter of a square or rectangular channel, µ f being the dynamic viscosity of the fluid. Rigid spheres are gathered into an annulus in a tube [9,10], while in a square channel, they tend to be focused at the four centers of the four faces [11–14], and at the four corners at higher Re [13,15,16]. The equilibrium positions in rectangular channels generally appear at the centers of the two long faces [16–22], and other positions at the centers of the short faces [18,21,22] or alongside the long faces [20] can emerge depending on Re and the aspect ratio of the channel cross section. The positions move closer to the centerline with increasing blockage ratio due to the size confinement, where is the ratio of the sphere diameter d to the tube diameter or to the smallest gap H of a 0.2 , their rectangular channel. For spheres with blockage ratio equilibrium positions are around 60% of the distance between the tube/ channel centerline and the wall [6,9–15,17,18,22–29]. When Re 100 , the positions move closer to the wall with increasing Re as the increasing flow inertia produces a stronger wall-directed shear-gradientinduced lift force [24,25]. However, when Re exceeds a critical value on the order of hundreds, with smaller values for higher blockage ratio , experimental [9,13,18] and numerical [14,22,28] studies show that the spheres begin to shift back towards the centerline, and concentrate to the so-called ’inner equilibrium position’ [9]. This reversal is surprising for it contradicts with predictions of the traditional asymptotic theories [9,24,25], who asserted that the spheres always migrate towards the wall with increasing Re for Re 3000 . Matas et al. [9] tried to explain this reversal through asymptotic theory but failed. Chun et al. [15] simply attribute it to the multibody interactions. Thus far no convincing explanations for this phenomenon have been given yet. Accordingly, we attempt to explore its underlying causes in this study. Theoretical study based on the matched asymptotic expansions by Ho and Leal [23] claimed that the inertial migration of rigid spheres in a Poiseuille flow results from the balance between a wall-directed shear-gradient-induced lift force Fs due to the parabolic velocity profile and a center-directed wall-induced lift force Fw due to the wall repulsion, and the lift was given as
FL = Fw + Fs =
2 4 f Um d 16H 2
( *2G1 +
G2),
s = p + 2 c is the sphere angular slip velocity, i.e. the difference between the sphere angular velocity p and the local fluid rotation 1 ( 2 c ), with c = 4Um z c / H being the flow shear rate at the sphere centroid. For a neutrally buoyant sphere in Poiseuille flow, it was found that the sphere always lags behind its surround fluid [23,31], and the sphere angular velocity is kind of less than the local fluid rotation [32], thus the lift induced by rotation directs towards the wall according to Eq. (2). However, no direct evidence supports this. Recently, Joseph and Ocando [32] claimed that the sphere angular slip velocity s is critical to the lift, which changes sign above and below the equilibrium position in Poiseuille flow, through their direct numerical simulations. Saffman [33] theoretically studied the effect of the slip-shear on small spheres in a linear shear flow, and an expression for the slipshear-induced lift force was given as 1
Fss = 1.615µ f Us d 2
8
fd
3
× Us,
(3)
where f = µ f / f is the kinematic viscosity of the fluid. He found that the lift induced by rotation [Eq. (2)] is smaller by one order of magnitude than that induced by slip-shear [Eq. (3)] for small Re s . Eq. (3) was then extended to the case in Poiseuille flow under conditions that Re 1, Re p = Re 2 1, and the sphere is not too close to either the tube axis or the wall. According to Eq. (3), the lift directs to the centerline when the sphere lags behind its ambient undisturbed fluid and directs to the wall when the sphere leads its ambient undisturbed fluid, which is thus known as Saffman effect. Above theory is qualitatively in agreement with experimental observations by Jeffrey and Pearson [34], who observed that, in a vertical tube, the spheres slightly denser than the fluid, i.e. with a slip relative to the surrounding fluid, drifted towards the wall for downward flow and moved towards the centerline for upward flow, though the spheres are with finite size ( = 0.05–0.1). Recenty, Asmolov et al. [35] generalized a theory to finite-size spheres at Re 20 , and found the slip-shear is crucial to the lift and determines the equilibrium positions of the spheres. Their theory is then validated 0.1 and 0.3. by lattice Boltzmann simulations for spheres with Though the effects of the slip-shear and the rotation on the inertial migration of neutrally buoyant spheres in Poiseuille flow have attracted some attentions in the past years, some issues are still mazy and poorly elucidated, especially their effects on the inward migration of the spheres with finite size at high Re. Early theoretical studies [9,23–25,29] cannot cope with the fluid disturbance triggered by the sphere with finite size. In experiments, it is difficult or even impossible to accurately measure the slip and the rotation velocity, especially at high Re. Thanks to the advance of computing capabilities, numerical simulation has become a powerful way to investigate hydrodynamic topics with well-controlled parameters. Here we employ a particlebased simulation method, dissipative particle dynamics (DPD), to investigate the inertial migration of neutrally buoyant spheres, with particular consideration of the effects of the slip-shear and the rotation on their equilibrium positions. DPD is a versatile method well applicable to the simulations of simple and complex fluids at mesoscopic scale [36–38], and it has been verified to be capable of the simulation of inertial migration of particles in Poiseuille flow up to Reynolds number of 500 or so by our previous studies [39–41]. We will introduce it as well as our model in the next section, following by our results and discussion, and ending with our concluding remarks.
(1)
= 4 are dimensionless shear rate and shear = 4z and where gradient normalized by Um and H, respectively, with z = 2z /H being the relative lateral position with respect to the half of the channel gap and z being the lateral coordinate; G1 and G2 are functions of z for wallinduced lift and shear-gradient-induced lift, respectively. Noteworthily, asymptotic theories [9,23–25] are based on the point-particle assumption, and the flows are not disturbed by the particle. It has shown that the asymptotic theories [9,23–25] cannot account for the inward migration of spheres at high Re, and other factors should have effects. Slip-shear-induced lift force and rotation-induced lift force, are two other lift forces exerted on the particles when they flow with the fluid. These two lifts should be significant at high Re [5], however, they used to be deemed as negligible relative to the shear-gradient-induced and wall-induced lifts [23], and are less considered, especially for spheres with finite size flowing in Poiseuille flows. Rubinow and Keller [30] calculated the lift induced by rotation in linear shear flow by involving the Stokes and Oseen expansions for small but finite Re s , where Re s = f |U s|d/ µ f , with Us = Up Uf being the relative velocity between the sphere and the fluid at the sphere centroid, which was given as Fr =
/ f,
2. Numerical algorithm and model 2.1. Numerical algorithm DPD was originally proposed by Hoogerbrugge and Koelman [42] in 1992. In DPD, fluid is composed of discrete beads, with each bead coarse-grained of a group of fluid molecules, which allows larger temporal and spatial scales to be accessed relative to traditional Molecular Dynamics simulations [43]. The motion of each bead is
(2)
where = p , with p being the rotational angular velocity of the = s [5], where sphere. Subsequently, it was corrected that 2
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
governed by Newton’s Second Law
d ri dv = vi , mi i = dt dt
Fij + F iext ,
(4)
where bead i located at position ri and with mass mi and velocity vi , is subjected to the pairwise additive inter-particle force Fij from bead j and external force Fiext . Fij consists of three center-to-center directed forces, namely a conservative force FCij , a dissipative force FD ij , and a random force FRij :
Fij = FCij + FijD + FRij .
(5)
The conservative force is set as a soft repulsion, which allows a large time step in DPD simulations, and is given as
FCij =
Aij (1 0,
rij ) rij, rij rc rij > rc
(6)
Aij where is the bead-bead repulsion coefficient, rij = ri rj, rij = rij/|rij |, and rc is the cutoff radius for inter-particle interactions. Aij is related to the compressibility of the DPD fluid, and in multiphase systems can be used to differentiate fluids of different kinds from each other. A dissipative force is introduced to transfer momentum between beads more efficiently than is possible from the soft potential alone, while a random force is applied to compensate for the loss of kinetic energy resulting from the dissipation. The dissipative force and the random force are given by FD ij =
D (r )( r · v ) r , ij ij ij ij
FRij =
R (r ) ij ij
t
1 2 rij,
Fig. 1. Model for a rigid sphere (in yellow) flowing in plane Poiseuille flow (in blue) with rigid walls (in gray).
are densely packed with a density w = 61.35rc 3 in HCP lattice to construct the wall, which is placed at z = 10rc and parallel to the x-y plane; see Fig. 1. The fluid beads are randomly distributed within the box with a density f = 4rc initially. The high bead density in the wall can effectively restrain the penetration of the fluid beads across the wall, although its thickness is small. In addition, the periodic boundary condition of the box enables the single wall to model a pair of walls perpendicular to the z direction, giving the height of the channel as H = Lz = 40rc . Wall beads are frozen in position during the course of the simulations. Multiple DPD beads packed in FCC lattice in spherical shape are connected by harmonic bonds to construct the rigid sphere. Considering the large size (compared to the effective size of a fluid bead) of the sphere in three dimensions, a lower number density n p = 32rc 3 is chosen for the sphere than that of the wall, although this density is still high enough to effectively prohibit the penetration of fluid beads into the sphere, but not so high as to make computer memory cost unmanageable. Each bead in the sphere is connected with its nearest, next-nearest, and next-next-nearest beads through a bond with a potential energy
(7)
where is the friction coefficient, is the amplitude of noise, t is the integral time step, D (rij ) and R (rij ) are rij -dependent dimensionless weighting functions, and ij is the Gaussian white noise following Gaussian statistics: and ij (t ) = 0 ij (t ) kl (t ) = ( ik jl + il jk ) (t t ) . The fluctuation–dissipation theorem linking the dissipative and random forces gives D (r )
=[
D (r )]2
and
2
= 2 kB T
(8)
as explained by Español and Warren [44]. Thus, the momentum and the temperature are conserved. D (r ) and R (r ) conventionally take the form D (r )
=[
D (r )]2
(1
=
r ) 2 , r rc 0, r > rc
(9)
An unmodified version of velocity-Verlet algorithm is used to advance the set of positions and velocities in the DPD we employed, which is given as
ri (t + t ) = ri (t ) + t vi (t ) + vi (t + t ) = vi +
1 2
1 ( 2
V (r ) =
t Fi , 1 2
t (Fi (t ) + Fi (t + t ).
l 0)2
(11)
where kb is the force constant and is set with a high value kb = 1000m p to restrain the flexibility of the bond, with m p being the mass of the sphere bead, l is the length of the bond, and l 0 is the bond rest length, with l 0 = 2 b/2, b, 3 b for bonds connecting the nearest, next-nearest, and next-next-nearest beads, respectively, with b = 0.5rc . The sphere is initially placed at the position (0, 0, 10rc) , which is at the center of the channel. Following the simulation work for Poiseuille flow in a planar channel using DPD by Millan et al. [46], we adopt Aij values that they took for our purposes, which are
t )2Fi (t ),
Fi (t + t ) = Fi (r (t + t ), v (t + t )), vi (t + t ) = vi (t ) +
1 k b (l 2
(10)
In a ‘standard’ DPD simulation [45], m , rc and kB T are normalized to unity, i.e., m = rc = kB T = 1; the unit of time is defined as = rc m / kB T = 1; the density of the fluid is taken as f = 4rc 3; and the friction coefficient between the fluid beads is set as ff = 4.5 in DPD units. 2.2. Model
w f p 3 10 kB T w Aij = f 3 25 3 , rc p 10 3
A box, with size L x × L y × Lz = 60rc × 40rc × 40rc and periodic boundary condition in all directions, is employed. The length (L x ) has been verified to be long enough for the fully-developed plane Poiseuille flow, and the width (L y ) is wide enough to minimize interactions between periodic images of the spheres [39]. Three layers of DPD beads
where the subscripts “w”, “f”, and “p” represent the DPD beads forming the wall, the fluid, and the spherical particle, respectively. The fluid–fluid repulsion parameter Aff is not far from the standard DPD value of 18.75 [45], while the other parameters (along with the choice of densities discussed above) are chosen to minimize fluid penetration and 3
(12)
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
N gives the standard error used to obtain the error bars reported here in our figures. Simulations are performed with DPD in the Highly Optimized Object-oriented Molecular Dynamics (HOOMD)-Blue software package [49,50], a GPU-based simulation toolkit. 3. Results In our simulations, the spheres are chosen to have a diameter d = 8rc , and the blockage ratio is therefore = d/ H = 0.2 . A smaller sphere is found to suffer significant Brownian motion and is hard to be focused. A larger sphere is beyond the computational ability when analyzing the data due to the overmuch beads within it. A “standard” case for the inertial migration of neutrally buoyant sphere in plane Poiseuille flow has been simulated in our previous work [39] as an effort to assess the DPD fluid. In that case, the bead mass of the sphere is one-eighth of that of the fluid (i.e., m p = m f /8 = 1/8), as is the driving force applied on per sphere bead compared to that of the fluid (i.e., fx p = fx f /8), so that equal density (note the bead number density is 8fold higher than that of fluid, i.e., n p = 8n f ) and acceleration (induced by the driving force) are achieved between the fluid and the sphere. The friction coefficient between the sphere beads and the fluid beads is set as fp = 4.5, thus there is no-slip at the fluid-sphere interface. This corresponds to the common case of the motion of a neutrally buoyant sphere free drifting in Poiseuille flow without external force or torque. The equilibrium positions of above sphere are around at middle between the channel centerline and the wall, and they move closer to the centerline as Re increases, which are shown red symbols in Fig. 3. Our results have been compared with other experimental and numerical studies for spheres with the same (=0.2) in square [22,14], rectangular [17,18,22], and plane [22] channels, respectively. Our results are in agreement with other studies for spheres in rectangular channels [17,22] and a plane channel [22] at Re = 100–200 . Based on the results
Fig. 2. Example velocity profiles of pure DPD fluid with fw = 1.0 (blue) and 4.5 (red), respectively, and are compared with the Hagen-Poiseuille equation (the black solid curve).
slip against the spherical particle and the wall. The friction coefficient between the fluid beads and the wall beads is set as fw = 4.5, with which the no-slip boundary condition is achieved at the fluid-wall interface. One example of the velocity profile from our simulation of the Poiseuille flow of pure fluid in the plane channel is shown with red symbols in Fig. 2, which is compared to the Hagen-Poiseuille equation (the black curve in Fig. 2)
u x (z ) =
1 dp 2 z 2µ f dx
H 2
2
=
f fx f
2µ f
z2
H 2
2
,
(13)
where u x (z ) is the flow velocity in the x direction (flow direction) as a function of gap coordinate z , dp /dx = f fx f is the pressure gradient in the x direction, and fxf is the external driving force on per fluid bead along the x direction. The external driving forces in y and z directions are set to zero, i.e., F ext f = [fx f , f y f , fz f ] = [fx f , 0, 0]. The good match of our simulation results with the Hagen-Poiseuille equation demonstrates that our parameter choice yields the Hagen-Poiseuille flow with a noslip boundary condition, and thus, the viscosity of the fluid can be estimated by
µf =
f fx f H
8Um
2
.
|z*eq|
0.6
0.4 Our work (nmrcl, plane channel) Lashgari et al. (nmrcl, square channel) Yuan et al. (nmrcl, square channel) Masaeli et al. (expt, rectangular channel, AR ≈ 1.6) Ciftlik et al. (expt, rectangular channel, AR = 1.6) Lashgari et al. (nmrcl, rectangular channel, AR = 2) Lashgari et al. (nmrcl, plane channel)
(14)
Decreasing the fluid-wall friction coefficient fw to a smaller value will lead to a partial-slip at the fluid-wall interface (one example of velocity profile of flow with fw = 1.0 is shown with blue symbols in Fig. 1) and to zero will result in a full-slip (no example of velocity profile can be shown as the fluid will be accelerated all the time, but the profile keeps flat). Similarly, the fluid-sphere boundary condition can be adjusted through the friction coefficient fp . 0.01, and it satisfies t = 0.1 / to keep the The time step is t/ fluctuation of the system temperature within 2% [54,48]. The system is first equilibrated by running for one million time steps to establish an accurate model. After applying driving forces, all simulations are run for one million steps to achieve a nearly steady-state position of the particle, and then for three million additional steps to collect fully developed data. After discarding the initial part of the run, the fluctuations around the equilibrium position are autocorrelated to obtain the characteristic autocorrelation time. The length of the run divided by this autocorrelation time gives the effective number of statistically independent positions N sampled by the sphere during the course of the simulation. Dividing the standard deviation of the average position by
0.2
0.0
0
100
200
300 Re
400
500
600
Fig. 3. Relative equilibrium positions z eq of neutrally buoyant spheres with p = f and no-slip boundary condition (the standard case) in plane channel
through our DPD simulations, comparing results of other studies for neutrally buoyant spheres with the same blockage ratio ( = 0.2 ) in square (Lashgari et al. [22], Yuan et al. [14]), rectangular (Masaeli et al. [17], Ciftlik et al. [18], Lashgari et al. [22]), and plane (Lashgari et al. [22]) channels, respectively. Notably, Re in some publications have different definition, which we have transformed to our definition which is based on the maximum velocity Um and the smallest gap of rectangular channel H.
4
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
0.8
0.8
(a)
axp/axf = 8 axp/axf = 4 axp/axf = 1
200
0.10
300 Re
0.4
ρp/ρf = 4
400
0.0
500
100
200
300 Re
Us [rc/τ]
Faxen correction
0.05
0.00
0.00
-0.10
-0.15
-0.15
100
200
300 Re
0.10
-0.15
0.04
0.02
0.02 300 Re
400
100
200
300 Re
400
500
0.12 (f)
0.10
500
0.00
(i)
0.08
0.06
0.04
200
500
0.08
αc/2
100
400
0.12
(c)
0.06
0.00
500
500
ωp [rad/τ]
ωp [rad/τ]
400
400
-0.10
ωp [rad/τ]
0.12
300 Re
300 Re
-0.05
-0.20
200
200
(h)
0.05
-0.10
100
0.10
-0.05
-0.05
0.08
0.0
500
(e)
0.00
0.10
400
0.10
0.05
100
Full slip Partical slip No-slip
0.2
ρp/ρf = 1
(b)
-0.25
0.4
ρp/ρf = 8 0.2
Us [rc/τ]
100
axp = 0 axp/axf = -4
0.6
|z*eq|
|z*eq|
|z*eq|
0.4
0.2
Us [rc/τ]
(g)
0.6
0.6
0.0
0.8 (d)
0.06 0.04
αc/2 100
200
300 Re
400
αc/2
0.02 500
0.00
100
200
300 Re
400
500
Fig. 4. Equilibrium positions [(a), (d), (g)], slips [(b), (e), (h)], and rotational angular velocities [(c), (f), (i)] of spheres with = 0.2 in plane Poiseuille flow. (a)–(c) are for spheres with different driving forces, (d)–(f) are for spheres with different densities, and (g)–(i) are for spheres with different fluid-sphere boundary slip conditions. The dashed lines in (c), (f), and (i) are Tollert’s theory [51] that p = c /2 .
from other studies shown in Fig. 3, the equilibrium positions move closer to the centerline with increasing AR, where AR is the aspect ratio of the cross section of the rectangular channel, which for square channels AR = 1 and for plane channels AR = . It is therefore reasonable that the positions for our spheres in plane channel are a bit closer to the centerline. Also, the trends of the inward migration are similar to that in square [14] and rectangular [18] channels. The validity of our model is thus verified. Here, the slip and rotational angular velocity of above spheres are measured and shown with red symbols in Fig. 4 (b) and (c), respectively. The slip is an average of 100 data samples when the sphere has been focused around its equilibrium position; see Fig. A.1 in Appendix
A. It is found that the spheres always lag behind their surrounding fluid, and the lag increases with increasing Re. The magnitudes of the slips are within 4% of the velocity of the ambient fluid, and 1.2–1.8 times of 2U /3 [35], with smaller times for sphere Faxen correction Us = m farther away from the wall. Asmolov et al. [35] also observed slips about twice larger than the Faxen correction for finite-size spheres 0.1) which are not too close to the wall. The deviations from the ( Faxen correction attributes to the resistance effect from the wall, which diminishes as the sphere moves away from the wall. As the slip is related to the distance between the sphere and the wall, the slip for finitesize spheres which are not too close to the wall in plane Poiseuille flow 0.27 for the spheres can be fitted as Us = k s 2Um/(1 z ) , where k s 5
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
of the standard case based on our simulation data. The rotational angular velocity is measured by tracing two beads at the two poles of one diameter of the sphere along 2 × 103 after the sphere being well focused, and one example of rotating orbit in the shear-gradient plane (x-z plane) is shown in Fig. A.2 in Appendix A. The angular velocity of the sphere p [the red curve in Fig. 4 (c)] approximately increases linearly with increasing Re and p c/2 , where c is the local flow shear rate at the sphere centroid, in accord with Tollert’s theory [51] for a free moving sphere in linear shear flow. This should attribute to the intermediate equilibrium position of the sphere between the centerline and the wall, at which the Poiseuille flow velocity profile can be approximately regarded as linear, and the wall resistance effect is slight. Also, we note that the local fluid rotation ( c /2 ) is a bit greater than the sphere angular velocity p , in agreement with the finding by Joseph and Ocando [32]. To further investigate the effects of the slip-shear, three different means, i.e., varying the driving force applied on the spheres, varying the sphere density, and changing the fluid-sphere boundary slip condition, are employed, respectively; while other parameters are kept the same as the standard case. Four different driving forces, which are fx p = fx f , fx p = fx f /2, fx p = 0 , and fx p = fx f /2, are applied on per bead of the sphere, respectively, while keeping p = f and fp = 4.5. The acceleration induced by the driving force on the sphere axp are therefore 8, 4, 0, and 4 times that on the fluid ax f . It can be seen from Fig. 4 (b) that, for ax p /a x f = 8 and 4, the sphere leads its surrounding fluid, and the leading increases as Re and axp increase. For ax p /a x f = 1, 0, and 4 , the sphere lags behind its surrounding fluid, and the lag increases with increasing Re but decreasing axp . The equilibrium positions of the spheres are close to the wall when the spheres lead the surrounding fluid, and close to the centerline when the spheres lag behind the surrounding fluid, in accordance with Saffman’s prediction [33]. The trend is more pronounced as the absolute value of the slip |Us | becomes larger. For ax p /a x f = 4 , the spheres are focused to positions around the channel centerline. This is in agreement with recent studies which, by applying an external axial force on the rigid spheres in the opposite direction to induce a large lag relative to the flow, achieved the focusing of the spheres with finite size at the channel center in Stokes regime (Re 0.05) experimentally [52] and at moderate Re (10 80 ) numerically [53], respectively. The rotational angular velocity still linearly increases with increasing Re, and approximately follows p = c/2 , especially for spheres that are away from both the wall and the channel centerline. For spheres close to the wall, their rotation rates deviate from c /2 due to the growing wall resistance; and for spheres at the channel centerline, their rotations are restrained and cannot be measured accurately as the local shear rate approaches 0. The variation of the particle density is realized by changing the mass of the sphere bead simultaneously with the driving force, in which way keeps the acceleration of the sphere induced by the driving force equal to that of the fluid. Notably, change of the mass of the sphere bead leads to unequal density between the sphere and the fluid, i.e., non-neutrally buoyant for spheres in real world. However, no gravity or buoyancy is applied in our simulation system, and therefore the spheres are still neutrally buoyant. Equilibrium positions, slips, and rotational angular velocities of spheres with densities 4 and 8 times that of fluid are shown in Fig. 4 (d)–(f) and compared with the standard case, respectively. It is shown that the increase of particle density (or inertia) accelerates the sphere, as for the spheres with p / f = 4 , though they still slightly lag behind its surrounding fluid, their lags become smaller relative to the standard case ( p / f = 1), and the spheres with p / f = 8 become to lead the surrounding fluid. The relation between the equilibrium positions and the slips still follows Saffman’s prediction [33]. It should be noted that the inward migration of the sphere with p / f = 8 at high Re (400–500) is due to the increasing fluctuation of the sphere in lateral direction, which leads to a small average of the position; see Fig. A.3 in Appendix A. The rotational angular velocity follows p = c/2 with
exceptions when the spheres are close to the wall. The slip between the sphere and its surrounding fluid is also altered by changing the fluid-sphere boundary slip condition, specifically, changing the fluid-sphere friction coefficient fp . We have shown that decreasing fw to 1 brings off a partial-slip at the fluid-wall interface, and to 0 leads to a full-slip. Similarly, we set fp = 1 and 0 to implement the partial-slip and the full-slip at the fluid-sphere interface, respectively; and then compare the results with the standard case which with no-slip fluid-sphere boundary condition ( fp = 4.5); see Fig. 3. (g) - (i). It shows that all the spheres lag behind their surrounding fluid, and the lags of the spheres with partial- and full-slip are almost the same, and are kind of smaller than that with no-slip at the same Re. It seems the boundary slip condition has slight effect on the slip of the sphere while the particle inertia play a decisive role. However, the equilibrium positions of the spheres with full-slip move close to the wall with increasing Re, which differ significantly from the trend for spheres with no-slip, while the equilibrium positions of spheres with partial-slip are kind of closer to wall than that of no-slip at the same Re. The rotational angular velocities of the spheres with full-slip deviate from p = c/2 , though the deviations are small. At low Re, where the sphere is relatively away from the wall, the full-slip boundary condition reduces the friction between the sphere and the fluid, the rotation rate is kind of larger than c /2 . At high Re, where the sphere is close to the wall, the wall resistance restrains the rotation and the angular velocity is kind of smaller than c /2 . 4. Discussion Above results show that, for spheres with no-slip boundary condition, a larger forward driving force or a larger density (or inertia) on/of the spheres than that on/of the fluid is conductive to their forward motion, and leads to the leading relative to their surrounding fluid; while a backward driving force make the spheres lag behind the surrounding fluid. As Saffman’s theory [33] predicted, the lift induced by the slip-shear directs to the wall when the sphere leads its surrounding fluid and directs to the centerline when it lags behind its surrounding fluid. Our simulation results quantitatively agree with this. The sphere angular velocity p c/2 , and is kind of less than the local fluid rotation ( c /2 ), except when the sphere is close to the wall, in accordance with the observations by Joseph and Ocando [32]. Due to the sphere angular slip velocity s is small, the rotation-induced lift [according to Eq. (2)] is found one to two orders of magnitude less than the lift induced by slip-shear [according to Eq. (3)] based on our simulation data, as shown in Fig. 5. This also holds for the spheres much close to the wall. Thus the rotation-induced lift can be neglected. Meanwhile, the slip-shear-induced lift is found comparable to the wall-induced lift and the shear-gradient-induced lift, which should play a role in the inertial migration of the sphere. The non-dimensionalized lift arising from the competition between the wall-induced lift and shear-gradient-induced lift can be written as
Fw + Fs =
Fw + Fs = G1 z *2 + G2 z 2 4 2 f Um d /H
(15)
according to Eq. (1), which is shown with black curve in Fig. 6, re0.6 vealing that the equilibrium position of the sphere is fixed at z and has no relation with the Um , i.e., Re. However, our study and most other studies have shown that the position varies with Re. Whilst the slip-shear-induced lift proposed by Saffman [33] can be non-dimensionalized as
Fss =
Fss ks = 3.23 2 4 2 1 z U d / H m f
z Re
(16)
by substituting Us = k s 2Um/(1 z ) and = 4Um z / H into Eq. (3), which indicates Fss is relative to Re. Among these three dominant lifts on spheres, the slip-shear-induced lift is therefore critical to adjust the 6
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
spheres in plane Poiseuille flow at Re over the range of 100
30
Fw Fs Fss Fr
FL
20
FL = FL/(
0 -10 200
300 Re
400
500
Fig. 5. The lifts induced by wall, shear-gradient, slip-shear, and rotation, respectively, on a neutrally buoyant sphere with p = f and no-slip boundary condition (i.e., the standard case) in plane Poiseuille flow.
0.2
FL/( fU2md4/H2)
0.0 -0.2 -0.4
Ho and Leal Re = 100 Re = 200 Re = 300 Re = 400 Re = 500
-0.6 -0.8 -1.0 0.0
0.2
0.4
z*
= Fs + Fw + css (Re) Fss.
500 as (17)
css is chosen in the form of css = g Re + h for simplicity, where g and h are a coefficient and a constant, respectively, and are set by fitting our simulation results into Eq. (17) through the least square method. css (Re) 0.0073Re 0.255 for the neutrally buoyant sphere with p = f and no-slip boundary condition (i.e., the standard case), and the non-dimensionalized lifts FL are shown with colored curves in Fig. 6 for different Re. The inward migration of neutrally buoyant spheres with p = f and no-slip boundary condition in Poiseuille flow at Re = 100 500 is thus interpreted. css as well as ks for other cases are shown in Table B.1 in Appendix B, which account for the Re-dependent inward migration for spheres lagging behind the surrounding fluid and the Re-dependent outward migration for spheres leading the surrounding fluid through Eq. (17). For non-neutrally buoyant spheres, by applying an forward or backward axial driving force, the magnitude and the direction of Us changes, resulting in a redistribution of the equilibrium positions from the channel centerline (z eq = 0 ) to the wall (z eq = 0.8); see Fig. 4 (a) and (b), respectively. This suggests a promising application, which utilizes an external driving force, such as electric, magnetic, or gravitational force, to adjust the sphere equilibrium positions within the channel and then for particle separation. Actually, Kim and Yoo [52] have achieved the focusing of rigid spheres at the channel centerline by applying a backward axial electric force on the spheres, however, at Stokes regime. One intractable limitation for this application is the strong hydrodynamic force on the sphere at high Re, which requires a strong external force to compare or overwhelm it. Changing the density of the spheres changes the particle inertia, which in turn changes the slips and the equilibrium positions of the spheres. However, this needs the change of the sort of the particles in reality and will bring about sedimentation or floatation in horizontal channels. By contrast, altering the fluid-sphere boundary slip condition brings more inspiration to the potential application for inertial migration/focusing at high Re. According to our simulation results, for a sphere with full-slip boundary condition, though it lags behind its surrounding fluid, its equilibrium position moves closer to the wall as Re increases. The change of boundary slip condition changes the hydrodynamic interaction between the sphere and its surrounding fluid, in consequence, above theories built for spheres with no-slip boundary condition do not hold here anymore.
10
100
2 4 2 f Um d / H )
5. Conclusion
0.6
0.8
We have performed dissipative particle dynamics simulations to study the effect of the slip-shear as well as the rotation on the inertial migration of finite-size spheres with blockage ratio = 0.2 flowing in a plane Poiseuille flows at Re = 100 500 , by changing the driving force, the sphere density, and the boundary slip condition between the sphere and the fluid, respectively. For neutrally buoyant spheres with p = f and no-slip boundary condition, they are equilibrated at lateral positions around midway between the channel centerline and the wall, with the positions moving closer to the centerline with increasing Re. The spheres are found always lagging behind the surrounding fluid with a 0.27 , which is improved upon the slip Us = k s 2Um/(1 z ) with k s 2U /3 [35] by considering the classical Faxen correction that Us = m distance between the sphere and the wall. A forward driving force or a large density relative to the fluid facilitates the forward motion of the sphere and leads to its leading with respect to the surrounding fluid, resulting in an outward migration of the sphere towards the wall; while a backward driving force makes the sphere lag behind the surrounding fluid, giving rise to an inward migration of the sphere towards the channel centerline. This quantitatively agrees with Saffman’s prediction [33]. Meanwhile, the rotational angular velocity p c/2 , and is
1.0
Fig. 6. Non-dimensionalized lift for theory proposed by Ho and Leal [23] FL = Fs + Fw (the black curve), and for present theory FL = Fs + Fw + css (Re)F ss (the colored curves), where css (Re) 0.0073Re 0.255 is for neutrally buoyant spheres with p = f and noslip boundary condition.
equilibrium positions of the spheres with variation of Re. However, Eq. (16) shows that the magnitude of Fss decreases with increasing Re. According to our simulations, the impact of slip-shear-induced lift increases with increasing Re over the range of 100 500 , which leads to the progressive inward migration for spheres lagging behind the ambient fluid and the progressive outward migration for spheres leading the ambient fluid. Based on above analysis, we multiply Fss by a coefficient css , which is a function of Re, and propose a lift for finite-size 7
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
kind of less than the local fluid rotation ( c /2 ), except when the sphere is close to the wall, where the wall resistance is significant. Based on previous theories [23,30,33], we found the rotation-induced lift is less by one to two orders of magnitude than the slip-shear-induced lift, while the latter is comparable to the dominant wall-induced lift and shear-gradient-induced lift, and is critical to adjust the equilibrium positions within the channel when Re varies. A theory including the lifts induced by wall, shear-gradient, and slip-shear is built, which well explains the relation between the slip-shear and the inertial migration of spheres. Above theory holds for spheres with no-slip boundary condition. Whilst for spheres with full-slip boundary condition, as their hydrodynamic interactions with the fluid are changed, they migrate towards the wall as Re increases, though they always lag behind the surrounding fluid.
CRediT authorship contribution statement Yuanding Huang: Conceptualization, Investigation, Methodology, Software, Data curation, Formal analysis, Visualization, Writing - original draft. Xuezeng Zhao: Funding acquisition, Supervision, Writing review & editing. Yunlu Pan: Funding acquisition, Project administration, Resources, Supervision, Formal analysis, Validation, Writing review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment
Data availability
This research was supported by the National Natural Science Foundation of China (Grant Nos. 51505108 and 51475118) and Key Laboratory of Micro-systems and Micro-structures Manufacturing of Ministry of Education, Harbin Institute of Technology (No. 2017KM009).
The raw data required to reproduce these findings are available to down-load from [ https://data.mendeley.com/datasets/jcjhhps7tr/ draft?a=b00e7da0-6e90-4bfc-ae36-aeb1c6c43aa1 ].
Appendix A. Axial velocity, rotating orbit, and motion trajectories of spheres
Fig. A.1. Comparison of axial velocities of neutrally buoyant spheres which with positions with the fluid velocities at different Re.
p
=
p
and no-slip boundary condition (i.e., the standard case) at their equilibrium
Fig. A.2. Example of rotating orbit of a neutrally buoyant sphere of the standard case in the shear-gradient plane (x-z plane), where mapping of a chosen diameter in x-z plane with the x direction. 8
z
is the angle between the
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
Fig. A.3. Comparison of motion trajectories of spheres with different densities at Re
500 .
Appendix B. css (Re) and k s
Table B.1 css (Re) and k s for spheres with different driving forces and different densities over the Re range 100 Cases p
=
f
ax p = ax f
ks
css (Re) ax p/ax f =
0.0006Re
4
500 .
0.21
-1.67 -0.63
ax p/ax f = 0
0.0034Re + 0.20
ax p/ax f = 1
0.0073Re
0.26
ax p/ax f = 4
0.0027Re
0.15
ax p/ax f = 8
0.0030Re + 0.17
0.23
-0.27 0.09
p/ f
=4
0.0114Re
1.60
-0.05
p/ f
=8
0.0022Re + 1.79
0.16
[15] B. Chun, A. Ladd, Inertial migration of neutrally buoyant particles in a square duct: An investigation of multiple equilibrium positions, Phys. Fluids 18 (3) (2006) 031704, , https://doi.org/10.1063/1.2176587. [16] A.A.S. Bhagat, S.S. Kuntaegowdanahalli, I. Papautsky, Enhanced particle filtration in straight microchannels using shear-modulated inertial migration, Phys. Fluids 20 (10) (2008) 101702, , https://doi.org/10.1063/1.2998844. [17] M. Masaeli, E. Sollier, H. Amini, W. Mao, K. Camacho, N. Doshi, S. Mitragotri, A. Alexeev, D. Di Carlo, Continuous inertial focusing and separation of particles by shape, Phys. Rev. X 2 (2012) 031017, , https://doi.org/10.1103/PhysRevX.2. 031017. [18] A.T. Ciftlik, M. Ettori, M. Gijs, High throughput-per-footprint inertial focusing, Small 9 (16) (2013) 2764–2773, https://doi.org/10.1002/smll.201201770. [19] J. Zhou, I. Papautsky, Fundamentals of inertial focusing in microchannels, Lab Chip 13 (6) (2013) 1121–1132, https://doi.org/10.1039/C2LC41248A. [20] C. Liu, G. Hu, X. Jiang, J. Sun, Inertial focusing of spherical particles in rectangular microchannels over a wide range of reynolds numbers, Lab Chip 15 (4) (2015) 1168–1177, https://doi.org/10.1039/C4LC01216J. [21] C. Liu, C. Xue, J. Sun, G. Hu, A generalized formula for inertial lift on a sphere in microchannels, Lab Chip 16 (5) (2016) 884–892, https://doi.org/10.1039/ C5LC01522G. [22] I. Lashgari, M.N. Ardekani, I. Banerjee, Inertial migration of spherical and oblate particles in straight ducts, J. Fluid Mech. 819 (2017) 540–561, https://doi.org/10. 1017/jfm.2017.189. [23] B.P. Ho, L.G. Leal, Inertial migration of rigid spheres in two-dimensional unidirectional flows, J. Fluid Mech. 65 (02) (1974) 365–400, https://doi.org/10.1017/ s0022112074001431. [24] J.A. Schonberg, E.J. Hinch, Inertial migration of a sphere in Poiseuille flow, J. Fluid Mech. 203 (1989) 517–524, https://doi.org/10.1017/s0022112089001564. [25] E.S. Asmolov, The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number, J. Fluid Mech. 381 (1999) 63–87, https://doi.org/ 10.1017/s0022112098003474. [26] J. Feng, H.H. Hu, D.D. Joseph, Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 2. Couette and Poiseuille flows, J. Fluid Mech. 277 (1994) 271–301, https://doi.org/10.1017/s0022112094002764. [27] D. Di Carlo, D. Irimia, R.G. Tompkins, M. Toner, Continuous inertial focusing, ordering, and separation of particles in microchannels, Proc. Nat. Acad. Sci. 104 (48) (2007) 18892–18897, https://doi.org/10.1073/pnas.0704958104. [28] X. Shao, Z. Yu, B. Sun, Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers, Phys. Fluids 20 (10) (2008) 103307, , https://doi.org/10.1063/1.3005427. [29] J.-P. Matas, J.F. Morris, E. Guazzelli, Lateral force on a rigid sphere in large-inertia
References [1] D. Di Carlo, Inertial microfluidics, Lab Chip 9 (21) (2009) 3038–3046, https://doi. org/10.1039/b912547g. [2] P. Sajeesh, A.K. Sen, Particle separation and sorting in microfluidic devices: a review, Microfluid. Nanofluid. 17 (1) (2014) 1–52, https://doi.org/10.1007/s10404013-1291-9. [3] H. Amini, W. Lee, D. Di Carlo, Inertial microfluidic physics, Lab Chip 14 (15) (2014) 2739–2761, https://doi.org/10.1039/c4lc00128a. [4] J.M. Martel, M. Toner, Inertial focusing in microfluidics, Annu. Rev. Biomed. Eng. 16 (1) (2014) 371–396, https://doi.org/10.1146/annurev-bioeng-121813-120704. [5] J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M.E. Warkiani, W. Li, Fundamentals and applications of inertial microfluidics: a review, Lab Chip 16 (1) (2016) 10–34, https://doi.org/10.1039/c5lc01159k. [6] G. Segré, A. Silberberg, Radial particle displacements in poiseuille flow of suspensions, Nature 189 (4760) (1961) 209–210, https://doi.org/10.1038/189209a0. [7] G.Y. Lee, C.T. Lim, Biomechanics approaches to studying human diseases, Trends Biotechnol. 25 (3) (2007) 111–118, https://doi.org/10.1016/j.tibtech.2007.01. 005. [8] A.J. Mach, O.B. Adeyiga, D. Di Carlo, Microfluidic sample preparation for diagnostic cytopathology, Lab Chip 13 (6) (2013) 1011–1026, https://doi.org/10.1039/ c2lc41104k. [9] J.P. Matas, J.F. Morris, É. Guazzelli, Inertial migration of rigid spherical particles in Poiseuille flow, J. Fluid Mech. 515 (2004) 171–195, https://doi.org/10.1017/ s0022112004000254. [10] Y.-S. Choi, S.-J. Lee, Holographic analysis of three-dimensional inertial migration of spherical particles in micro-scale pipe flow, Microfluid. Nanofluid. 9 (4) (2010) 819–829, https://doi.org/10.1007/s10404-010-0601-8. [11] D. Di Carlo, J.F. Edd, K.J. Humphry, H.A. Stone, M. Toner, Particle segregation and dynamics in confined flows, Phys. Rev. Lett. 102 (2009) 094503, , https://doi.org/ 10.1103/PhysRevLett. 102.094503. [12] Y.-S. Choi, K.-W. Seo, S.-J. Lee, Lateral and cross-lateral focusing of spherical particles in a square microchannel, Lab Chip 11 (3) (2011) 460–465, https://doi.org/ 10.1039/C0LC00212G. [13] K. Miura, T. Itano, M. Sugihara, Inertial migration of neutrally buoyant spheres in a pressure-driven flow through square channels, J. Fluid Mech. 749 (2014) 320–330, https://doi.org/10.1017/jfm.2014.232. [14] C. Yuan, Z. Pan, H. Wu, Inertial migration of single particle in a square microchannel over wide ranges of re and particle sizes, Microfluid. Nanofluid. 22 (9) (2018) 102, https://doi.org/10.1007/s10404-018-2120-y.
9
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
[30] [31] [32] [33] [34] [35] [36] [37]
[38] [39] [40]
[41]
laminar pipe flow, J. Fluid Mech. 621 (2009) 59–67, https://doi.org/10.1017/ S0022112008004977. S. Rubinow, J.B. Keller, The transverse force on a spinning sphere moving in a viscous fluid, J. Fluid Mech. 11 (3) (1961) 447–459, https://doi.org/10.1017/ S0022112061000640. R. Simha, Untersuchungen über die viskosität von suspensionen und lösungen, Kolloid-Zeitschrift 76 (1) (1936) 16–19, https://doi.org/10.1007/BF01432457. D. Joseph, D. Ocando, Slip velocity and lift, J. Fluid Mech. 454 (2002) 263–286, https://doi.org/10.1017/S0022112001007145. P. Saffman, The lift on a small sphere in a slow shear flow, J. Fluid Mech. 22 (2) (1965) 385–400, https://doi.org/10.1017/S0022112065000824. R.C. Jeffrey, J. Pearson, Particle motion in laminar vertical tube flow, J. Fluid Mech. 22 (4) (1965) 721–735, https://doi.org/10.1017/S0022112065001106. E.S. Asmolov, A.L. Dubov, T.V. Nizkaya, J. Harting, O.I. Vinogradova, Inertial focusing of finite-size particles in microchannels, J. Fluid Mech. 840 (2018) 613–630, https://doi.org/10.1017/jfm.2018.95. A. Yazdani, M. Deng, B. Caswell, G.E. Karniadakis, Flow in complex domains simulated by dissipative particle dynamics driven by geometry-specific body-forces, J. Comput. Phys. 305 (2016) 906–920, https://doi.org/10.1016/j.jcp.2015.11.001. B. Zhou, W. Luo, J. Yang, X. Duan, Y. Wen, H. Zhou, R. Chen, B. Shan, Simulation of dispersion and alignment of carbon nanotubes in polymer flow using dissipative particle dynamics, Comput. Mater. Sci. 126 (2017) 35–42, https://doi.org/10. 1016/j.commatsci.2016.09.012. A.L. Blumers, Y.-H. Tang, Z. Li, X. Li, G.E. Karniadakis, Gpu-accelerated red blood cells simulations with transport dissipative particle dynamics, Computer Phys. Commun. 217 (2017) 171–179, https://doi.org/10.1016/j.cpc.2017.03.016. Y. Huang, R.L. Marson, R.G. Larson, Inertial migration of a rigid sphere in plane Poiseuille flow as a test of dissipative particle dynamics simulations, J. Chem. Phys. 149 (16) (2018) 164912, , https://doi.org/10.1063/1.5047923. R.L. Marson, Y. Huang, M. Huang, T. Fu, R.G. Larson, Inertio-capillary crossstreamline drift of droplets in Poiseuille flow using dissipative particle dynamics simulations, Soft Matter 14 (12) (2018) 2267–2280, https://doi.org/10.1039/ c7sm02294h. Y. Huang, R.L. Marson, R.G. Larson, Inertial migration of neutrally buoyant prolate and oblate spheroids in plane poiseuille flow using dissipative particle dynamics simulations, Comput. Mater. Sci. 162 (2019) 178–185, https://doi.org/10.1016/j. commatsci.2019.02.048.
[42] P. Hoogerbrugge, J. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, EPL (Europhysics Letters) 19 (3) (1992) 155, https://doi.org/10.1209/0295-5075/19/3/001. [43] M. Liu, G. Liu, L. Zhou, J. Chang, Dissipative particle dynamics (dpd): an overview and recent developments, Arch. Comput. Methods Eng. 22 (4) (2015) 529–556, https://doi.org/10.1007/s11831-014-9124-x. [44] E. Pep, P. Warren, Statistical mechanics of dissipative particle dynamics, EPL (Europhysics Letters) 30 (4) (1995) 191, https://doi.org/10.1209/0295-5075/30/ 4/001. [45] R.D. Groot, P.B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (11) (1997) 4423–4435, https://doi.org/10.1063/1.474784. [46] J.A. Millan, W. Jiang, M. Laradji, Y. Wang, Pressure driven flow of polymer solutions in nanoscale slit pores, J. Chem. Phys. 126 (12) (2007) 03B617, https://doi. org/10.1063/1.2711435. [47] C.A. Marsh, J.M. Yeomans, Dissipative particle dynamics: the equilibrium for finite time steps, Europhys. Lett. (EPL) 37 (8) (1997) 511–516, https://doi.org/10.1209/ epl/i1997-00183-2. [48] D.C. Visser, H.C.J. Hoefsloot, P.D. Iedema, Modelling multi-viscosity systems with dissipative particle dynamics, J. Comput. Phys. 214 (2) (2006) 491–504, https:// doi.org/10.1016/j.jcp.2005.09.022. [49] J.A. Anderson, C.D. Lorenz, A. Travesset, General purpose molecular dynamics simulations fully implemented on graphics processing units, J. Comput. Phys. 227 (10) (2008) 5342–5359, https://doi.org/10.1016/j.jcp.2008.01.047. [50] J. Glaser, T.D. Nguyen, J.A. Anderson, P. Lui, F. Spiga, J.A. Millan, D.C. Morse, S.C. Glotzer, Strong scaling of general-purpose molecular dynamics simulations on GPUs, Comput. Phys. Commun. 192 (2015) 97–107, https://doi.org/10.1016/j.cpc. 2015.02.028. [51] H. Tollert, Die wirkung der magnus-kraft in laminaren strömungen. i. teil. der quertrieb in sedimentierenden teilchenschüttungen, Chem. Ing. Tech. 26(3) (1954,) 141–150, https://doi.org/10.1002/cite.330260306. [52] Y.W. Kim, J.Y. Yoo, Axisymmetric flow focusing of particles in a single microchannel, Lab Chip 9 (8) (2009) 1043–1045, https://doi.org/10.1039/B815286A. [53] C. Prohm, H. Stark, Feedback control of inertial microfluidics using axial control forces, Lab Chip 14 (12) (2014) 2115–2123, https://doi.org/10.1039/ C4LC00145A.
10
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Slip-shear and inertial migration of finite-size spheres in plane Poiseuille flow
T
Yuanding Huang, Xuezeng Zhao, Yunlu Pan
⁎
Key laboratory of Micro-Systems and Micro-Structures Manufacturing Ministry of Education and School of Mechatronics Engineering, Harbin Institute of Technology, Harbin, Heilongjiang Province 150001, PR China
ARTICLE INFO
ABSTRACT
Keywords: Inertial migration Sphere Slip-shear Rotation Plane Poiseuille flow
The effect of the slip-shear as well as the rotation on the inertial migration of finite-size spheres with blockage ratio (ratio of the sphere diameter to the channel gap) 0.2 are investigated in a plane Poiseuille flow at Re = 100–500 using dissipative particle dynamics simulations, by means of changing the driving force, the sphere density, and the fluid-sphere boundary slip condition, respectively. Neutrally buoyant spheres with equal density with the fluid and with no-slip boundary condition are observed to equilibrated at lateral positions around midway between the channel centerline and the wall, and the positions move closer to the channel centerline with increasing Re. The spheres always lag behind the surrounding fluid. A correction to the slip is proposed by considering the effect of the distance between the sphere and the wall. A forward driving force or a larger density facilitates the forward motion of the spheres and leads to their leading relative to the surrounding fluid, which drives the spheres to lateral equilibrium positions close to the wall; and vice versa. All spheres rotate with angular velocities approximate to but a bit less than the half of the local fluid shear rate at the sphere centroid, no mater the spheres are neutrally or non-neutrally buoyant, expect the spheres much close to the wall. The lift induced by rotation is found one to two orders of magnitude less than the lift induced by slip-shear due to the small angular slip velocity between the sphere and the surrounding fluid; while the slip-shear-induced lift is found critical to determine the lateral equilibrium positions of finite-size spheres cooperating with the wallinduced lift and the shear-gradient-induced lift as Re varies. An empirical theory considering the last three lifts above mentioned is built, which interprets the relation between the slip-shear and the equilibrium positions of the spheres. Above relation holds for spheres with no-slip boundary condition, whilst for neutrally buoyant spheres with full-slip boundary condition, though they always lag behind the surrounding fluid, they migrate towards the wall as Re increases, as their hydrodynamic interactions with the fluid are changed.
1. Introduction
throughput cytometry, enrichment of target cells, and detection of infected or diseased cells (e.g., mutant red blood cells, cancer cells, etc.) [1,7,8]. For inertial migration, neutrally buoyant solid spheres flowing in a straight channel or tube Poiseuille flow, is the most common and to some extent, the simplest case, as it avoids influences such as the noncentrosymmetric particle shape, deformation of deformable particles, viscoelasticity of vesicles, cells, or non-Newtonian medium fluid, and second flow induced by channel curvature or disturbance structure, etc. [5]. However, even for this simple case, complicated phenomena have been observed by substantial subsequent studies after the pioneering work of Segré and Silberberg [6]. It was found that rigid spherical particles are effectively focused when fluid inertia is in the intermediate 1) and the range (~1 < Re < ~100 ) between the Stokes regime (Re
Neutrally buoyant (external force- and torque-free) solid particles, droplets, vesicles, and cells, can migrate laterally across streamlines and even be focused at well-defined off-center positions when flowing with the medium fluid in a channel or a tube [1–5]. This counterintuitive phenomenon was first reported by Segré and Silberberg [6] in 1961, who observed that initially uniformly dispersed rigid spheres with millimeter size in a tube Poiseuille flow concentrated into an annulus at a position around 60% of the tube radius away from the tube centerline. They attributed it to the inertia of the moving fluid, thereafter this phenomenon is referred to as inertial migration or inertial focusing. Microfluidic devices utilizing inertial migration have been used to focus, separate, and sort particles, enabling applications such as high-
⁎
Corresponding author. E-mail address: [email protected] (Y. Pan).
https://doi.org/10.1016/j.commatsci.2020.109542 Received 18 November 2019; Received in revised form 13 January 2020; Accepted 15 January 2020 0927-0256/ © 2020 Elsevier B.V. All rights reserved.
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
turbulent regime (Re~2000 ) [5], where Re = f U m Lc / µ f is the channel Reynolds number, with f being the density of the fluid, Um being the maximum fluid velocity, Lc being the channel characteristic dimension, or rather, the inner diameter of a tube or the hydraulic diameter of a square or rectangular channel, µ f being the dynamic viscosity of the fluid. Rigid spheres are gathered into an annulus in a tube [9,10], while in a square channel, they tend to be focused at the four centers of the four faces [11–14], and at the four corners at higher Re [13,15,16]. The equilibrium positions in rectangular channels generally appear at the centers of the two long faces [16–22], and other positions at the centers of the short faces [18,21,22] or alongside the long faces [20] can emerge depending on Re and the aspect ratio of the channel cross section. The positions move closer to the centerline with increasing blockage ratio due to the size confinement, where is the ratio of the sphere diameter d to the tube diameter or to the smallest gap H of a 0.2 , their rectangular channel. For spheres with blockage ratio equilibrium positions are around 60% of the distance between the tube/ channel centerline and the wall [6,9–15,17,18,22–29]. When Re 100 , the positions move closer to the wall with increasing Re as the increasing flow inertia produces a stronger wall-directed shear-gradientinduced lift force [24,25]. However, when Re exceeds a critical value on the order of hundreds, with smaller values for higher blockage ratio , experimental [9,13,18] and numerical [14,22,28] studies show that the spheres begin to shift back towards the centerline, and concentrate to the so-called ’inner equilibrium position’ [9]. This reversal is surprising for it contradicts with predictions of the traditional asymptotic theories [9,24,25], who asserted that the spheres always migrate towards the wall with increasing Re for Re 3000 . Matas et al. [9] tried to explain this reversal through asymptotic theory but failed. Chun et al. [15] simply attribute it to the multibody interactions. Thus far no convincing explanations for this phenomenon have been given yet. Accordingly, we attempt to explore its underlying causes in this study. Theoretical study based on the matched asymptotic expansions by Ho and Leal [23] claimed that the inertial migration of rigid spheres in a Poiseuille flow results from the balance between a wall-directed shear-gradient-induced lift force Fs due to the parabolic velocity profile and a center-directed wall-induced lift force Fw due to the wall repulsion, and the lift was given as
FL = Fw + Fs =
2 4 f Um d 16H 2
( *2G1 +
G2),
s = p + 2 c is the sphere angular slip velocity, i.e. the difference between the sphere angular velocity p and the local fluid rotation 1 ( 2 c ), with c = 4Um z c / H being the flow shear rate at the sphere centroid. For a neutrally buoyant sphere in Poiseuille flow, it was found that the sphere always lags behind its surround fluid [23,31], and the sphere angular velocity is kind of less than the local fluid rotation [32], thus the lift induced by rotation directs towards the wall according to Eq. (2). However, no direct evidence supports this. Recently, Joseph and Ocando [32] claimed that the sphere angular slip velocity s is critical to the lift, which changes sign above and below the equilibrium position in Poiseuille flow, through their direct numerical simulations. Saffman [33] theoretically studied the effect of the slip-shear on small spheres in a linear shear flow, and an expression for the slipshear-induced lift force was given as 1
Fss = 1.615µ f Us d 2
8
fd
3
× Us,
(3)
where f = µ f / f is the kinematic viscosity of the fluid. He found that the lift induced by rotation [Eq. (2)] is smaller by one order of magnitude than that induced by slip-shear [Eq. (3)] for small Re s . Eq. (3) was then extended to the case in Poiseuille flow under conditions that Re 1, Re p = Re 2 1, and the sphere is not too close to either the tube axis or the wall. According to Eq. (3), the lift directs to the centerline when the sphere lags behind its ambient undisturbed fluid and directs to the wall when the sphere leads its ambient undisturbed fluid, which is thus known as Saffman effect. Above theory is qualitatively in agreement with experimental observations by Jeffrey and Pearson [34], who observed that, in a vertical tube, the spheres slightly denser than the fluid, i.e. with a slip relative to the surrounding fluid, drifted towards the wall for downward flow and moved towards the centerline for upward flow, though the spheres are with finite size ( = 0.05–0.1). Recenty, Asmolov et al. [35] generalized a theory to finite-size spheres at Re 20 , and found the slip-shear is crucial to the lift and determines the equilibrium positions of the spheres. Their theory is then validated 0.1 and 0.3. by lattice Boltzmann simulations for spheres with Though the effects of the slip-shear and the rotation on the inertial migration of neutrally buoyant spheres in Poiseuille flow have attracted some attentions in the past years, some issues are still mazy and poorly elucidated, especially their effects on the inward migration of the spheres with finite size at high Re. Early theoretical studies [9,23–25,29] cannot cope with the fluid disturbance triggered by the sphere with finite size. In experiments, it is difficult or even impossible to accurately measure the slip and the rotation velocity, especially at high Re. Thanks to the advance of computing capabilities, numerical simulation has become a powerful way to investigate hydrodynamic topics with well-controlled parameters. Here we employ a particlebased simulation method, dissipative particle dynamics (DPD), to investigate the inertial migration of neutrally buoyant spheres, with particular consideration of the effects of the slip-shear and the rotation on their equilibrium positions. DPD is a versatile method well applicable to the simulations of simple and complex fluids at mesoscopic scale [36–38], and it has been verified to be capable of the simulation of inertial migration of particles in Poiseuille flow up to Reynolds number of 500 or so by our previous studies [39–41]. We will introduce it as well as our model in the next section, following by our results and discussion, and ending with our concluding remarks.
(1)
= 4 are dimensionless shear rate and shear = 4z and where gradient normalized by Um and H, respectively, with z = 2z /H being the relative lateral position with respect to the half of the channel gap and z being the lateral coordinate; G1 and G2 are functions of z for wallinduced lift and shear-gradient-induced lift, respectively. Noteworthily, asymptotic theories [9,23–25] are based on the point-particle assumption, and the flows are not disturbed by the particle. It has shown that the asymptotic theories [9,23–25] cannot account for the inward migration of spheres at high Re, and other factors should have effects. Slip-shear-induced lift force and rotation-induced lift force, are two other lift forces exerted on the particles when they flow with the fluid. These two lifts should be significant at high Re [5], however, they used to be deemed as negligible relative to the shear-gradient-induced and wall-induced lifts [23], and are less considered, especially for spheres with finite size flowing in Poiseuille flows. Rubinow and Keller [30] calculated the lift induced by rotation in linear shear flow by involving the Stokes and Oseen expansions for small but finite Re s , where Re s = f |U s|d/ µ f , with Us = Up Uf being the relative velocity between the sphere and the fluid at the sphere centroid, which was given as Fr =
/ f,
2. Numerical algorithm and model 2.1. Numerical algorithm DPD was originally proposed by Hoogerbrugge and Koelman [42] in 1992. In DPD, fluid is composed of discrete beads, with each bead coarse-grained of a group of fluid molecules, which allows larger temporal and spatial scales to be accessed relative to traditional Molecular Dynamics simulations [43]. The motion of each bead is
(2)
where = p , with p being the rotational angular velocity of the = s [5], where sphere. Subsequently, it was corrected that 2
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
governed by Newton’s Second Law
d ri dv = vi , mi i = dt dt
Fij + F iext ,
(4)
where bead i located at position ri and with mass mi and velocity vi , is subjected to the pairwise additive inter-particle force Fij from bead j and external force Fiext . Fij consists of three center-to-center directed forces, namely a conservative force FCij , a dissipative force FD ij , and a random force FRij :
Fij = FCij + FijD + FRij .
(5)
The conservative force is set as a soft repulsion, which allows a large time step in DPD simulations, and is given as
FCij =
Aij (1 0,
rij ) rij, rij rc rij > rc
(6)
Aij where is the bead-bead repulsion coefficient, rij = ri rj, rij = rij/|rij |, and rc is the cutoff radius for inter-particle interactions. Aij is related to the compressibility of the DPD fluid, and in multiphase systems can be used to differentiate fluids of different kinds from each other. A dissipative force is introduced to transfer momentum between beads more efficiently than is possible from the soft potential alone, while a random force is applied to compensate for the loss of kinetic energy resulting from the dissipation. The dissipative force and the random force are given by FD ij =
D (r )( r · v ) r , ij ij ij ij
FRij =
R (r ) ij ij
t
1 2 rij,
Fig. 1. Model for a rigid sphere (in yellow) flowing in plane Poiseuille flow (in blue) with rigid walls (in gray).
are densely packed with a density w = 61.35rc 3 in HCP lattice to construct the wall, which is placed at z = 10rc and parallel to the x-y plane; see Fig. 1. The fluid beads are randomly distributed within the box with a density f = 4rc initially. The high bead density in the wall can effectively restrain the penetration of the fluid beads across the wall, although its thickness is small. In addition, the periodic boundary condition of the box enables the single wall to model a pair of walls perpendicular to the z direction, giving the height of the channel as H = Lz = 40rc . Wall beads are frozen in position during the course of the simulations. Multiple DPD beads packed in FCC lattice in spherical shape are connected by harmonic bonds to construct the rigid sphere. Considering the large size (compared to the effective size of a fluid bead) of the sphere in three dimensions, a lower number density n p = 32rc 3 is chosen for the sphere than that of the wall, although this density is still high enough to effectively prohibit the penetration of fluid beads into the sphere, but not so high as to make computer memory cost unmanageable. Each bead in the sphere is connected with its nearest, next-nearest, and next-next-nearest beads through a bond with a potential energy
(7)
where is the friction coefficient, is the amplitude of noise, t is the integral time step, D (rij ) and R (rij ) are rij -dependent dimensionless weighting functions, and ij is the Gaussian white noise following Gaussian statistics: and ij (t ) = 0 ij (t ) kl (t ) = ( ik jl + il jk ) (t t ) . The fluctuation–dissipation theorem linking the dissipative and random forces gives D (r )
=[
D (r )]2
and
2
= 2 kB T
(8)
as explained by Español and Warren [44]. Thus, the momentum and the temperature are conserved. D (r ) and R (r ) conventionally take the form D (r )
=[
D (r )]2
(1
=
r ) 2 , r rc 0, r > rc
(9)
An unmodified version of velocity-Verlet algorithm is used to advance the set of positions and velocities in the DPD we employed, which is given as
ri (t + t ) = ri (t ) + t vi (t ) + vi (t + t ) = vi +
1 2
1 ( 2
V (r ) =
t Fi , 1 2
t (Fi (t ) + Fi (t + t ).
l 0)2
(11)
where kb is the force constant and is set with a high value kb = 1000m p to restrain the flexibility of the bond, with m p being the mass of the sphere bead, l is the length of the bond, and l 0 is the bond rest length, with l 0 = 2 b/2, b, 3 b for bonds connecting the nearest, next-nearest, and next-next-nearest beads, respectively, with b = 0.5rc . The sphere is initially placed at the position (0, 0, 10rc) , which is at the center of the channel. Following the simulation work for Poiseuille flow in a planar channel using DPD by Millan et al. [46], we adopt Aij values that they took for our purposes, which are
t )2Fi (t ),
Fi (t + t ) = Fi (r (t + t ), v (t + t )), vi (t + t ) = vi (t ) +
1 k b (l 2
(10)
In a ‘standard’ DPD simulation [45], m , rc and kB T are normalized to unity, i.e., m = rc = kB T = 1; the unit of time is defined as = rc m / kB T = 1; the density of the fluid is taken as f = 4rc 3; and the friction coefficient between the fluid beads is set as ff = 4.5 in DPD units. 2.2. Model
w f p 3 10 kB T w Aij = f 3 25 3 , rc p 10 3
A box, with size L x × L y × Lz = 60rc × 40rc × 40rc and periodic boundary condition in all directions, is employed. The length (L x ) has been verified to be long enough for the fully-developed plane Poiseuille flow, and the width (L y ) is wide enough to minimize interactions between periodic images of the spheres [39]. Three layers of DPD beads
where the subscripts “w”, “f”, and “p” represent the DPD beads forming the wall, the fluid, and the spherical particle, respectively. The fluid–fluid repulsion parameter Aff is not far from the standard DPD value of 18.75 [45], while the other parameters (along with the choice of densities discussed above) are chosen to minimize fluid penetration and 3
(12)
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
N gives the standard error used to obtain the error bars reported here in our figures. Simulations are performed with DPD in the Highly Optimized Object-oriented Molecular Dynamics (HOOMD)-Blue software package [49,50], a GPU-based simulation toolkit. 3. Results In our simulations, the spheres are chosen to have a diameter d = 8rc , and the blockage ratio is therefore = d/ H = 0.2 . A smaller sphere is found to suffer significant Brownian motion and is hard to be focused. A larger sphere is beyond the computational ability when analyzing the data due to the overmuch beads within it. A “standard” case for the inertial migration of neutrally buoyant sphere in plane Poiseuille flow has been simulated in our previous work [39] as an effort to assess the DPD fluid. In that case, the bead mass of the sphere is one-eighth of that of the fluid (i.e., m p = m f /8 = 1/8), as is the driving force applied on per sphere bead compared to that of the fluid (i.e., fx p = fx f /8), so that equal density (note the bead number density is 8fold higher than that of fluid, i.e., n p = 8n f ) and acceleration (induced by the driving force) are achieved between the fluid and the sphere. The friction coefficient between the sphere beads and the fluid beads is set as fp = 4.5, thus there is no-slip at the fluid-sphere interface. This corresponds to the common case of the motion of a neutrally buoyant sphere free drifting in Poiseuille flow without external force or torque. The equilibrium positions of above sphere are around at middle between the channel centerline and the wall, and they move closer to the centerline as Re increases, which are shown red symbols in Fig. 3. Our results have been compared with other experimental and numerical studies for spheres with the same (=0.2) in square [22,14], rectangular [17,18,22], and plane [22] channels, respectively. Our results are in agreement with other studies for spheres in rectangular channels [17,22] and a plane channel [22] at Re = 100–200 . Based on the results
Fig. 2. Example velocity profiles of pure DPD fluid with fw = 1.0 (blue) and 4.5 (red), respectively, and are compared with the Hagen-Poiseuille equation (the black solid curve).
slip against the spherical particle and the wall. The friction coefficient between the fluid beads and the wall beads is set as fw = 4.5, with which the no-slip boundary condition is achieved at the fluid-wall interface. One example of the velocity profile from our simulation of the Poiseuille flow of pure fluid in the plane channel is shown with red symbols in Fig. 2, which is compared to the Hagen-Poiseuille equation (the black curve in Fig. 2)
u x (z ) =
1 dp 2 z 2µ f dx
H 2
2
=
f fx f
2µ f
z2
H 2
2
,
(13)
where u x (z ) is the flow velocity in the x direction (flow direction) as a function of gap coordinate z , dp /dx = f fx f is the pressure gradient in the x direction, and fxf is the external driving force on per fluid bead along the x direction. The external driving forces in y and z directions are set to zero, i.e., F ext f = [fx f , f y f , fz f ] = [fx f , 0, 0]. The good match of our simulation results with the Hagen-Poiseuille equation demonstrates that our parameter choice yields the Hagen-Poiseuille flow with a noslip boundary condition, and thus, the viscosity of the fluid can be estimated by
µf =
f fx f H
8Um
2
.
|z*eq|
0.6
0.4 Our work (nmrcl, plane channel) Lashgari et al. (nmrcl, square channel) Yuan et al. (nmrcl, square channel) Masaeli et al. (expt, rectangular channel, AR ≈ 1.6) Ciftlik et al. (expt, rectangular channel, AR = 1.6) Lashgari et al. (nmrcl, rectangular channel, AR = 2) Lashgari et al. (nmrcl, plane channel)
(14)
Decreasing the fluid-wall friction coefficient fw to a smaller value will lead to a partial-slip at the fluid-wall interface (one example of velocity profile of flow with fw = 1.0 is shown with blue symbols in Fig. 1) and to zero will result in a full-slip (no example of velocity profile can be shown as the fluid will be accelerated all the time, but the profile keeps flat). Similarly, the fluid-sphere boundary condition can be adjusted through the friction coefficient fp . 0.01, and it satisfies t = 0.1 / to keep the The time step is t/ fluctuation of the system temperature within 2% [54,48]. The system is first equilibrated by running for one million time steps to establish an accurate model. After applying driving forces, all simulations are run for one million steps to achieve a nearly steady-state position of the particle, and then for three million additional steps to collect fully developed data. After discarding the initial part of the run, the fluctuations around the equilibrium position are autocorrelated to obtain the characteristic autocorrelation time. The length of the run divided by this autocorrelation time gives the effective number of statistically independent positions N sampled by the sphere during the course of the simulation. Dividing the standard deviation of the average position by
0.2
0.0
0
100
200
300 Re
400
500
600
Fig. 3. Relative equilibrium positions z eq of neutrally buoyant spheres with p = f and no-slip boundary condition (the standard case) in plane channel
through our DPD simulations, comparing results of other studies for neutrally buoyant spheres with the same blockage ratio ( = 0.2 ) in square (Lashgari et al. [22], Yuan et al. [14]), rectangular (Masaeli et al. [17], Ciftlik et al. [18], Lashgari et al. [22]), and plane (Lashgari et al. [22]) channels, respectively. Notably, Re in some publications have different definition, which we have transformed to our definition which is based on the maximum velocity Um and the smallest gap of rectangular channel H.
4
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
0.8
0.8
(a)
axp/axf = 8 axp/axf = 4 axp/axf = 1
200
0.10
300 Re
0.4
ρp/ρf = 4
400
0.0
500
100
200
300 Re
Us [rc/τ]
Faxen correction
0.05
0.00
0.00
-0.10
-0.15
-0.15
100
200
300 Re
0.10
-0.15
0.04
0.02
0.02 300 Re
400
100
200
300 Re
400
500
0.12 (f)
0.10
500
0.00
(i)
0.08
0.06
0.04
200
500
0.08
αc/2
100
400
0.12
(c)
0.06
0.00
500
500
ωp [rad/τ]
ωp [rad/τ]
400
400
-0.10
ωp [rad/τ]
0.12
300 Re
300 Re
-0.05
-0.20
200
200
(h)
0.05
-0.10
100
0.10
-0.05
-0.05
0.08
0.0
500
(e)
0.00
0.10
400
0.10
0.05
100
Full slip Partical slip No-slip
0.2
ρp/ρf = 1
(b)
-0.25
0.4
ρp/ρf = 8 0.2
Us [rc/τ]
100
axp = 0 axp/axf = -4
0.6
|z*eq|
|z*eq|
|z*eq|
0.4
0.2
Us [rc/τ]
(g)
0.6
0.6
0.0
0.8 (d)
0.06 0.04
αc/2 100
200
300 Re
400
αc/2
0.02 500
0.00
100
200
300 Re
400
500
Fig. 4. Equilibrium positions [(a), (d), (g)], slips [(b), (e), (h)], and rotational angular velocities [(c), (f), (i)] of spheres with = 0.2 in plane Poiseuille flow. (a)–(c) are for spheres with different driving forces, (d)–(f) are for spheres with different densities, and (g)–(i) are for spheres with different fluid-sphere boundary slip conditions. The dashed lines in (c), (f), and (i) are Tollert’s theory [51] that p = c /2 .
from other studies shown in Fig. 3, the equilibrium positions move closer to the centerline with increasing AR, where AR is the aspect ratio of the cross section of the rectangular channel, which for square channels AR = 1 and for plane channels AR = . It is therefore reasonable that the positions for our spheres in plane channel are a bit closer to the centerline. Also, the trends of the inward migration are similar to that in square [14] and rectangular [18] channels. The validity of our model is thus verified. Here, the slip and rotational angular velocity of above spheres are measured and shown with red symbols in Fig. 4 (b) and (c), respectively. The slip is an average of 100 data samples when the sphere has been focused around its equilibrium position; see Fig. A.1 in Appendix
A. It is found that the spheres always lag behind their surrounding fluid, and the lag increases with increasing Re. The magnitudes of the slips are within 4% of the velocity of the ambient fluid, and 1.2–1.8 times of 2U /3 [35], with smaller times for sphere Faxen correction Us = m farther away from the wall. Asmolov et al. [35] also observed slips about twice larger than the Faxen correction for finite-size spheres 0.1) which are not too close to the wall. The deviations from the ( Faxen correction attributes to the resistance effect from the wall, which diminishes as the sphere moves away from the wall. As the slip is related to the distance between the sphere and the wall, the slip for finitesize spheres which are not too close to the wall in plane Poiseuille flow 0.27 for the spheres can be fitted as Us = k s 2Um/(1 z ) , where k s 5
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
of the standard case based on our simulation data. The rotational angular velocity is measured by tracing two beads at the two poles of one diameter of the sphere along 2 × 103 after the sphere being well focused, and one example of rotating orbit in the shear-gradient plane (x-z plane) is shown in Fig. A.2 in Appendix A. The angular velocity of the sphere p [the red curve in Fig. 4 (c)] approximately increases linearly with increasing Re and p c/2 , where c is the local flow shear rate at the sphere centroid, in accord with Tollert’s theory [51] for a free moving sphere in linear shear flow. This should attribute to the intermediate equilibrium position of the sphere between the centerline and the wall, at which the Poiseuille flow velocity profile can be approximately regarded as linear, and the wall resistance effect is slight. Also, we note that the local fluid rotation ( c /2 ) is a bit greater than the sphere angular velocity p , in agreement with the finding by Joseph and Ocando [32]. To further investigate the effects of the slip-shear, three different means, i.e., varying the driving force applied on the spheres, varying the sphere density, and changing the fluid-sphere boundary slip condition, are employed, respectively; while other parameters are kept the same as the standard case. Four different driving forces, which are fx p = fx f , fx p = fx f /2, fx p = 0 , and fx p = fx f /2, are applied on per bead of the sphere, respectively, while keeping p = f and fp = 4.5. The acceleration induced by the driving force on the sphere axp are therefore 8, 4, 0, and 4 times that on the fluid ax f . It can be seen from Fig. 4 (b) that, for ax p /a x f = 8 and 4, the sphere leads its surrounding fluid, and the leading increases as Re and axp increase. For ax p /a x f = 1, 0, and 4 , the sphere lags behind its surrounding fluid, and the lag increases with increasing Re but decreasing axp . The equilibrium positions of the spheres are close to the wall when the spheres lead the surrounding fluid, and close to the centerline when the spheres lag behind the surrounding fluid, in accordance with Saffman’s prediction [33]. The trend is more pronounced as the absolute value of the slip |Us | becomes larger. For ax p /a x f = 4 , the spheres are focused to positions around the channel centerline. This is in agreement with recent studies which, by applying an external axial force on the rigid spheres in the opposite direction to induce a large lag relative to the flow, achieved the focusing of the spheres with finite size at the channel center in Stokes regime (Re 0.05) experimentally [52] and at moderate Re (10 80 ) numerically [53], respectively. The rotational angular velocity still linearly increases with increasing Re, and approximately follows p = c/2 , especially for spheres that are away from both the wall and the channel centerline. For spheres close to the wall, their rotation rates deviate from c /2 due to the growing wall resistance; and for spheres at the channel centerline, their rotations are restrained and cannot be measured accurately as the local shear rate approaches 0. The variation of the particle density is realized by changing the mass of the sphere bead simultaneously with the driving force, in which way keeps the acceleration of the sphere induced by the driving force equal to that of the fluid. Notably, change of the mass of the sphere bead leads to unequal density between the sphere and the fluid, i.e., non-neutrally buoyant for spheres in real world. However, no gravity or buoyancy is applied in our simulation system, and therefore the spheres are still neutrally buoyant. Equilibrium positions, slips, and rotational angular velocities of spheres with densities 4 and 8 times that of fluid are shown in Fig. 4 (d)–(f) and compared with the standard case, respectively. It is shown that the increase of particle density (or inertia) accelerates the sphere, as for the spheres with p / f = 4 , though they still slightly lag behind its surrounding fluid, their lags become smaller relative to the standard case ( p / f = 1), and the spheres with p / f = 8 become to lead the surrounding fluid. The relation between the equilibrium positions and the slips still follows Saffman’s prediction [33]. It should be noted that the inward migration of the sphere with p / f = 8 at high Re (400–500) is due to the increasing fluctuation of the sphere in lateral direction, which leads to a small average of the position; see Fig. A.3 in Appendix A. The rotational angular velocity follows p = c/2 with
exceptions when the spheres are close to the wall. The slip between the sphere and its surrounding fluid is also altered by changing the fluid-sphere boundary slip condition, specifically, changing the fluid-sphere friction coefficient fp . We have shown that decreasing fw to 1 brings off a partial-slip at the fluid-wall interface, and to 0 leads to a full-slip. Similarly, we set fp = 1 and 0 to implement the partial-slip and the full-slip at the fluid-sphere interface, respectively; and then compare the results with the standard case which with no-slip fluid-sphere boundary condition ( fp = 4.5); see Fig. 3. (g) - (i). It shows that all the spheres lag behind their surrounding fluid, and the lags of the spheres with partial- and full-slip are almost the same, and are kind of smaller than that with no-slip at the same Re. It seems the boundary slip condition has slight effect on the slip of the sphere while the particle inertia play a decisive role. However, the equilibrium positions of the spheres with full-slip move close to the wall with increasing Re, which differ significantly from the trend for spheres with no-slip, while the equilibrium positions of spheres with partial-slip are kind of closer to wall than that of no-slip at the same Re. The rotational angular velocities of the spheres with full-slip deviate from p = c/2 , though the deviations are small. At low Re, where the sphere is relatively away from the wall, the full-slip boundary condition reduces the friction between the sphere and the fluid, the rotation rate is kind of larger than c /2 . At high Re, where the sphere is close to the wall, the wall resistance restrains the rotation and the angular velocity is kind of smaller than c /2 . 4. Discussion Above results show that, for spheres with no-slip boundary condition, a larger forward driving force or a larger density (or inertia) on/of the spheres than that on/of the fluid is conductive to their forward motion, and leads to the leading relative to their surrounding fluid; while a backward driving force make the spheres lag behind the surrounding fluid. As Saffman’s theory [33] predicted, the lift induced by the slip-shear directs to the wall when the sphere leads its surrounding fluid and directs to the centerline when it lags behind its surrounding fluid. Our simulation results quantitatively agree with this. The sphere angular velocity p c/2 , and is kind of less than the local fluid rotation ( c /2 ), except when the sphere is close to the wall, in accordance with the observations by Joseph and Ocando [32]. Due to the sphere angular slip velocity s is small, the rotation-induced lift [according to Eq. (2)] is found one to two orders of magnitude less than the lift induced by slip-shear [according to Eq. (3)] based on our simulation data, as shown in Fig. 5. This also holds for the spheres much close to the wall. Thus the rotation-induced lift can be neglected. Meanwhile, the slip-shear-induced lift is found comparable to the wall-induced lift and the shear-gradient-induced lift, which should play a role in the inertial migration of the sphere. The non-dimensionalized lift arising from the competition between the wall-induced lift and shear-gradient-induced lift can be written as
Fw + Fs =
Fw + Fs = G1 z *2 + G2 z 2 4 2 f Um d /H
(15)
according to Eq. (1), which is shown with black curve in Fig. 6, re0.6 vealing that the equilibrium position of the sphere is fixed at z and has no relation with the Um , i.e., Re. However, our study and most other studies have shown that the position varies with Re. Whilst the slip-shear-induced lift proposed by Saffman [33] can be non-dimensionalized as
Fss =
Fss ks = 3.23 2 4 2 1 z U d / H m f
z Re
(16)
by substituting Us = k s 2Um/(1 z ) and = 4Um z / H into Eq. (3), which indicates Fss is relative to Re. Among these three dominant lifts on spheres, the slip-shear-induced lift is therefore critical to adjust the 6
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
spheres in plane Poiseuille flow at Re over the range of 100
30
Fw Fs Fss Fr
FL
20
FL = FL/(
0 -10 200
300 Re
400
500
Fig. 5. The lifts induced by wall, shear-gradient, slip-shear, and rotation, respectively, on a neutrally buoyant sphere with p = f and no-slip boundary condition (i.e., the standard case) in plane Poiseuille flow.
0.2
FL/( fU2md4/H2)
0.0 -0.2 -0.4
Ho and Leal Re = 100 Re = 200 Re = 300 Re = 400 Re = 500
-0.6 -0.8 -1.0 0.0
0.2
0.4
z*
= Fs + Fw + css (Re) Fss.
500 as (17)
css is chosen in the form of css = g Re + h for simplicity, where g and h are a coefficient and a constant, respectively, and are set by fitting our simulation results into Eq. (17) through the least square method. css (Re) 0.0073Re 0.255 for the neutrally buoyant sphere with p = f and no-slip boundary condition (i.e., the standard case), and the non-dimensionalized lifts FL are shown with colored curves in Fig. 6 for different Re. The inward migration of neutrally buoyant spheres with p = f and no-slip boundary condition in Poiseuille flow at Re = 100 500 is thus interpreted. css as well as ks for other cases are shown in Table B.1 in Appendix B, which account for the Re-dependent inward migration for spheres lagging behind the surrounding fluid and the Re-dependent outward migration for spheres leading the surrounding fluid through Eq. (17). For non-neutrally buoyant spheres, by applying an forward or backward axial driving force, the magnitude and the direction of Us changes, resulting in a redistribution of the equilibrium positions from the channel centerline (z eq = 0 ) to the wall (z eq = 0.8); see Fig. 4 (a) and (b), respectively. This suggests a promising application, which utilizes an external driving force, such as electric, magnetic, or gravitational force, to adjust the sphere equilibrium positions within the channel and then for particle separation. Actually, Kim and Yoo [52] have achieved the focusing of rigid spheres at the channel centerline by applying a backward axial electric force on the spheres, however, at Stokes regime. One intractable limitation for this application is the strong hydrodynamic force on the sphere at high Re, which requires a strong external force to compare or overwhelm it. Changing the density of the spheres changes the particle inertia, which in turn changes the slips and the equilibrium positions of the spheres. However, this needs the change of the sort of the particles in reality and will bring about sedimentation or floatation in horizontal channels. By contrast, altering the fluid-sphere boundary slip condition brings more inspiration to the potential application for inertial migration/focusing at high Re. According to our simulation results, for a sphere with full-slip boundary condition, though it lags behind its surrounding fluid, its equilibrium position moves closer to the wall as Re increases. The change of boundary slip condition changes the hydrodynamic interaction between the sphere and its surrounding fluid, in consequence, above theories built for spheres with no-slip boundary condition do not hold here anymore.
10
100
2 4 2 f Um d / H )
5. Conclusion
0.6
0.8
We have performed dissipative particle dynamics simulations to study the effect of the slip-shear as well as the rotation on the inertial migration of finite-size spheres with blockage ratio = 0.2 flowing in a plane Poiseuille flows at Re = 100 500 , by changing the driving force, the sphere density, and the boundary slip condition between the sphere and the fluid, respectively. For neutrally buoyant spheres with p = f and no-slip boundary condition, they are equilibrated at lateral positions around midway between the channel centerline and the wall, with the positions moving closer to the centerline with increasing Re. The spheres are found always lagging behind the surrounding fluid with a 0.27 , which is improved upon the slip Us = k s 2Um/(1 z ) with k s 2U /3 [35] by considering the classical Faxen correction that Us = m distance between the sphere and the wall. A forward driving force or a large density relative to the fluid facilitates the forward motion of the sphere and leads to its leading with respect to the surrounding fluid, resulting in an outward migration of the sphere towards the wall; while a backward driving force makes the sphere lag behind the surrounding fluid, giving rise to an inward migration of the sphere towards the channel centerline. This quantitatively agrees with Saffman’s prediction [33]. Meanwhile, the rotational angular velocity p c/2 , and is
1.0
Fig. 6. Non-dimensionalized lift for theory proposed by Ho and Leal [23] FL = Fs + Fw (the black curve), and for present theory FL = Fs + Fw + css (Re)F ss (the colored curves), where css (Re) 0.0073Re 0.255 is for neutrally buoyant spheres with p = f and noslip boundary condition.
equilibrium positions of the spheres with variation of Re. However, Eq. (16) shows that the magnitude of Fss decreases with increasing Re. According to our simulations, the impact of slip-shear-induced lift increases with increasing Re over the range of 100 500 , which leads to the progressive inward migration for spheres lagging behind the ambient fluid and the progressive outward migration for spheres leading the ambient fluid. Based on above analysis, we multiply Fss by a coefficient css , which is a function of Re, and propose a lift for finite-size 7
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
kind of less than the local fluid rotation ( c /2 ), except when the sphere is close to the wall, where the wall resistance is significant. Based on previous theories [23,30,33], we found the rotation-induced lift is less by one to two orders of magnitude than the slip-shear-induced lift, while the latter is comparable to the dominant wall-induced lift and shear-gradient-induced lift, and is critical to adjust the equilibrium positions within the channel when Re varies. A theory including the lifts induced by wall, shear-gradient, and slip-shear is built, which well explains the relation between the slip-shear and the inertial migration of spheres. Above theory holds for spheres with no-slip boundary condition. Whilst for spheres with full-slip boundary condition, as their hydrodynamic interactions with the fluid are changed, they migrate towards the wall as Re increases, though they always lag behind the surrounding fluid.
CRediT authorship contribution statement Yuanding Huang: Conceptualization, Investigation, Methodology, Software, Data curation, Formal analysis, Visualization, Writing - original draft. Xuezeng Zhao: Funding acquisition, Supervision, Writing review & editing. Yunlu Pan: Funding acquisition, Project administration, Resources, Supervision, Formal analysis, Validation, Writing review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment
Data availability
This research was supported by the National Natural Science Foundation of China (Grant Nos. 51505108 and 51475118) and Key Laboratory of Micro-systems and Micro-structures Manufacturing of Ministry of Education, Harbin Institute of Technology (No. 2017KM009).
The raw data required to reproduce these findings are available to down-load from [ https://data.mendeley.com/datasets/jcjhhps7tr/ draft?a=b00e7da0-6e90-4bfc-ae36-aeb1c6c43aa1 ].
Appendix A. Axial velocity, rotating orbit, and motion trajectories of spheres
Fig. A.1. Comparison of axial velocities of neutrally buoyant spheres which with positions with the fluid velocities at different Re.
p
=
p
and no-slip boundary condition (i.e., the standard case) at their equilibrium
Fig. A.2. Example of rotating orbit of a neutrally buoyant sphere of the standard case in the shear-gradient plane (x-z plane), where mapping of a chosen diameter in x-z plane with the x direction. 8
z
is the angle between the
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
Fig. A.3. Comparison of motion trajectories of spheres with different densities at Re
500 .
Appendix B. css (Re) and k s
Table B.1 css (Re) and k s for spheres with different driving forces and different densities over the Re range 100 Cases p
=
f
ax p = ax f
ks
css (Re) ax p/ax f =
0.0006Re
4
500 .
0.21
-1.67 -0.63
ax p/ax f = 0
0.0034Re + 0.20
ax p/ax f = 1
0.0073Re
0.26
ax p/ax f = 4
0.0027Re
0.15
ax p/ax f = 8
0.0030Re + 0.17
0.23
-0.27 0.09
p/ f
=4
0.0114Re
1.60
-0.05
p/ f
=8
0.0022Re + 1.79
0.16
[15] B. Chun, A. Ladd, Inertial migration of neutrally buoyant particles in a square duct: An investigation of multiple equilibrium positions, Phys. Fluids 18 (3) (2006) 031704, , https://doi.org/10.1063/1.2176587. [16] A.A.S. Bhagat, S.S. Kuntaegowdanahalli, I. Papautsky, Enhanced particle filtration in straight microchannels using shear-modulated inertial migration, Phys. Fluids 20 (10) (2008) 101702, , https://doi.org/10.1063/1.2998844. [17] M. Masaeli, E. Sollier, H. Amini, W. Mao, K. Camacho, N. Doshi, S. Mitragotri, A. Alexeev, D. Di Carlo, Continuous inertial focusing and separation of particles by shape, Phys. Rev. X 2 (2012) 031017, , https://doi.org/10.1103/PhysRevX.2. 031017. [18] A.T. Ciftlik, M. Ettori, M. Gijs, High throughput-per-footprint inertial focusing, Small 9 (16) (2013) 2764–2773, https://doi.org/10.1002/smll.201201770. [19] J. Zhou, I. Papautsky, Fundamentals of inertial focusing in microchannels, Lab Chip 13 (6) (2013) 1121–1132, https://doi.org/10.1039/C2LC41248A. [20] C. Liu, G. Hu, X. Jiang, J. Sun, Inertial focusing of spherical particles in rectangular microchannels over a wide range of reynolds numbers, Lab Chip 15 (4) (2015) 1168–1177, https://doi.org/10.1039/C4LC01216J. [21] C. Liu, C. Xue, J. Sun, G. Hu, A generalized formula for inertial lift on a sphere in microchannels, Lab Chip 16 (5) (2016) 884–892, https://doi.org/10.1039/ C5LC01522G. [22] I. Lashgari, M.N. Ardekani, I. Banerjee, Inertial migration of spherical and oblate particles in straight ducts, J. Fluid Mech. 819 (2017) 540–561, https://doi.org/10. 1017/jfm.2017.189. [23] B.P. Ho, L.G. Leal, Inertial migration of rigid spheres in two-dimensional unidirectional flows, J. Fluid Mech. 65 (02) (1974) 365–400, https://doi.org/10.1017/ s0022112074001431. [24] J.A. Schonberg, E.J. Hinch, Inertial migration of a sphere in Poiseuille flow, J. Fluid Mech. 203 (1989) 517–524, https://doi.org/10.1017/s0022112089001564. [25] E.S. Asmolov, The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number, J. Fluid Mech. 381 (1999) 63–87, https://doi.org/ 10.1017/s0022112098003474. [26] J. Feng, H.H. Hu, D.D. Joseph, Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 2. Couette and Poiseuille flows, J. Fluid Mech. 277 (1994) 271–301, https://doi.org/10.1017/s0022112094002764. [27] D. Di Carlo, D. Irimia, R.G. Tompkins, M. Toner, Continuous inertial focusing, ordering, and separation of particles in microchannels, Proc. Nat. Acad. Sci. 104 (48) (2007) 18892–18897, https://doi.org/10.1073/pnas.0704958104. [28] X. Shao, Z. Yu, B. Sun, Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers, Phys. Fluids 20 (10) (2008) 103307, , https://doi.org/10.1063/1.3005427. [29] J.-P. Matas, J.F. Morris, E. Guazzelli, Lateral force on a rigid sphere in large-inertia
References [1] D. Di Carlo, Inertial microfluidics, Lab Chip 9 (21) (2009) 3038–3046, https://doi. org/10.1039/b912547g. [2] P. Sajeesh, A.K. Sen, Particle separation and sorting in microfluidic devices: a review, Microfluid. Nanofluid. 17 (1) (2014) 1–52, https://doi.org/10.1007/s10404013-1291-9. [3] H. Amini, W. Lee, D. Di Carlo, Inertial microfluidic physics, Lab Chip 14 (15) (2014) 2739–2761, https://doi.org/10.1039/c4lc00128a. [4] J.M. Martel, M. Toner, Inertial focusing in microfluidics, Annu. Rev. Biomed. Eng. 16 (1) (2014) 371–396, https://doi.org/10.1146/annurev-bioeng-121813-120704. [5] J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M.E. Warkiani, W. Li, Fundamentals and applications of inertial microfluidics: a review, Lab Chip 16 (1) (2016) 10–34, https://doi.org/10.1039/c5lc01159k. [6] G. Segré, A. Silberberg, Radial particle displacements in poiseuille flow of suspensions, Nature 189 (4760) (1961) 209–210, https://doi.org/10.1038/189209a0. [7] G.Y. Lee, C.T. Lim, Biomechanics approaches to studying human diseases, Trends Biotechnol. 25 (3) (2007) 111–118, https://doi.org/10.1016/j.tibtech.2007.01. 005. [8] A.J. Mach, O.B. Adeyiga, D. Di Carlo, Microfluidic sample preparation for diagnostic cytopathology, Lab Chip 13 (6) (2013) 1011–1026, https://doi.org/10.1039/ c2lc41104k. [9] J.P. Matas, J.F. Morris, É. Guazzelli, Inertial migration of rigid spherical particles in Poiseuille flow, J. Fluid Mech. 515 (2004) 171–195, https://doi.org/10.1017/ s0022112004000254. [10] Y.-S. Choi, S.-J. Lee, Holographic analysis of three-dimensional inertial migration of spherical particles in micro-scale pipe flow, Microfluid. Nanofluid. 9 (4) (2010) 819–829, https://doi.org/10.1007/s10404-010-0601-8. [11] D. Di Carlo, J.F. Edd, K.J. Humphry, H.A. Stone, M. Toner, Particle segregation and dynamics in confined flows, Phys. Rev. Lett. 102 (2009) 094503, , https://doi.org/ 10.1103/PhysRevLett. 102.094503. [12] Y.-S. Choi, K.-W. Seo, S.-J. Lee, Lateral and cross-lateral focusing of spherical particles in a square microchannel, Lab Chip 11 (3) (2011) 460–465, https://doi.org/ 10.1039/C0LC00212G. [13] K. Miura, T. Itano, M. Sugihara, Inertial migration of neutrally buoyant spheres in a pressure-driven flow through square channels, J. Fluid Mech. 749 (2014) 320–330, https://doi.org/10.1017/jfm.2014.232. [14] C. Yuan, Z. Pan, H. Wu, Inertial migration of single particle in a square microchannel over wide ranges of re and particle sizes, Microfluid. Nanofluid. 22 (9) (2018) 102, https://doi.org/10.1007/s10404-018-2120-y.
9
Computational Materials Science 176 (2020) 109542
Y. Huang, et al.
[30] [31] [32] [33] [34] [35] [36] [37]
[38] [39] [40]
[41]
laminar pipe flow, J. Fluid Mech. 621 (2009) 59–67, https://doi.org/10.1017/ S0022112008004977. S. Rubinow, J.B. Keller, The transverse force on a spinning sphere moving in a viscous fluid, J. Fluid Mech. 11 (3) (1961) 447–459, https://doi.org/10.1017/ S0022112061000640. R. Simha, Untersuchungen über die viskosität von suspensionen und lösungen, Kolloid-Zeitschrift 76 (1) (1936) 16–19, https://doi.org/10.1007/BF01432457. D. Joseph, D. Ocando, Slip velocity and lift, J. Fluid Mech. 454 (2002) 263–286, https://doi.org/10.1017/S0022112001007145. P. Saffman, The lift on a small sphere in a slow shear flow, J. Fluid Mech. 22 (2) (1965) 385–400, https://doi.org/10.1017/S0022112065000824. R.C. Jeffrey, J. Pearson, Particle motion in laminar vertical tube flow, J. Fluid Mech. 22 (4) (1965) 721–735, https://doi.org/10.1017/S0022112065001106. E.S. Asmolov, A.L. Dubov, T.V. Nizkaya, J. Harting, O.I. Vinogradova, Inertial focusing of finite-size particles in microchannels, J. Fluid Mech. 840 (2018) 613–630, https://doi.org/10.1017/jfm.2018.95. A. Yazdani, M. Deng, B. Caswell, G.E. Karniadakis, Flow in complex domains simulated by dissipative particle dynamics driven by geometry-specific body-forces, J. Comput. Phys. 305 (2016) 906–920, https://doi.org/10.1016/j.jcp.2015.11.001. B. Zhou, W. Luo, J. Yang, X. Duan, Y. Wen, H. Zhou, R. Chen, B. Shan, Simulation of dispersion and alignment of carbon nanotubes in polymer flow using dissipative particle dynamics, Comput. Mater. Sci. 126 (2017) 35–42, https://doi.org/10. 1016/j.commatsci.2016.09.012. A.L. Blumers, Y.-H. Tang, Z. Li, X. Li, G.E. Karniadakis, Gpu-accelerated red blood cells simulations with transport dissipative particle dynamics, Computer Phys. Commun. 217 (2017) 171–179, https://doi.org/10.1016/j.cpc.2017.03.016. Y. Huang, R.L. Marson, R.G. Larson, Inertial migration of a rigid sphere in plane Poiseuille flow as a test of dissipative particle dynamics simulations, J. Chem. Phys. 149 (16) (2018) 164912, , https://doi.org/10.1063/1.5047923. R.L. Marson, Y. Huang, M. Huang, T. Fu, R.G. Larson, Inertio-capillary crossstreamline drift of droplets in Poiseuille flow using dissipative particle dynamics simulations, Soft Matter 14 (12) (2018) 2267–2280, https://doi.org/10.1039/ c7sm02294h. Y. Huang, R.L. Marson, R.G. Larson, Inertial migration of neutrally buoyant prolate and oblate spheroids in plane poiseuille flow using dissipative particle dynamics simulations, Comput. Mater. Sci. 162 (2019) 178–185, https://doi.org/10.1016/j. commatsci.2019.02.048.
[42] P. Hoogerbrugge, J. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, EPL (Europhysics Letters) 19 (3) (1992) 155, https://doi.org/10.1209/0295-5075/19/3/001. [43] M. Liu, G. Liu, L. Zhou, J. Chang, Dissipative particle dynamics (dpd): an overview and recent developments, Arch. Comput. Methods Eng. 22 (4) (2015) 529–556, https://doi.org/10.1007/s11831-014-9124-x. [44] E. Pep, P. Warren, Statistical mechanics of dissipative particle dynamics, EPL (Europhysics Letters) 30 (4) (1995) 191, https://doi.org/10.1209/0295-5075/30/ 4/001. [45] R.D. Groot, P.B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (11) (1997) 4423–4435, https://doi.org/10.1063/1.474784. [46] J.A. Millan, W. Jiang, M. Laradji, Y. Wang, Pressure driven flow of polymer solutions in nanoscale slit pores, J. Chem. Phys. 126 (12) (2007) 03B617, https://doi. org/10.1063/1.2711435. [47] C.A. Marsh, J.M. Yeomans, Dissipative particle dynamics: the equilibrium for finite time steps, Europhys. Lett. (EPL) 37 (8) (1997) 511–516, https://doi.org/10.1209/ epl/i1997-00183-2. [48] D.C. Visser, H.C.J. Hoefsloot, P.D. Iedema, Modelling multi-viscosity systems with dissipative particle dynamics, J. Comput. Phys. 214 (2) (2006) 491–504, https:// doi.org/10.1016/j.jcp.2005.09.022. [49] J.A. Anderson, C.D. Lorenz, A. Travesset, General purpose molecular dynamics simulations fully implemented on graphics processing units, J. Comput. Phys. 227 (10) (2008) 5342–5359, https://doi.org/10.1016/j.jcp.2008.01.047. [50] J. Glaser, T.D. Nguyen, J.A. Anderson, P. Lui, F. Spiga, J.A. Millan, D.C. Morse, S.C. Glotzer, Strong scaling of general-purpose molecular dynamics simulations on GPUs, Comput. Phys. Commun. 192 (2015) 97–107, https://doi.org/10.1016/j.cpc. 2015.02.028. [51] H. Tollert, Die wirkung der magnus-kraft in laminaren strömungen. i. teil. der quertrieb in sedimentierenden teilchenschüttungen, Chem. Ing. Tech. 26(3) (1954,) 141–150, https://doi.org/10.1002/cite.330260306. [52] Y.W. Kim, J.Y. Yoo, Axisymmetric flow focusing of particles in a single microchannel, Lab Chip 9 (8) (2009) 1043–1045, https://doi.org/10.1039/B815286A. [53] C. Prohm, H. Stark, Feedback control of inertial microfluidics using axial control forces, Lab Chip 14 (12) (2014) 2115–2123, https://doi.org/10.1039/ C4LC00145A.
10