Pressure Gradients Fail to Predict Diffusio-Osmosis
aa r X i v : . [ c ond - m a t . s o f t ] J a n Pressure Gradients Fail to Predict Diffusio-Osmosis
Yawei Liu
Beijing Advanced Innovation Center for Soft Matter Science and Engineering,Beijing University of Chemical Technology, Beijing 100029, China
Raman Ganti and Daan Frenkel ∗ Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, UK (Dated: July 25, 2018)We present numerical simulations of diffusio-osmotic flow, i.e. the fluid flow generated by aconcentration gradient along a solid-fluid interface. In our study, we compare a number of distinctapproaches that have been proposed for computing such flows and compare them with a referencecalculation based on direct, non-equilibrium Molecular Dynamics simulations. As alternatives, weconsider schemes that compute diffusio-osmotic flow from the gradient of the chemical potentialsof the constituent species and from the gradient of the component of the stress tensor parallel tothe interface. We find that the approach based on treating chemical potential gradients as externalforces acting on various species agrees with the direct simulations, thereby supporting the approachof Marbach et al. (J Chem Phys 146, 194701 (2017)). In contrast, an approach based on computingthe gradients of the microscopic pressure tensor does not reproduce the direct non-equilibriumresults.
INTRODUCTION
Flow in macroscopic channels is driven by pressuregradients or body forces, such as gravity. In contrast,flow in nano-structured materials (e.g. nano-channels),is usually dominated by phoretic effects, where transportis caused by thermodynamic gradients acting near inter-faces. The key point to note is that the gradients respon-sible for phoretic flow (e.g. electrical fields, concentrationgradients or thermal gradients), cannot cause bulk flow:they only act in narrow interfacial layers where the fluidexperiences surface-specific interactions. In view of theincreasing importance of micro- and nano-fluidic devicesand self-propelling particles (see, e.g. 1–7) it becomes im-portant to be able to predict the strength of phoretic flowphenomena based on knowledge of the microscopic inter-actions between the atoms or molecules in the system. Inthe present paper, we focus on the numerical predictionof diffusio-osmosis, where a chemical potential gradientalong a solid surface drives the flow [8–10]. Tradition-ally, such flows have been described using a continuumpicture, where the material properties were characterisedby (local) thermodynamic quantities and the flow wascomputed, assuming that the (Navier-)Stokes equationholds near the surface. Such a macroscopic perspectiveis appropriate in systems where the fluid-interface in-teraction acts over a range that is much larger than atypical molecular size, as is the case for electro-osmosisof dilute electrolytes. However, as most inter-molecularforces have a range comparable to a molecular diameter,the continuum picture is not expected to hold for mostdiffusio-osmotic phenomena.In what follows, we focus on a situation where thesurface-fluid interaction is short-ranged. Specifically, weconsider MD simulations of diffusio-osmosis of a neutral solvent containing a neutral solute, near a crystallinesolid surface that is, on average, flat.Thus, to determine the flow velocity, it is critical toaccurately calculate the surface force induced by the con-centration gradient. The flow can then be obtained fromMD simulations by applying the force to fluid particles.In order to validate our method, we developed an in-genious non-equilibrium simulation technique to directlycompute the flow.
THERMODYNAMIC GRADIENTS
We first consider an atomically flat wall at z ( x, y ) ≤ z >
0. Thefluid contains a majority component ( A : ‘solvent’) and aminority component ( B : ‘solute’). The fluid, as a whole,is maintained at constant bulk pressure and constanttemperature. When a concentration gradient in B (and,via the Gibbs-Duhem relation, also in A ) is imposed inthe fluid along x , flow occurs due to a pressure gradi-ent at the interface. As the concentration gradients in A and B are not independent, we will focus on the phoreticeffect of the concentration gradient of B . It should beborne in mind that this effect, includes the effect of theconcentration gradient in A . The most intuitive (but, aswe will show, incorrect) method to obtain the force act-ing on fluid particles near an interface is to calculate theforce due to the pressure gradient acting on a small vol-ume element, and then obtain the force per particle bydividing the force per volume by the local number den-sity. The transverse component of the pressure tensor at z ( p xx ( z )) depends on x only through its dependence onthe spatial variation of the bulk concentration ( ρ B ) or,equivalently, of the chemical potential ( µ B ) of the speciessubject to a concentration gradient: f V ( z ) = − ∂p xx ( z ) ∂ρ B ∂ρ B ∂x . (1)In the case of a sufficiently small concentration gradi-ent, we can assume local thermodynamic equilibrium(LTE) – deviations from LTE are expected to be of higherthan linear order in the concentration gradient. Thus,rather than computing the local pressure tensor in a non-equilibrium system, we can compute it in equilibrium as afunction of concentration. The pressure gradient is thencomputed using f V ( z ) = − ∂p xx ( z ) ∂ρ B ∂ρ B ∂x ≈ − p xxρ B +∆ ρ B ( z ) − p xxρ B − ∆ ρ B ( z )2∆ ρ B · ∇ ρ B , (2)where ∆ ρ B is a small change in the bulk concentrationof species B . Note that we will always assume that thebulk pressure is constant.The local pressure tensor at position r is defined as theensemble average of the negative of the stress tensor at r [11, 12]. It contains two terms: a kinetic term, arisingfrom the change in momentum due to particles crossingthe boundaries of an elemental volume at r , and a con-figurational term, related to the change in momentumdue to intermolecular interactions between the particles.There is, however, a problem: for inhomogeneous systems(e.g. fluid near a solid wall), the configurational compo-nent of the pressure tensor cannot be defined uniquely.For the computation of surface tension, this ambiguityhas no effect [12], but for the computation of phoreticflow, the problem does not go away.In the present work, we explore the predicted phoreticflow for two different definitions of the local pressuretensor [13, 14] and use these definitions to compute thetransverse pressure as a function of z in a series of slabsparallel to the wall. We employed both the Irving-Kirkwood definition, in which an intermolecular forcecontributes to the local pressure in every slab betweentwo molecules [11], and the Virial definition in which anintermolecular force contributes to the local pressure inthe slab(s) where the two molecules are located [15, 16].Alternatively, we can use local thermodynamics tocompute the force driving flow. Consider an n -component fluid mixture at constant temperature, T .The Gibbs-Duhem relation can be written as V dp = P ni =1 N i dµ i where N i is the number of particles of species i in volume V , p the pressure and µ i the chemical poten-tial of species i . Let us denote the number density ofspecies i in the mixture by ρ i . Then, dp = P ni =1 ρ i dµ i .A concentration gradient of component i along x will leadto a chemical potential gradient ∂µ i /∂x . As the pressureremains constant in the bulk, Gibbs-Duhem relation re-duces to 0 = P ni =1 ρ bulk i ( x ) ( ∂µ i /∂x ). At a position z FIG. 1. (a) Simulation box used to compute the force andflow profiles via the forces obtained from Eq. 2 and Eq. 4.(b) Simulation box used in the direct non-equilibrium MDsimulations with explicitly imposed concentration gradients.The blue particles represent the solvent ( A ), the green par-ticles represent the solute ( B ), the red and yellow particlesrepresent the solid particles in the top wall, and the blackand silver particles represent the solid particles in the bottomwall. near the interface, a pressure gradient remains, giving aforce per unit volume f V ( z ) = (cid:18) − ∂p ( z, x ) ∂x (cid:19) = n X i =1 (cid:0) ρ i ( z, x ) − ρ bulk i ( x ) (cid:1) (cid:18) − ∂µ i ∂x (cid:19) . (3)We can interpret ( − ∂µ i /∂x ) as the force per-particle act-ing on the particles of species i . This expression is conve-nient, because the imposed chemical potential gradientsare constant throughout the system. In the bulk, thecomposition is such that the forces balance (the bulkpressure equilibrates rapidly). Upon approaching thewall, the concentration of different components changes,leading to non-zero net forces. In other words, particlesof a given species experience the same force regardlessof their distance from the interface. The force acting onspecies i is then f i = (cid:18) − ∂µ bulk i ∂x (cid:19) = (cid:18) − ∂µ bulk i ∂ρ i (cid:19) P · ∇ ρ i . (4)We now have two approaches (Eq. 2 and 3) for computingthe force driving diffusio-osmotic flow. Eq. 2 is a mechan-ical expression while Eq. 3 is thermodynamic. Of course,in steady state, these forces are exactly balanced by theforce due to the gradient in the shear stress in the mov-ing fluid. To test which, if either, of these microscopicexpressions is correct, we performed MD simulations onthe simple model system mentioned above. MOLECULAR DYNAMICS SIMULATIONS
We performed MD simulations of a fluid mixture com-posed of solvent ( A ) + solute ( B ) particles, confined be-tween two parallel solid walls (Fig. 1). All particles have FIG. 2. (a) The solute mole fraction profile along z from an equilibrium simulation at ρ B = 0 .
04. (b) The bulk concentrationprofiles in the non-equilibrium simulations. (c) The average density profiles in the diffusio-osmosis region in the direct non-equilibrium simulation at ∇ ρ B = 0 . ρ B = 0 . the same molecular diameter σ . In what follows, we use σ as our unit of length. Interactions between these parti-cles are given by a Lennard-Jones potential truncatedand shifted at 4 σ , U αβ ( r ) = 4 ǫ αβ (cid:2) ( σ/r ) − ( σ/r ) (cid:3) ( α, β ∈ { A, B, top , bottom } ), with ǫ interaction energy.To narrow our exploration, we focused on the ideal solu-tion composed of identical solvent and solute particles inthe bulk, but different values of ǫ with the bottom wall.We chose ǫ AA = ǫ BB = ǫ AB = ǫ A,top = ǫ B,top ≡ . ǫ and ǫ B,bottom = 2 ǫ A,bottom = 1 . ǫ . Thus, the pressure differ-ence only appears in the fluid near the bottom wall whena concentration gradient is imposed. In the remainder ofthis paper, we use ǫ as our unit of energy and the mass m of the fluid particles as our unit of mass.All simulations were carried out using the LAMMPSpackage [17] in an isothermal, isobaric ( N P zz T ) ensem-ble. The top wall was used as a barostat and was other-wise free to move in the horizontal and vertical direction.A constant external force along z was exerted on the topwall to maintain the bulk pressure at a value p ex = 0 . x and y directions. The solid walls were composed of solid par-ticles placed on a FCC lattice with lattice spacing of 1 . ǫ/σ equilibrium length 1 . / √ σ . During the simulations,the layer of solid particles in the bottom wall furthestremoved from the interface was rigidly anchored. Thevelocity-Verlet algorithm with a time step of 0 .
001 wasused to integrate the equations of motion, and a Nos´e-Hoover thermostat with a time constant of 0 . T = 0 . × − × steps to obtainsufficient statistics.The calculation of ∂p xx ( z ) /∂ρ B required several equi-librium simulations at different uniform bulk concentra-tions. These simulations could be carried out in a rel-atively small simulation box as shown in Figure 1(a).The box dimensions were L x = 16 . L y = 9 .
86, and h L z i = 29 . L z fluctuates, and depends on the soluteconcentration). The system contained 2640 fluid parti- cles. To compute the composition-dependence of p xx ( z ),we performed simulations where we varied the concen-tration of the solute B , while keeping the total num-ber of particles fixed. From the numerical estimate of ∂p xx ( z ) /∂ρ B , we computed the corresponding force us-ing Eq. 2 with ρ B = 0 .
02, ∆ ρ B = 0 .
01, and various ∇ ρ B . We verified that our estimate of the pressure gra-dient did not depend on the size of ∆ ρ B . Subsequently,we converted the force per unit volume to a force perparticle, by dividing by the total number density profile ρ ( z ). These per-particle forces were then applied in anon-equilibrium simulation of the fluid with solute den-sity ρ B to determine the flow profile for a given ∇ ρ B [Fig. 1(a)]. Similarly, starting from Eq. 4, we can com-pute the forces that would result from the gradient ofthe chemical potentials. These per-particle forces wereapplied to the solute and solvent particles.To validate our approach, we performed direct non-equilibrium simulations where a concentration gradientwas explicitly imposed. Figure 1(b) shows the simula-tion box in which the fluid mixture has a constant bulkconcentration gradient along x . Non-equilibrium simu-lations were carried out to measure the flow profile atdifferent values of ∇ ρ B . In this case, we employed boxesthat contain several regions: two source regions withthe width of 4 σ , a diffusio-osmosis region with variouswidths, and a transition region of width ∼ σ betweenthe two source regions. During the simulations, every500 steps, the identities of the fluid particles in thesesource regions were reset to maintain constant concen-trations and a steady gradient along x . In the low con-centration source region, ρ B = 0 so that all fluid particleswere reset to the solvent type. In the high concentrationsource region, ρ B = 0 .
04. As the concentration variesonly close to the wall, the fraction of selected particlesreset to the solute type in each slab parallel to the wallis equal to the local solute mole fraction calculated viaan equilibrium simulation at ρ B = 0 .
04 [Fig.2(a)]. In thetransition region, we set ǫ A,bottom = ǫ B,bottom = 0 .
55 toprevent diffusio-osmosis in this region. A steady flow canbe achieved with an ingenious design: During all simula-tions, the bottom wall is fixed by freezing particles in thelast layer of the bottom wall. As the top wall is free, itmoves at the same velocity as the bulk fluid. We trackedthe position of the top wall (i.e. the position of yellowparticles in the top wall) and redefined the position ofsource regions each time the identities of fluid particlesare reset. The concentration gradient along x dependson the box size. In this work, the box size is L x = 36 . .
6, and 92 . L y = 9 .
86, and h L z i = 29 . x for the non-equilibrium simulations. The figureshows that the concentration profile is linear between thelimiting values imposed in the source regions. The con-centration gradients in the three independent simulationswere ∇ ρ B = 0 . . . z from the simulation of ∇ ρ B = 0 . .
02, the density profiles fromthe equilibrium simulation at ρ B = 0 .
02 were also plot-ted for comparison. The results show good agreementfor the local densities from the two simulations, indicat-ing that in the non-equilibrium simulation, fluid statesare still close to equilibrium, validating the assumptionsunderlying our calculation of the surface force via Eqs. 1and 3.
RESULTS AND DISCUSSION
In order to calculate the surface force at ρ B = 0 . ρ B = 0 .
01 and ρ B = 0 .
03. Figures 3(a) and (b) show thepressure profiles along z near the wall at ρ B = 0 .
01 us-ing the Irving-Kirkwood and virial definitions. In thebulk, where fluid is homogeneous and far away fromthe wall, all definitions lead to the same value since p zz = p xx = p ex = 0 . p zz from the Irving-Kirkwood definition is, as expected,equivalent to the bulk pressure, reflecting mechanicalequilibrium along z [Fig. 3 (a)], while the virial expressionfor p zz is not constant [Fig. 3 (b)]. For p xx , the differ-ent expressions for pressure result in different, oscillatingprofiles near the wall.The chemical potential for component i is given by µ i = µ i + k B T ln ρ bulki + µ exci , with k B the Boltzmann con-stant. µ i denotes a (constant) reference value and µ exci denotes the excess chemical potential due to intermolec-ular interactions. Because the bulk solutions are ideal, µ exci does not depend on the concentration of B . Thus, at ρ B = 0 .
02, with ρ bulkA = 0 .
74 and ρ bulkB = 0 .
02 [Fig. 2(b)],if ∇ ρ B = 1 .
0, we obtain f A = 1 .
14 and f B = − . ρ B = 0 .
02 and ∇ ρ B = 1 . f ave ( z ) = [ ρ A ( z ) f A + ρ B ( z ) f B ] /ρ ( z ). As shown in thefigure, the two expressions for the surface force (Eq. 1and 3) produce significantly different results near the in-terface. Unsurprisingly, the forces calculated from thechemical potential gradients (Eq. 3) are concentratedwhere there is an excess of solute [Fig. 2(b)]. However,the forces calculated via the local pressure tensors i.e.the Irving-Kirkwood and virial definitions (Eq. 1), ex-tend over larger distances and fluctuate strongly.The fact that pressure tensor and chemical poten-tial routes yield different force profiles implies that theywould predict different flow profiles. At most, one ofthese can be correct. To test whether the computedflow profiles are correct, we applied the force profilesthat we computed to the fluid mixture at ρ B = 0 . z forfixed ∇ ρ B . Figure 4(a) shows the predicted velocity pro-files at ∇ ρ B = 0 . ∇ ρ B wasalso computed directly in a non-equilibrium simulation[Fig. 1(b)]. The result is shown in Fig. 4(b). We see thatthe velocity profile that follows from the direct simulationdiffers markedly from the one obtained from the pressuregradients. However, it agrees quite well with the pre-dictions based on the chemical-potential gradients. Thelatter agreement was also observed for two other concen-tration gradients ( ∇ ρ B = 0 . . v s ) obtained fromdifferent methods. The results are shown in Fig. 4(d).We see that prediction of v s based on the chemical-potential gradients is in excellent agreement with thoseobtained from direct simulations. However, v s predictedbased on the pressure gradients is significantly different.The failure of mechanical expressions near the inter-face is consistent with our recent calculations on thesolutal Marangoni effect [18] and thermo-osmosis [19].Therefore, it is not entirely surprising that the methodalso fails here. Yet, where the present result differsfrom our earlier studies is that the pressure gradientmethod also predicts an incorrect value of the bulk ve-locity [Fig. 4(d)]. More interestingly, the flow velocityprofile computed with the (presumably correct) chemical-potential-gradient method shows an overshoot near thewall (the effect is clearest for larger solute concentra-tions). This observation is interesting because if the flowprofile could be computed from the force profile usingthe Stokes equation (i.e. assuming a position indepen-dent viscosity) then an overshoot is not possible if theforce always has the same sign (as it does – see Fig. 3(c), FIG. 3. (a) The Irving-Kirkwood pressure profiles at ρ B = 0 .
01. (b) The virial pressure profiles at ρ B = 0 .
01. (c) The averageper-particle force profiles at ρ B = 0 .
02 and ∇ ρ B = 1 . ρ B = 0 .
02 and ∇ ρ B = 0 . green dotted curve).To understand the failure of the pressure-gradient ap-proach, it is useful to revisit the thermodynamic descrip-tion of diffusio-osmotic transport. On a macroscopiclevel, it is the gradient in the surface free-energy density γ (which for fluid-fluid interfaces is equal to the surfacetension) that determines the flow: in the case of fluidinterfaces, this is the well-known Marangoni effect [18].If we consider the variation of γ with the chemical po-tential of the species in a binary mixture, we can write (cid:18) ∂γ∂x (cid:19) P,T = (cid:18) ∂γ∂µ B (cid:19) (cid:18) ∂µ B ∂x (cid:19) + (cid:18) ∂γ∂µ A (cid:19) (cid:18) ∂µ A ∂x (cid:19) = − (cid:20) Γ B (cid:18) ∂µ B ∂x (cid:19) + Γ A (cid:18) ∂µ A ∂x (cid:19)(cid:21) (5)where Γ i ( i = A or B ) is the Gibbs adsorption of species i at the interface. It is clear that this thermodynamic ex-pression immediately yields a relation between the driv-ing force of the Marangoni flow and the chemical poten- tial gradients, which is in agreement with our observa-tion that the microscopic expression [Eq. (3)] couples thechemical potential gradient to the local excess density ofthe corresponding species.For a liquid-liquid (or a liquid-flat wall) interface, wecan use the Kirkwood and Buff expression to relate thesurface tension to the pressure tensor γ = Z ∞−∞ p − p xx ( z ) d z (6)where p is the hydrostatic pressure. DifferentiatingEq. (6) at constant pressure and temperature gives (cid:18) ∂γ∂x (cid:19) P,T = − Z ∞−∞ (cid:18) ∂p xx ( z ) ∂x (cid:19) d z. (7)And hence, the driving force for flow is, in that case, re-lated to the gradient of the pressure tensor. Yet, whenthe solid surface is not flat, but atomistically structured,Eq. (6) no longer holds. The latter integral would yieldthe surface stress . The gradient of the surface stress isnot the driving force for particle transport, and hencefor structured walls, we should not expect to obtain thedriving force for diffusio-osmotic flows from stress gradi-ents.The previous arguments can be numerically testedby multiplying the per-particle force profiles shown inFig. 3(c) by ρ ( z ) and integrating from the surface into thebulk. The surface tension gradient predicted by Eq. (7)using the virial and Irving-Kirkwood pressure expressionsis 1.07 as opposed to -5.81 predicted by Eq. (5). Thesenumerical results confirm that for a structured surface,the stress gradient is not the driving force.For an ideally flat wall, however, we would expectEq. (7) to predict accurately the driving force. We inves-tigate the latter case by applying our methods to fluidinteracting with a surface via specular boundary con-ditions. Fluid-fluid interactions remain unchanged, butfluid-wall interactions are governed by U fw ( z ) = 4 ǫ fw (cid:2) ( σ/z ) − ( σ/z ) (cid:3) (8)where ǫ B,w = 2 ǫ A,w = 1 . ǫ . In this system, the wallis simply a flat surface located at z = 0. Therefore, FIG. 5. (a) Density profiles for species A and B at ρ B = 0 . ρ B = 0 .
023 and ρ B = 0 . ∇ ρ B = 1 . fluid atoms at a position z near the surface will experi-ence the same wall force for all x and y . Changing thewall potential required increasing the external pressureto p ex = 0 .
122 to maintain a stable system.In order to achieve sufficient signal via the pressuregradient approach, ∆ ρ B was increased to 0 .
03 in Eq. (2).Fig 5(a) and (b) show density and pressure profiles forfluid interacting with the specular wall. Comparing withFig 3(a) and (b), it is clear that removing transverseforce contributions from the wall significantly changesthe behavior of fluid near the interface.The surface tension gradient can be computed by in-tegrating the per-volume force profiles shown in Fig 5(c)from z = 0 to z = 5, beyond which none of the methodspredict any effect. The results were − . ± .
01 usingEq. (5), − . ± . − . ± . ∗ Corresponding author: [email protected][1] H. Stone, A. Stroock, and A. Ajdari, Annu. Rev. FluidMech. , 381 (2004).[2] T. M. Squires and S. R. Quake, Rev. Mod. Phys. , 977(2005).[3] L. Bocquet and E. Charlaix,Chem. Soc. Rev. , 1073 (2010).[4] L. Bocquet and P. Tabeling, Lab Chip , 3143 (2014).[5] R. Golestanian, T. B. Liverpool, and A. Ajdari, Phys.Rev. Lett. , 220801 (2005).[6] W. F. Paxton, S. Sundararajan, T. E. Mallouk, andA. Sen, Angew. Chemie Int. Ed. , 5420 (2006).[7] M. Guix, C. C. Mayorga-Martinez, and A. Merko¸ci,Chem. Rev. , 6285 (2014).[8] J. L. Anderson, M. E. Lowell, and D. C. Prieve, J. FluidMech. , 107 (1982).[9] J. Anderson, Annu. Rev. Fluid Mech. , 61 (1989). [10] J. F. Brady, J. Fluid Mech. , 216 (2011).[11] J. H. Irving and J. G. Kirkwood, J. Chem. Phys. , 817(1950).[12] P. Schofield and J. R. Henderson, Proc. R. Soc. A Math.Phys. Eng. Sci. , 231 (1982).[13] J. Kirkwood and F. Buff, J. Chem. Phys. , 338 (1949).[14] A. Harasima, J. Phys. Soc. Japan , 343 (1953).[15] J.-G. Weng, S. Park, J. R. Lukes, and C.-L. Tien, J.Chem. Phys. , 5917 (2000).[16] J. Cormier, J. M. Rickman, and T. J. Delph, J. Appl. Phys. , 99 (2001).[17] S. Plimpton, J. Comput. Phys. , 1 (1995).[18] Y. Liu, R. Ganti, H. G. A. Burton, X. Zhang, W. Wang,and D. Frenkel, Phys. Rev. Lett. , 224502 (2017).[19] R. Ganti, Y. Liu, and D. Frenkel, arXiv preprintarXiv:1710.01657 (2017).[20] P. Schofield and J. Henderson, Proc. R. Soc. A379