Molecular Dynamics Simulations of Janus Particle Dynamics in Uniform Flow
Aurelien Y. M. Archereau, Shaun C. Hendy, Geoff. R. Willmott
MMolecular Dynamics Simulations of Janus Particle Dynamics in Uniform Flow
Aurelien Y. M. Archereau ∗ Ecole Normale Superieure and Pierre and Marie Curie University, France
Shaun C. Hendy † Te P¯unaha Matatini, Department of Physics,University of Auckland, New Zealand
Geoff R. Willmott ‡ The MacDiarmid Institute for Advanced Materials and Nanotechnology,Departments of Physics and Chemistry,University of Auckland, New Zealand (Dated: March 28, 2018)We use molecular dynamics simulations to study the dynamics of Janus particles, micro- ornanoparticles which are not spherically symmetric, in the uniform flow of a simple liquid. In partic-ular we consider spheres with an asymmetry in the solid-liquid interaction over their surfaces andcalculate the forces and torques experienced by the particles as a function of their orientation withrespect to the flow. We also examine particles that are deformed slightly from a spherical shape. Wecompare the simulation results to the predictions of a previously introduced theoretical approach,which computes the forces and torques on particles with variable slip lengths or aspherical deforma-tions that are much smaller than the particle radius. We find that there is good agreement betweenthe forces and torques computed from our simulations and the theoretical predictions, when the slipcondition is applied to the first layer of liquid molecules adjacent to the surface.
I. INTRODUCTION
Janus particles [1, 2], named after the Roman godwith two faces, are micro- or nanoparticles which are notspherically symmetric. Although many types of asym-metry are possible [3], the typical example is a spheri-cal particle consisting of two hemispheres with distinctsurface properties. Many different types of Janus parti-cles have been reported [1, 2], and possible applicationsoften draw upon the interesting surface wettability, self-assembly, or multifunctional properties of these particles.Multifunctional Janus particles are being used to controlaffinity for human endothelial cells [4] or breast cancercells [5] in biomedical applications. Examples of the im-portance of wettability are provided by the use of am-phiphilic Janus particles for the stabilization of water-in-oil or oil-in-water emulsions [1]. Much as amphiphilicmolecules can form various patterned micelles and vesi-cles in water [6], collections of Janus particles can alsoself-assemble into similar phases [7, 8].Links between surface wettability and micro- tonanoscale fluid dynamics have become of significant in-terest in recent times, particularly in relation to a slipboundary condition [9–13]. At the solid-liquid interface,the usual hydrodynamic assumption is the Dirichlet (non- ∗ [email protected] † [email protected] ‡ [email protected] slip) boundary condition, in which the fluid adjacent tothe surface has zero speed. When this assumption is vi-olated, slip occurs, and non-zero slip has been observedexperimentally for atomically smooth, hydrophobic sur-faces [9, 11]. The non-slip boundary condition is mostcommonly described in the rest frame of the wall usingthe slip length b , defined according to Navier’s approach[10, 14–16] in which the velocity ( v ) of an incompressiblefluid has a component tangential to the surface ( v (cid:107) ) thatis proportional to the shear rate: v (cid:107) = b n · (cid:16) ∇ v + ( ∇ v ) T (cid:17) · ( I − nn ) . (1)Here n is the unit normal to the surface. In a simple shearflow over a planar surface, the slip length represents thedepth at which a linear extrapolation of the tangentialflow velocity at the surface becomes zero.A Janus particle may have an asymmetric hydrody-namic boundary condition, because different areas onthe particle surface can have different slip lengths. Suchan asymmetry can generate interesting dynamics, andthe force and torque on a slip-asymmetric Janus par-ticle in uniform Newtonian flow have been analyticallycalculated at the continuum level [17, 18]. These cal-culations suggest that a Janus particle consisting of twohemispheres with differing boundary conditions will ex-perience a torque, and have a reduced linear drag co-efficient compared to a sphere with a no-slip boundarycondition [17]. The analytical results apply to particles ofradius R in the limit b (cid:28) R , in which case the slip asym-metry is mathematically equivalent to slight deformation a r X i v : . [ phy s i c s . f l u - dyn ] J un of a sphere [19].Slip-induced dynamics of Janus particles are of prac-tical interest for several reasons. Firstly, experimentalevidence for molecular scale slip has been obtained usinga limited range of techniques [9, 11, 20]. The most suc-cessful of these (AFM colloidal probe, and surface forceapparatus) measure the force between two surfaces dueto drainage, rather than directly measuring the shearforce couple between the fluid and the surface. There-fore, studies of slip on micro- and nanospheres could ver-ify and further elucidate the physical nature of molecularscale slip. Secondly, the dynamics provide a mechanismfor manipulation and orientation of particles, without theapplication of external fields, suggesting possible applica-tions as multifunctional particles. Thirdly, the collectivedynamics of Janus particles could be affected. Simula-tions of Janus particle self-assembly [7, 8, 21] have not sofar considered slip.This paper presents molecular dynamics (MD) simula-tions of Janus particles with asymmetric wettability (andconsequently an asymmetry in the slip length) in a sim-ple fluid flow. The values of the force and torque on theparticles are calculated as a function of the slip length ofthe two hemispheres, and the angle between the asym-metry and the flow. Results are compared with the con-tinuum theoretical studies of the slip-asymmetric Janussphere, and of a slightly deformed sphere with homo-geneous wettability. Presently, MD simulations are es-sential for studying slip-asymmetric Janus particles, be-cause the predicted dynamics have not been experimen-tally confirmed, mostly due to competition with Brown-ian motion [17]. Here we make explicit predictions of theforces and torques on particles, whereas previous studieshave been largely descriptive [22, 23]. II. THEORETICAL APPROACHA. The Janus particle
The system studied is a Janus particle modelled as asphere of radius R evenly divided into hemispheres re-ferred to as 1 and 2, with the respective slip lengths b and b . The fluid (viscosity η ) moves at a velocity U relative to the Janus particle. The angle θ between thesphere and the flow is defined in Fig. 1.In general, the effective slip length measured for flowaround a sphere b (cid:107) will not be equal to the effective sliplength measured for flow over a plane b (cid:63) for the sameliquid-solid interaction strength. As shown by Chen etal. [24], when Eq. 1 is used to define slip, for a surfacewith radius of curvature R we have the relation:1 b (cid:107) = 1 b (cid:63) + 1 R (2)where b (cid:107) ≤ b (cid:63) . Note that the sign of the curvature R differs from some previous accounts [14, 17, 18] because θ R U ,η xz FIG. 1. Geometry of the system : cross-section through thecentre of a Janus sphere in the x - z plane. The sphere of radius R is divided into two hemispheres 1 and 2. The unboundeduniform flow U is parallel to the x axis, and the sphere isrotated by an angle θ about the y axis. we define the radius of curvature of the sphere to bepositive (and therefore equal to the sphere radius).We consider the drag force acting parallel to the inci-dent flow, which is in the positive x direction, and thetorque about the y axis. If U has a component in the y direction, there are other force and torque componentsnot considered here. It was previously found [17] that theforce on a Janus particle with asymmetric slip can be ex-panded as a first order approximation to the well-knownStokes drag force: F Stokes = 6 πηRU. (3)For this Janus particle (Fig. 1), the force and torque aregiven by F Janus = 6 πηRU (1 − b (cid:107) + b (cid:107) R ) T Janus = − πηR U ( b (cid:107) − b (cid:107) )2 R cos( θ ) (4)where b (cid:107) (cid:28) R , b (cid:107) (cid:28) R . B. The slightly deformed sphere
We also consider a slightly deformed sphere where thedeformation is aligned with the slip asymmetry. Theforce and torque on a slightly deformed sphere [19] can bedetermined using the spherical polar co-ordinate angles θ (see Fig. 1) and φ to define a point on the surface ofthe sphere which we associate with the radius R d ( θ, φ ).The sphere surface takes the locus of points: R d ( θ, φ ) = R (1 + (cid:15)f ( θ, φ )) (5)where (cid:15) is the amplitude of the deformation on one hemi-sphere ( (cid:15) (cid:28) γ(cid:15) is the amplitude of the deforma-tion on the other hemisphere ( (cid:15)γ (cid:28) f ( θ, φ ) = (1 + γ ) + ( γ − cos ( θ ) is the first order spherical har-monic expansion describing the geometry.Calculation of force and torque on the slightly de-formed sphere with a non-slip boundary condition leadsto the set of equations F Janus = 6 πηRU (1 − (cid:15)/ γ )) T Janus = − πηR U (cid:15) (1 − γ ) cos ( θ ) . (6)The effects of asymmetric slip (Eq. 4) and deformation(Eq. 6) are both first order corrections to the Stokes drag(Eq. 3), so we can add them to describe a deformed par-ticle with slip asymmetry, leading to the set of equations: F Janus = 6 πηRU (1 − (cid:15)/ γ ) − b (cid:107) + b (cid:107) R ) T Janus = − πηR U ( (cid:15) (1 − γ ) + b (cid:107) − b (cid:107) R ) cos ( θ ) . (7)By choosing a deformed sphere with the same slip con-dition b (cid:107) = b (cid:107) = b (cid:107) on both hemispheres, the slip-inducedtorque vanishes, while the first order correction remainsfor the force equation: F Janus = 6 πηRU (1 − b (cid:107) R − (cid:15)/ γ )) T Janus = − πηR U (cid:15) (1 − γ ) cos ( θ ) . (8) III. METHODS AND SIMULATIONSA. Inter-molecular interactions
The molecular dynamics simulations were carried outusing an explicit solvent with inter-molecular interactionsdescribed by a Lennard-Jones pair potential. For exam-ple, the interactions between two fluid monomers ( i , j )follow the equation: U ij ( r ) = 4 E [( σr ) − ( σr ) )] (9)where E represents the depth of the potential well, r isthe distance between the centres of the monomers, and σ is the distance where U ij = 0. The potential is truncatedat 2.5 σ for computational efficiency, at which separation U ij (2 . σ ) = 0 . E .The surface of the particle is divided evenly in twohemispheres referred to as 1 and 2, each interacting withthe fluid monomers by a different Lennard Jones poten-tial, with a different potential well. For example, theinteraction between the monomers of hemisphere 1 andthe fluid monomers follow the equation: U ij ( r ) = 4 E [( σr ) − ( σr ) )] E = A E (10) where A is a parameter allowing adjustments of thestrength of the pair interaction. As we will see below,the fluid-Janus monomer pair interaction is related tothe slip length in that hemisphere, so the slip length b (cid:107) can be controlled by altering A . The Janus monomersof hemispheres 1 and 2 do not interact with each other. B. The Janus Particle
The Janus particle was constructed from 270monomers arranged on the surface of a sphere (Fig. 2,left). The radius of a sphere containing the monomers is4 . σ . The theoretical predictions (above) are expectedto be valid when b (cid:28) R , and our simulations extend tothe upper bound of where we might expect this theoryto apply.The particle behaves as a rigid body : the relativedistance from each monomer on the Janus particle isfixed, and the Janus particle is static in the flow (notranslational nor angular motion). The positions of themonomers on the Janus particle were calculated using aMonte Carlo Scheme so that the monomers were evenlydistributed on the surface. Each monomer ( i , j ) on thesurface is defined by the two angles ( φ i , φ i ). We designedan arbitrary score function which strongly disadvantagesconfigurations in which two monomers are close, definedas S = (cid:88) i (cid:54) = j ( 1 d ij ) , (11)where ( i , j ) refer to the monomers of the Janus parti-cle. The score function can be seen as a repulsive forcebetween close monomers. Convergence was assessed bystudying the first peak of the pair distribution function.A new Monte Carlo design of the Janus particle was donebefore each simulation in order to account for the irreg-ularities of the particle.The deformed homegeneous sphere (Fig. 2, right) wasdesigned using the same Monte Carlo scheme, with a uni-form pair interaction potential ( A = A ). C. Molecular dynamics simulations
The equations of motion are integrated using a Verletmethod with the timestep set at ∆ t = 0 . τ LJ , where τ LJ is the Lennard-Jones timescale ( m σ/E ) / and m is the monomer mass. All simulations are run in Lennard-Jones units, where σ is the unit of distance, E is theunit of energy, m the unit of mass and τ LJ the unitof time, without loss of generality. The simulations werecarried out with LAMMPS [25] in the canonical ensemble(constant NVT). The temperature was set at 1 . (cid:15)/k B (where k B is the Boltzmann constant) using a Langevinthermostat with a damping factor (coupling parameter) FIG. 2. Snapshot of the Janus particle (left) and the deformedsphere (right). Both particles are made of 270 Lennard Jonesmonomers. Note that the deformation of the sphere has beenexaggerated here for the sake of visibility. set at 20 τ LJ . Initial velocities for the monomers weredrawn from a Maxwell-Boltzmann distribution.For the study of the Janus particle, the system con-sisted of 46288 fluid monomers of mass m in a cubic boxwith sides of length 39.62 σ ( ρ fluid = 0 . σ − ). Periodicboundary conditions are used in the three dimensions.The flow is created by moving the Janus particle at aconstant speed ( U = 0 . σ/τ ) in a static fluid. The ef-fect of the finite box size was assessed by monitoring theaverage fluid speed at the boundary of the box, whichwas found to be less than 3% of the speed of the Janusparticle in all simulations reported here.We computed the fluid’s shear viscosity η using theGreen-Kubo relation, which relates the viscosity to theautocorrelation function of a diagonal component of thepressure tensor. This was calculated using a system withonly fluid monomers in a box with periodic boundaryconditions and with the fluid at the same fluid densityas the fluid in the Janus particle simulations. The cal-culation converged in 2 × MD steps to a value of η = 3 . ± . m /στ LJ .Slip lengths b (cid:63)P for flow over a planar surface were com-puted as a function of the fluid-Janus particle interactionenergy E = AE by simulating a Poiseuille flow betweentwo solid walls in the yz plane. The walls were composedof (111) planes of face-centered cubic lattice (i.e. a close-packed surface) with a surface density corresponding tothat of the Janus particle. The Poiseuille flow was cre-ated by adding a constant force F x = 0 . E /σ to allthe fluid monomers. The velocity profile across the chan-nel was then fitted with a parabolic curve to determinethe slip length at the walls for each value of the fluid-wallmonomer interaction energy E = AE . The exact posi-tion of the walls in this calculation is discussed furtherbelow.The number of monomers in our structures varied be-tween around 10 000 and 60 000. Simulations were al-lowed 2 × time steps for equilibration, followed by8 × time steps for production runs. For each type ofsimulation 10 simulations were run to estimate the meanand standard deviation of computed quantities. IV. RESULTSA. Effect of thermostat on the slip length
Although, as discussed above, we use a Langevin ther-mostat for most of the work reported here, it is impor-tant to check whether this has an impact on computedslip lengths. Figure 3 shows the slip length for a pla-nar surface b (cid:63)P as a function of E = AE for both theLangevin thermostat (Fig. 3, blue line), a Nos´e-Hooverthermostat (Fig. 3, green line) and several values fromRef [22]. The Nos´e-Hoover thermostat used a couplingparameter of 1.0 τ LJ , which is the same value used inRef [22]. The differences are relatively small, but ourchoice of thermostat enabled us to avoid the ‘flying icecube’ effect [26] over a wider range of simulation param-eters. This choice of thermostat leads to slightly largervalues of slip length for small E than the simulations inRef [22]. s li p l eng t h pair interation strength E Ref [22]Nosé−HooverLangevin FIG. 3. Slip length in σ units, as a function of E = AE onthe surface. Literature values of b (cid:63)P from Ref [22] are shownin black, we were able to replicate these results (green). Theresult with the Langevin thermostat is shown in blue. B. Slip for a homogeneous sphere
We can compute the slip length for a homogeneoussphere in order to check the applicability of Eq. 2. Us-ing Eq. 4 with b (cid:107) = b (cid:107) , an effective slip length b (cid:107) S canbe extracted from the measured force on a simulated ho-mogeneous sphere for various values of the interactionstrength E : b (cid:107) S = R (1 − F πηRU ) (12)The calculations to find the Stokes-extracted slip length b (cid:107) S and the corrected slip length b (cid:107) C both use the radiusof the sphere R . There are two different, but reasonabledefinitions for the radius of the sphere (see Fig. 4) : First-layerradius First-layerslip lengthWallradiusradius Wall slip lengthInter-atomicdistance Inter-atomicdistance
FIG. 4. Both the effective radius R of the simulated sphere(left) and the slip length (right) can be defined using eitherthe position of the solid wall (surface), or the first layer offluid monomers. The latter definition was consistently usedhere.FIG. 5. The corrected values of the slip lengths b (cid:107) C (purple),along with the slip length values b (cid:107) S computed from the sim-ulation of the homogeneous sphere (red). the distance from the center of the sphere to the surfacemonomers ( R = 4 . σ ), or the distance from the centerof the sphere to the first layer of fluid ( R = 5 . σ ). Thedifference is 20%, and therefore not negligible.Analysis of the simulations shows that the latter defini-tion has the better agreement, and is used in what followsthroughout. In particular, we find that the slip lengthcalculated from the homogeneous particle b (cid:107) S is in goodagreement with the values calculated from the Poiseuilleflow if the effective radius is chosen to be R = 5 . σ (seeFig. 5). This makes physical sense, as it is in the firstlayer of fluid where the boundary condition (Eq. 1) is de-fined. At low slip length the agreement is especially goodwhereas at high slip length a difference is visible, whichis consistent with the fact that the theory used the the-oretical assumption b (cid:28) R . This is not the case here, asthe maximum slip length is nearly half the radius. In thefollowing, the slip length used to calculate the expectedforce and torque is the slip length b (cid:107) S extracted from thehomogeneous sphere simulations. C. Force and torque for a Janus particle
To assess the agreement of the simulations with Eqs. 4,we calculated the force and torque on simulated Janusparticles with two different values of the interactionstrength R as a function of the angle θ , and the twoslip lengths b (cid:107) and b (cid:107) imputed from the calculations inthe previous section. The results, shown in Fig. 6, are ingood agreement with the the expected values from thetheory. The functional agreement is high (the fits havePearson coefficients larger than 0.994) and the magni-tude of the simulated values are always within 4% of theexpected ones.One can detect a small variation of frequency π in thesubplot (a) , representing a variation in the force as afunction of the angle. This variation is likely to be theconsequence of a second order term, as the values of sliplength are of the same order of magnitude as the radiusof the sphere. One other possibility is the design of thesphere, which is not regular along the equator betweenthe two poles (see Fig. 2). The amplitude of this effect is,however, small : 0.27% of the force value. We concludethat the theory does an excellent job of describing theforces and torques on our simulated Janus particle. D. Force and torque for a deformed sphere
Finally we calculated the force and torque using simu-lations for a deformed sphere, which should be describedby Eqs. 13, as a function of the angle θ , and the defor-mation parameters (cid:15) and γ (see Fig. 7). We designed adeformed sphere with homogeneous interaction strength E where the slip length b (cid:107) is as small as possible ( A = 1 . b (cid:107) = 1 . σ ). Again, wefind that the simulations are in good agreement with theexpected values. The differences are larger than in thecase of the Janus particle simulations, quite possibly dueto the combination of the slip and the deformation.In particular, the simulation of the force as a functionof the angle (subplot (a) ) exhibits a sinusoidal deviationof frequency π from the expected behavior. Again, thismight be a second order effect as in the Janus particle setof simulations, or a combined effect of the slip and the de-formation. The amplitude of this effect is however small: less than 2% of the value of the force. A theoreticalanalysis to include second order terms is not trivial.However, when fitting the data with a free multiplica-tive parameter B as follows T = − B ∗ πηR U (cid:15) (1 − γ ) cos ( θ ) (13)we find that B = 0 .
72, which is of the order of (1 − b (cid:107) R ) =0 .
69 in these simulations. This correction can be physi-cally interpreted as the first order torque acting from anadjusted distance R adjusted = R (1 − b (cid:107) R ).
62 63 64 65 (a) F o r c e ( un i t s o f f o r c e )
58 59 60 61 62 63 64 65 (b) 60 61 62 63 64 65 66 67 68 (c)−20−15−10−5 0 5 10 15 20 0 0.5 1 1.5 2(d) T o r que ( un i t s o f f o r c e ) θ (radians) −25−20−15−10−5 0 5 10 15 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3(e)b (units of sigma) 0 5 10 15 20 25 30 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3(f)b (units of sigma) FIG. 6. Force and torque values as a function of θ , b (cid:107) and b (cid:107) for the Janus particle simulation. The values calculated fromsimulations are plotted in purple, and the expected values (Eq. 4) are shown in green. Fixed parameters : for a,d : A = 1 . A = 0 .
7. For b,e : A = 0 . θ = 0. For c,f : A = 1 . θ = 0.
63 64 65 66 (a) F o r c e ( un i t s o f f o r c e )
63 64 65 66 67 68 (b) 62 63 64 65 (c)−15−10−5 0 5 10 15 0 0.5 1 1.5 2(d) T o r que ( un i t s o f f o r c e ) θ (radians) −15−10−5 0 5 0 0.02 0.04 0.06 0.08 0.1(e) ε (units of sigma) −25−20−15−10−5 0 5 0 0.2 0.4 0.6 0.8 1(f) γ (units of sigma) FIG. 7. Force and torque values as a function of θ , (cid:15) and γ for the deformed sphere simulation. The values are plotted inpurple, and the expected values (Eqs. 13) are shown in green. On subplot (a) a fit with a cosine is presented and shows asecond-order effect (red). On subplots (d,e,f ) a fit is shown in red. Fixed parameters : for a,d : (cid:15) = 0 . γ = 0 .
3. For b,e : θ = 0, γ = 0 .
3. For c,f : θ = 0, (cid:15) = 0 . V. CONCLUSION
For small Janus particles in a simple liquid flow, wefind good agreement between the forces and torquescomputed from molecular dynamics simulations and acontinuum theoretical description based on inhomoge-neous perturbations to a Stokes flow about a homogenous sphere. We have considered particles with an asymmetryin the solid-liquid interaction strength over their surfacesand particles that are deformed slightly from a spheri-cal shape but have a uniform interaction strength. Thetheory gives excellent quantitative agreement with thesimulations, particularly in the limit of small perturba-tions or slip lengths as used in the theory, and describesthe functional dependence of forces and torques on anglewell. The results demonstrate the importance of apply-ing the slip condition at the first layer of fluid molecules,and of adjusting the slip length to account for surfacecurvature. This work suggests that molecular dynamicssimulations could be used to study related problems withmore complicated geometries and flows.
ACKNOWLEDGMENTS
The molecular dynamics simulations were run usingthe LAMMPS software (Large-Scale Molecular Massively Parallel Simulator [25]). The calculations were run onNeSI Pan Cluster, part of the Centre for eResearchhosted by the University of Auckland. [1] A. Perro, S. Reculusa, S. Ravaine, E. Bourgeat-Lami,and E. Duguet, J. Mater. Chem. , 3745 (2005).[2] A. Walther and A. H. E. Muller, Soft Matter , 663(2008).[3] S. C. Glotzer and M. J. Solomon, Nature Materials ,557 (2007).[4] M. Yoshida, K.-H. Roh, S. Mandal, S. Bhaskar, D. Lim,H. Nandivada, X. Deng, and J. Lahann, Advanced ma-terials , 4920 (2009).[5] L. Y. Wu, B. M. Ross, S. Hong, and L. P. Lee, Small ,503 (2010).[6] J. N. Israelachvili, D. J. Mitchell, and B. W. Ninham,Journal of the Chemical Society, Faraday Transactions 2:Molecular and Chemical Physics , 1525 (1976).[7] G. Rosenthal, K. E. Gubbins, and S. H. L. Klapp, TheJournal of Chemical Physics , 174901 (2012).[8] Z. W. Li, Z. Y. Lu, Z.-Y. Sun, and L. J. An, Soft Matter , 6693 (2012).[9] C. Neto, D. R. Evans, E. Bonaccurso, H.-J. Butt, andV. S. J. Craig, Rep. Prog. Phys. , 2859 (2005).[10] E. Lauga, M. P. Brenner, and H. A. Stone, “Handbookof experimental fluid dynamics,” (Springer, New York,2005) Chap. 15: Microfluidics: The No-Slip BoundaryCondition.[11] L. Bocquet and E. Charlaix, Chem. Soc. Rev. , 1073(2010).[12] L. Bocquet and J.-L. Barrat, Soft Matter , 685 (2007). [13] S. C. Hendy and N. J. Lund, Phys. Rev. E , 066313(2007).[14] D. Einzel, P. Panzer, and M. Liu, Phys. Rev. Lett. ,2269 (1990).[15] C. L. M. H. Navier, Mem. Acad. Sci. Inst. Fr. , 839(1827).[16] N. J. Lund, X. P. Zhang, K. Mahelona, and S. C. Hendy,Phys. Rev. E , 046303 (2012).[17] G. R. Willmott, Physical Review E , 055302 (2008).[18] G. R. Willmott, Physical Review E , 066309 (2009).[19] J. Happel and H. Brenner, Low Reynolds Number Hy-drodynamics (Noordhoff International Publishing, TheNetherlands, 1973).[20] C. I. Bouzigues, L. Bocquet, E. Charlaix, C. Cottin-Bizonne, B. Cross, L. Joly, A. Steinberger, C. Ybert, andP. Tabeling, Phil. Trans. R. Soc. A , 1455 (2008).[21] M. M. Moghani and B. Khomami, Soft Matter, 2013, 9, , 4815 (2013).[22] A. Kharazmi and N. V. Priezjev, Journal Of ChemicalPhysics , 234503 (2015).[23] E. Bianchi, A. Z. Panagiotopoulos, and A. Nikoubash-man, Soft matter , 3767 (2015).[24] W. Chen, R. Zhang, and J. Koplik, Phys. Rev. E ,023005 (2014).[25] S. Plimpton, Journal of computational physics , 1(1995).[26] S. C. Harvey, R. K.-Z. Tan, and T. E. Cheatham, Journalof Computational Chemistry19