Entrance Effects in Concentration-Gradient-Driven Flow Through an Ultrathin Porous Membrane
EEntrance effects in concentration-gradient-driven flow through an ultrathin porousmembrane
Daniel J. Rankin, Lyd´eric Bocquet, and David M. Huang a)1) Department of Chemistry, School of Physical Sciences, The University of Adelaide,Australia Laboratoire de Physique Statistique, CNRS UMR 8550,Ecole Normale Sup´erieure, PSL Research University, Paris,France
Transport of liquid mixtures through porous membranes is central to processes such asdesalination, chemical separations and energy harvesting, with ultrathin membranesmade from novel 2D nanomaterials showing exceptional promise. Here we derive,for the first time, general equations for the solution and solute fluxes through acircular pore in an ultrathin planar membrane induced by a solute concentrationgradient. We show that the equations accurately capture the fluid fluxes measured infinite-element numerical simulations for weak solute–membrane interactions. We alsoderive scaling laws for these fluxes as a function of the pore size and the strength andrange of solute–membrane interactions. These scaling relationships differ markedlyfrom those for concentration-gradient-driven flow through a long cylindrical pore orfor flow induced by a pressure gradient or electric field through a pore in an ultrathinmembrane. These results have broad implications for transport of liquid mixturesthrough membranes with a thickness on the order of the characteristic pore size. a) Electronic mail: [email protected] a r X i v : . [ c ond - m a t . s o f t ] J un . INTRODUCTION Fluid transport through pores and porous membranes plays a key role in many processesof fundamental and practical interest, including cellular homeostasis in biological systems, chemical separations, desalination, and energy conversion. Thus, a general theoreticalunderstanding of the parameters that control these transport phenomena has broad im-plications for a variety of domains. Many theoretical models of fluid transport in porousmembranes have considered flows only within the pores and have neglected the effectof transport between the membrane pores and the fluid outside the membrane. These so-called entrance or access effects can dominate fluid transport processes when the membranethickness approaches the characteristic pore size or when the fluid–solid friction becomessmall.
The most extreme examples of this situation are membranes of atomic thicknessmade from 2D materials such as graphene and its derivatives or molybdenum sulfide.
Such 2D membranes are of great interest, as they have been shown to exhibit exceptionalproperties compared with conventional membranes for applications such as desalination and electrical energy harvesting from salinity gradients. Fluid fluxes across a membrane can be induced by a variety of driving forces, includinggradients of pressure, electrical potential, or solute concentration. Equations have previouslybeen derived to quantify entrance effects on fluid flow driven by a pressure gradient and on fluid flow and ionic electrical currents induced by an electric field acting ona electrolyte solution. However, to date, no theory has been developed to describe theentrance effects on fluid fluxes driven by concentration gradients and how they vary withrelevant parameters.Fluid fluxes driven by concentration gradients are of particular relevance in applicationssuch as chemical separations, desalination, and salinity-gradient-driven energyharvesting. The work presented here focuses specificially on entrance effects onthe concentration-gradient-driven process of diffusioosmosis, in which flow of a solute-containing solution is driven by an osmotic-pressure gradient that develops in the inhomo-geneous interfacial fluid layer induced by interactions of the fluid with the solid surfaces.Diffusioosmosis has been shown to play a key role in astonishing energy densities measuredfor salinity gradient energy harvesting in a nanotube membrane. Thus, entrance effects onthis phenomenon are of considerable interest.2ere we derive, for the first time, general equations to quantify the diffusioosmotic solu-tion flux and solute flux of a dilute solution through a circular aperture in a 2D membraneas a function of the aperture size and the strength and range of the interactions betweenthe solute and membrane surface. We verify the accuracy of the equations by compari-son with finite-element numerical simulations. We go on to compare the scaling behaviorpredicted for concentration-gradient-driven flow through a circular aperture with those forother membrane geometries and driving forces and discuss the implications of these resultsfor real systems.
II. THEORYA. Diffusioosmotic flow
Consider the flow of a solution containing a single solute type through a circular apertureof radius a in an infinitesimally thin planar wall, as illustrated in Fig. 1. Assuming that thefluid flows can be described by continuum hydrodynamic equations for low-Reynolds-numbersteady-state flow of a dilute solution of an incompressible Newtonian fluid, the governingequations are −∇ p − c ∇ U + η ∇ u = 0 , (1) ∇ · j = ∇ · ( − D ∇ c − µc ∇ U + u c ) = 0 , (2) ∇ · u = 0 , (3)where u , j , p , and c are the solution velocity, solute current density, pressure, and soluteconcentration, respectively, η is the solution shear viscosity, and U is the solute–membraneinteraction potential. D and µ are the solute diffusivity and mobility respectively, which weassume are related by the Einstein relation, µ = Dk B T , where k B is the Boltzmann constantand T is the temperature. U is the interaction potential per solute molecule and so − c ∇ U is a body force per unit volume acting on the fluid due to the solute–membrane interactions. U is assumed to depend on the position in the fluid relative to the membrane surface. For aneutral solute, U typically depends on the distance from the surface and has a range onthe order of the solute molecular diameter. Further assuming that advection of the solute isnegligible compared with diffusion (i.e. low P´eclet number flow), Eq. (2) for the solute flux3 IG. 1. Schematic of flow of a solution through a circular aperture of radius a in an infinitesimallythin planar wall. The origin in cylindrical ( r, z ) coordinates is at the center of the aperture, asindicated, and the axis of symmetry is the depicted z axis. The solute concentration c and solutionpressure p far from the membrane are c H and p H , respectively, in the upper half-plane and c L and p L , respectively, in the lower half-plane. Q and J are the solution and solute fluxes, respectively.Contours of constant ζ and ν in oblate–spheroidal ( ζ, ν ) coordinates are also shown as dashed lines,with unit vectors shown at one point in space. simplifies to ∇ · j = ∇ · (cid:18) − D ∇ c − D ck B T ∇ U (cid:19) = 0 . (4)The solution velocity and solute flux are assumed to satisfy the usual no-slip and zero fluxboundary conditions at the membrane surface, i.e. u = 0 and ˆ n · j = 0, where ˆ n is the unitvector normal to the membrane surface.We note that a similar approach based on the widely used Poisson–Nernst–Planck–Stokesequations for electrolytes, in which the electric potential energy plays an analogous role forthe electrolyte that the interaction potential U does for a neutral solute, could be used toextend this study to concentration-gradient-driven electrolyte transport. However, such anextension is non-trivial, as the electric potential must be determined by solving an additionalcoupled differential equation (the Poisson equation) that depends on the solute (electrolyte)concentration, rather than being specified. Thus, we leave this extension to electrolytes tofuture work. 4onsider the fluid flow induced by a concentration difference, ∆ c = c H − c L , between thetwo sides of the membrane, with the pressure far from the membrane the same on both sidesof the membrane, i.e. p H = p L = p ∞ , as shown in Fig. 1.Our derivation uses a combination of cylindrical ( r, z, φ ) and oblate spheroidal ( ζ, ν, φ )coordinates, where z = aνζ , r = a p (1 + ν )(1 − ζ ), and θ is the angle about the z axis(0 ≤ ν < ∞ , − ≤ ζ ≤
1, and 0 ≤ φ < π ). Writing the solute concentration in thepresence of a concentration gradient as c ( ζ, ν ) ≡ c s ( ζ, ν ) e − U ( ζ,ν ) /k B T , (5)where c s ( ζ, ν ) is to be determined, and inserting this expression into Eq. (4) and the bound-ary condition for the solute current density j gives ∇ c s − ∇ (cid:18) Uk B T (cid:19) · ∇ c s = 0 , (6)with the boundary condition ˆ n · ∇ c s = 0 at the membrane surface. If the solute–membraneinteraction potential is small relative to the thermal energy ( U (cid:28) k B T ), Eq. (6) reduces to ∇ c s = 0 . (7)Solving this equation, subject to the boundary conditions on c s at the membrane surfaceand far from the membrane ( c s → c H = c ∞ + ∆ c and c L = c ∞ − ∆ c in the upper and lowerhalf-planes, respectively, where U → c s = c ∞ + ∆ cπ tan − ν. (8)We have verified using finite-element numerical simulations (see Fig. S3 of the supplementarymaterial) that Eq. (8) with Eq. (5) accurately describe the solute concentration even when U is several times k B T . A possible reason why Eq. (8) appears to be accurate outside theregime for which it was derived is that the second term in Eq. (6) that was neglected toarrive at Eq. (7) and (8) can be small even if the magnitude of the potential U is large: forexample, for c s given by Eq. (8) this term is zero for a potential that is a function only of ζ or in the pore mouth (at ν = z = 0) for a potential that is a function only of distance fromthe membrane, due to the orthogonality of ∇ U and ∇ c s in these cases.The fluid flow induced by the concentration gradient can be obtained from the reciprocaltheorem for steady incompressible creeping flow, which allows the fluid flow due to a5ody force F acting on the fluid to be related to the pressure-driven flow in the same poregeometry, for which an analytical solution exists for the fluid velocity through a circularaperture. As shown by Mao et al. for the related problem of electroosmosis through acircular aperture, Q = − p Z Z Z V d V ¯ u · F , (9)where ¯ u is the fluid velocity induced by a pressure difference ∆ p = p H − p L for ∆ c = 0 in thesystem geometry in Fig. 1 and the integral is over the volume V occupied by the fluid. Forconcentration-gradient-driven flow described by Eq. (1), F = − c ∇ U . The pressure-drivenflow velocity can be obtained from the stream function ψ = − a πη (1 − ζ ) ∆ p for the flow using ¯ u = r ˆ φ × ∇ ψ , where ˆ φ is the unit vector in the φ direction, as¯ u = − aζ ∆ p πη p (1 + ν )( ν + ζ ) ˆ ν . (10)Inserting this expression for ¯ u and F = − c ∇ U into Eq. (9) (with c given by Eqs. (5) and (8))and making use of ∂c s ∂ν = ∆ cπ (1 + ν ) − yields Q = − κ DO ∆ c, (11)with κ DO = − k B T a πη Z d ζ ζ Z ∞ d ν (cid:18) e − U/k B T −
11 + ν (cid:19) . (12)Equation (12) is the main result of this work.Furthermore, the solute flux density can be obtained, using Eqs. (4) and (5), as j = − De − U/k B T ∇ c s . (13)The solute flux across the membrane is J = Z Z S d s j · ˆ n , (14)where the unit vector normal to the pore mouth is ˆ n = ˆ ν = ˆ z and the surface integral isover the pore aperture. Evaluating the solute flux at the pore mouth (at ν = z = 0) usingEqs (8) and (13) yields J = − Da ∆ c Z d ζ e − U/k B T , (15)= − D ∆ c Z a d r re − U/k B T √ a − r . (16)6 . Limiting cases and scaling behavior The diffusioosmotic mobility predicted in Eq. (12) depends crucially on the range ofthe interaction potential U , which we will call λ . However, the term exp( − U/k B T ) − ζ and ν . As a consequence, the mobility,and its scaling with the pore radius a and interaction range λ , may depend on specific detailsof the geometry dependence of the interaction. Therefore, we consider various limiting casesfor the geometry of the potential and the consequences for the dependence of the scalingwith the pore radius and interaction range.
1. Case of potential that depends only on the distance to the membranesurface and/or pore edge
In most circumstances, the potential U is expected to be a function of the distance d from the membrane surface. However, the integral in Eq. (12) for the mobility cannot besimplified due to the geometrical interplay between the various variables.Let us assume for simplicity that the solute excess/depletion at the membrane surfacecan be represented by a step function as a function of the distance d from the surface, i.e. e − U ( d ) /k B T = α, d ≤ λ , d > λ , (17)where α characterizes the solute excess close to the membrane surface ( α > α < λ (cid:29) a ), e − U ( d ) /k B T canbe approximated as a constant independent of the coordinates, and we find from Eq. (12)in this limit that κ DO ’ − k B T ∆ c η ( α − a , λ (cid:29) a. (18)On the other hand, if λ (cid:28) a , the integral in Eq. (12) can be approximated as κ DO ’ − k B Tη (cid:18)
14 + 1 π (cid:19) ( α − a λ . (19)Details of the derivation of this equation can be found in the supplementary material.7 similar calculation can be performed for the case in which the interaction originatesonly from the pore edge, and U depends on the distance to the edge. As detailed in thesupplementary material, the mobility in this case is κ DO ’ − k B T η ( α − a λ . (20) Solute flux.
Similarly, using Eq. (17) in Eq. (16) and noting that d = a − r in the poremouth at z = 0, the solute flux across the membrane for the step-function potential can beshown to be J step = − D ∆ c h a + ( α − p λ (2 a − λ ) i (21)for λ ≤ a and J step ’ − D ∆ cαa (22)for λ (cid:29) a .
2. Case of potential that depends only on ζ In the case in which the potential U is a function of the ζ coordinate only (see Fig. 1),the diffusioosmotic mobility in Eq. (12) simplifies to κ DO = − k B T a ∆ cη Z d ζ ζ (cid:0) e − U ( ζ ) /k B T − (cid:1) . (23)Assuming that U only depends on ζ is a stringent condition in terms of symmetry, butinterestingly such a result is expected for a circular aperture in a planar membrane at afixed electrostatic potential (however, in the absence of screening).To determine how Q scales with λ , particularly in the limit λ (cid:28) a , we can use therelationship between ζ and the distances, d and d , from a point with this ζ value to thetwo foci of the hyperbolae or ellipses of constant ζ or ν that are shown in Fig. 1. For aninfinitesimally thin membrane, these foci are located at the pore edge (at r = a ), and thus1 − ζ = [( d − d ) / (2 a )] . So a potential U ( ζ ) that depends only on ζ depends only onthe relative distance, d − d . Furthermore, assuming that the potential has a distancerange λ implies that U ( ζ ) = U [( a − | d − d | / /λ ]. In addition, U is only non-zero when a − | d − d | / ∼ λ (cid:28) a , so a − | d − d | / a (1 − p − ζ ) ’ aζ / U ( ζ ) ∼ U [ aζ / (2 λ )] and the integral in Eq. (23)8ecomes Z d ζ ζ (cid:16) e − U [ aζ / (2 λ ) ] /k B T − (cid:17) ’ (2 λ/a ) / Z ∞ d x x (cid:16) e − U ( x ) /k B T − (cid:17) , (24)where x ≡ [ a/ (2 λ )] / ζ . So the mobility in this case scales as κ DO ∼ a ( λ/a ) / ∼ ( aλ ) / .Equation (23) can also be rewritten in radial coordinates, by focusing on the pore mouth( ν = z = 0), as κ DO = − k B Tη Z a d rr √ a − r (cid:0) e − U/k B T − (cid:1) , (25)where r is the distance to the center of the pore in the membrane plane ( z = 0). As analternative approach to predicting the dependence of the mobility on λ , we again restrictourselves for simplicity to a step-function interaction versus distance d = a − r from thepore edge, i.e. e − U ( r ) /k B T = α, a − r ≤ λ , a − r > λ . (26)Inserting Eq. (26) into Eq. (25) gives the diffusioosmotic mobility for the step-functionpotential as κ DO = − k B T η ( α −
1) [ λ (2 a − λ ))] / . (27)For a small interaction range λ compared with the pore radius, the predicted scaling of themobility is κ DO ’ − k B T η ( α −
1) ( aλ ) / , (28)which is identical to the scaling predicted directly from Eq. (23).
3. Case of potential that depends only on ν In the case in which the potential U is a function of the ν coordinate only (see Fig. 1),the diffusioosmotic mobility in Eq. (12) simplifies to κ DO = − k B T a πη Z ∞ d ν (cid:18) e − U ( ν ) /k B T −
11 + ν (cid:19) . (29)Following similar reasoning to the previous section, in terms of the distances, d and d ,from a point with a given ν value to the foci of the ellipses or hyperbolae in Fig. 1, 1 + ν =[( d + d ) / (2 a )] . So a potential U ( ν ) that depends only on ν depends only on the averagedistance, ( d + d ) /
2. Therefore, assuming that the potential has a distance range λ entails9 ( ν ) = U ([( d + d ) / − a ] /λ ). For λ (cid:28) a , ( d + d ) / − a = a ( √ ν − ’ aν / U ( ν ) ∼ U [ aν / (2 λ )] and the integral in Eq. (29) becomes Z ∞ d ν e − U [ aν / (2 λ )] /k B T −
11 + ν ! = (2 λ/a ) / Z ∞ d x e − U ( x ) /k B T −
11 + (2 λ/a ) x ! , (30) ’ (2 λ/a ) / Z ∞ d x (cid:16) e − U ( x ) /k B T − (cid:17) , (31)where x ≡ [ a/ (2 λ )] / ν . So the mobility in this case scales as κ DO ∼ a ( λ/a ) / ∼ a / λ / . III. NUMERICAL RESULTS
To validate the theory, finite-element method (FEM) simulations of concentration-gradient-driven flow were carried out using Comsol Multiphysics (version 4.3a) throughpores with various radii and solute–membrane interactions. Here we consider that the soluteinteracts with the membrane via a potential that depends on the (shortest) distance d to themembrane surface. Accordingly, the solute–membrane interaction potential was modelledusing a hyperbolic tangent function, U ( d ) = (cid:15) (cid:20) − tanh (cid:18) d − λλ (cid:19)(cid:21) , (32)defined by parameters (cid:15) and λ describing the strength and range of the potential. In allsimulations, the average solute concentration c ∞ was 10 − σ − and the solute diffusivity D was σ /τ , and unless otherwise stated the aperture radius a was 10 σ and λ was σ ,where σ is the unit of length ( σ can be regarded as the diameter of a fluid molecule) and τ = ησ / ( k B T ) is the unit of time. Details of the FEM simulations, which all correspond tolow-P´eclet number flow, are given in the supplementary material.We have quantified the concentration-gradient-driven solution and solute fluxes measuredin the numerical simulations and predicted by our theory in terms of the diffusioosmoticmobility κ DO defined by κ DO ≡ − Q ∆ c (33)and the solute permeance P s defined by P s ≡ − J ∆ c . (34)The equations that we have derived for the solution and solute fluxes (Eqs. (12) and (15))predict that the fluxes are linearly related to the concentration difference ∆ c and thus that10he conductances and resistances defined in Eqs. (33) and (34) are independent of ∆ c . Wehave verified that this is indeed the case for the range of concentration differences studiedin the numerical simulations (∆ c = 10 − to 3 × − σ − = 10 − to 0 . c ∞ ), as shown inFig. S5 of the supplementary material.Figure 2 shows the diffusioosmotic mobility κ DO from the simulations and theory of acircular aperture as a function of aperture radius a and solute–membrane interaction range λ for two different values of the solute–membrane interaction strength (cid:15) , k B T /
10 and k B T ,with all other parameters kept constant. The theory curves were calculated by numericallyevaluating the integral in Eq. (12) with the solute–membrane potential U in Eq. (32). Thesign of κ DO has been defined so that a positive and negative values correspond to fluidflow in opposite and same direction, respectively, to the applied concentration gradient.Hence, for solute depletion at the membrane surface ( (cid:15) >
0) the flow is towards highersolute concentration ( κ DO < κ DO > Figure 2 shows good quantitative agreement between the theory and simulations for thevariation of κ DO with all relevant parameters for (cid:15) . k B T . Since the theory assumes a weakpotential in deriving Eq. (8) for the solute concentration, we indeed find deviations betweenthe prediction and the simulations as the magnitude of the solute–membrane potential in-creases. Nevertheless, the agreement is reasonable well beyond the regime of validity of thisapproximation. Note that for the values of (cid:15) in Fig. 2, α − ≈ e − (cid:15)/k B T − ≈ − (cid:15)/k B T ,and so from Eq. (18) or (19) κ DO is approximately proportional to (cid:15) , but this scaling is notexpected in general and already starts to break down for (cid:15) = k B T .Figure 2 also compares the simulation results with the approximate scaling with the poreradius a and solute–membrane interaction range λ predicted by the theory. For a (cid:28) λ ,Eq. (18) predicts that κ DO is proportional to a and independent of λ , which is evident inthe scaling for small a in Fig. 2(a) and in the saturation at large λ in Fig. 2(b). On theother hand, for a (cid:28) λ , Eq. (19) predicts scaling of κ DO with aλ , which is seen to hold inthe large- a regime in Fig. 2(a) and in the small- λ regime in Fig. 2(b).Figure 3 shows the analogous comparison between the FEM simulations and theory forthe solute permeance P s . The theory curves were calculated by numerically evaluatingthe integral in Eq. (16) with the solute–membrane potential U in Eq. (32). As for thediffusioosmotic conductance, the theory accurately captures the simulated solute permeance11 IG. 2. Diffusioosmotic mobility κ DO versus (a) pore radius a (with λ = σ ) and (b) solute–membrane interaction range λ (with a = 10 σ ) for solute–membrane interaction strength (cid:15) = k B T /
10 or k B T from FEM simulations (points) and theory (solid lines). The dashed lines showscaling with various powers of a and λ . for (cid:15) . k B T , with deviations between the theory and simulations becoming evident formagnitudes of the solute–membrane potential greater than k B T . For the parameters used inFig. 3(a), the first term in Eq. (21) dominates and so the permeance P s shows the expectedlinear scaling with the pore radius a . For λ (cid:28) a , Eq. (21) predicts that P s varies from itsvalue at λ = 0 with a scaling as √ λ , which is evident in Fig. 3(b).12 IG. 3. Solute permeance P s versus (a) pore radius a (with λ = σ ) and (b) solute–membraneinteraction range λ (with a = 10 σ ) for solute–membrane interaction strength (cid:15) = k B T /
10 or k B T from FEM simulations (points) and theory (solid lines). IV. DISCUSSION
An interesting outcome of the previous results is that the diffusioosmotic mobility κ DO = − Q/ ∆ c across a circular aperture in ultrathin membrane is strongly dependent on the detailsof the interaction of the solute with the membrane. We have shown in particular that themobility scales with the pore radius a and interaction range λ as a − γ λ γ , with an exponent γ that depends on the underlying symmetries of the potential: when λ (cid:28) a , for a potentialthat depends only on the ζ coordinate (in the oblate-spheroidal system, see Fig. 1) anexponent γ = 3 / ν coordinate an13 ABLE I. Scaling of the diffusioomotic mobility κ DO with pore radius a and solute–membraneinteraction range λ for a circular aperture in a 2D membrane, as well as with pore length L fora long cylindrical pore, for different functional forms of the solute–membrane potential U in thelimits λ (cid:28) a and λ (cid:29) a . System Limit Potential ScalingCircular aperture λ (cid:28) a U ( d ) aλ λ (cid:28) a U ( ζ ) a / λ / λ (cid:28) a U ( ν ) a / λ / λ (cid:29) a any a Cylinder λ (cid:28) a U ( d ) a λ /Lλ (cid:29) a any a /L exponent γ = 1 / γ = 2 is found.It is also interesting to compare the results for the circular aperture with those obtainedin long cylindrical pores (e.g. as a model for nanotubes). As derived in detail in thesupplementary material, the diffusioosmotic mobility of a long cylindrical pore of length L is proportional to ( aλ ) /L for λ (cid:28) a . When compared to the case leading to an exponent γ =2, the scaling of κ DO with pore size and interaction range λ is therefore recovered by replacingthe length of the nanopore L with the pore size a , which is indeed expected for entranceeffects. However, as shown with the case leading to an exponent γ = 3 / / κ DO for a circular aperture and a long cylindrical pore are summarized in Table I.These result are relevant for transport through finite-length pores, for which the totalresistance to flow can often be accurately given by the sum of the resistance due to the poreinterior and that due to the pore ends, which can be approximated by that of a circularaperture. The predicted scaling behavior of the diffusioosmotic mobility and solute permeance ofa circular aperture for concentration-gradient-driven flow differ markedly from the scalingbehavior derived previously for other types of flows in the same system geometry. Forexample, the hydraulic conductance (solution flux per unit pressure difference) in pressure-14riven fluid flow through a circular aperture has been shown to be proportional to a in contrast to the proportionalty with a − γ where γ > λ (cid:28) a .On the other hand, for λ (cid:29) a , the diffusioosmotic mobility shows the same a scal-ing as the hydraulic conductance. The equivalent scaling of the hydraulic conductanceand diffusioosmotic mobility in this limit is because concentration-driven-flow in which thesolute–membrane interaction range is larger than the aperture radius is equivalent to os-motic transport through a semipermeable membrane, with the osmotic pressure gradientdue to the concentration gradient playing an equivalent role to the pressure gradient inpressure-driven flow. Likewise, the scaling behavior predicted here for the diffusioosmoticmobility differs from the electroosmotic conductance for electric-field-driven fluid flow of anelectrolyte, which has been shown to be proportional to a /λ D for a (cid:28) λ D , where λ D isthe Debye length characterizing the electric double layer width, equivalent to λ here; in thesame limit, scaling of the diffusioosmotic mobility with a is predicted here.The electrical conductance across a circular aperture in electric-field-driven transportof an electrolyte has been shown to be proportional to a for an uncharged membrane, with the addition of surface charge to the membrane only changing the length scale in thescaling relationship to an effective radius a eff that is the sum of the aperture radius andthe Dukhin length characterizing the ratio of the surface electrical conductivity to the bulkelectrical conductivity, without changing the scaling exponent. This scaling differs fromthat derived for the solute permeance for λ ≤ a , for which the equivalent effective radius is a eff = h a + ( α − p λ (2 a − λ ) i , in which the second term in the sum depends on both thesolute–membrane interaction range λ and pore radius a .We can consider the implications of our theory for realistic systems, in particular fordiffusioosmotic transport of electrolyte solutions. Extending the theory to electrolytesis desirable for applications, such as salinity-gradient-driven energy conversion anddesalination, but it is technically difficult, so we leave this derivation for future work.However, as a rule of thumb, one may expect that this would amount to replacing the solute–membrane interaction range λ by the salt-concentration-dependent Debye length λ D ∝ c − / .A counterintuitve outcome of the non-universal dependence of the diffusioosmotic mobilityas a function of the interaction range λ is a possible impact on the dependence of the mo-bility on the salt concentration. In a long pore under a salinity difference ∆ c salt , the solvent15ux is predicted to behave as Q = K DO ∆ log c salt , with the mobility K DO ∝ c πa /L fora long pore of length L , i.e. no dependence on salt concentration. This is due to thedependence of the diffusioosmotic mobility κ DO ∼ λ D ∼ c − , so that Q ∼ ∆ log c salt . Now,for a circular aperture, we have shown that the dependence of the diffusioosmotic mobility κ DO on the interaction range scales as κ DO ∼ λ γ with a non-universal exponent γ . Thus, thecase of γ = 2 (occurring for potentials that depend on the distance to the membrane) willlead to K DO ∝ c as for the long-pore case; but cases with γ = 2 (as highlighted above intwo cases), will lead to K DO ∝ c − γ/ , exhibiting therefore a curious dependence on c salt . V. CONCLUSIONS
In summary, we have derived general equations and scaling relationships as a function ofpore radius and solute–membrane interaction strength and range for the solution and solutefluxes induced by a solute concentration gradient through a circular aperture in an ultrathinplanar membrane. We have shown, by comparing with finite-element numerical simulations,that the equations accurately quantify the fluid fluxes when the solute–membrane interactionstrength is small compared with the thermal energy k B T . In the limit of a solute–membraneinteraction range much smaller than the pore radius, the theory predicts a non-universaldependence of the fluid fluxes on the pore radius and interaction range. These results havesignificant implications for applications involving concencentration-gradient-driven flow inmembranes in which the thickness is on the order of the pore size, such as those made from2D nanomaterials, notably in the context of blue energy harvesting. SUPPLEMENTARY MATERIAL
The supplementary material contains derivations of scaling laws for diffusioosmosisthrough a circular aperture in an ultrathin planar membrane, the theory of concentration-gradient-driven flow in a long cylindrical pore, and further details and supplementary resultsof finite-element numerical simulations. 16
CKNOWLEDGMENTS
D.J.R. acknowledges the support of an Australian Government Research Training Pro-gram Scolarship and a University of Adelaide Faculty of Sciences Divisional Scholarship.L.B. acknowledges funding from the EU H2020 Framework Programme/ERC AdvancedGrant agreement number 785911–Shadoks and ANR project Neptune.
REFERENCES H. X. Sui, B. G. Han, J. K. Lee, P. Walian, and B. K. Jap, “Structural basis of water-specific transport through the AQP1 water channel,” Nature , 872–878 (2001). J. R. Werber, C. O. Osuji, and M. Elimelech, “Materials for next-generation desalinationand water purification membranes,” Nat. Rev. Mater. , 16018 (2016). M. Elimelech and W. A. Phillip, “The future of seawater desalination: Energy, technology,and the environment,” Science , 712–717 (2011). B. E. Logan and M. Elimelech, “Membrane-based processes for sustainable power genera-tion using water,” Nature , 313–319 (2012). A. Siria, M.-L. Bocquet, and L. Bocquet, “New avenues for the large-scale harvesting ofblue energy,” Nat. Rev. Chem. , 0091 (2017). J. C. Fair and J. F. Osterle, “Reverse electrodialysis in charged capillary membranes,” J.Chem. Phys. , 3307–3316 (1971). D. J. Bonthuis, K. F. Rinne, K. Falk, C. N. Kaplan, D. Horinek, A. N. Berker, L. Boc-quet, and R. R. Netz, “Theory and simulations of water flow through carbon nanotubes:prospects and pitfalls,” J. Phys.: Condens. Matter , 184110 (2011). S. Balme, F. Picaud, M. Manghi, J. Palmeri, M. Bechelany, S. Cabello-Aguilar, A. Abou-Chaaya, P. Miele, E. Balanzat, and J. M. Janot, “Ionic transport through sub-10 nmdiameter hydrophobic high-aspect ratio nanopores: experiment, theory and simulation,”Sci. Rep. , 10135 (2015). P. B. Peters, R. van Roij, M. Z. Bazant, and P. M. Biesheuvel, “Analysis of electrolytetransport through charged nanopores,” Phys. Rev. E , 053108 (2016). J. D. Sherwood, M. Mao, and S. Ghosal, “Electroosmosis in a finite cylindrical pore:Simple models of end effects,” Langmuir , 9261–9272 (2014).17 D. V. Melnikov, Z. K. Hulings, and M. E. Gracheva, “Electro-osmotic flow throughnanopores in thin and ultrathin membranes,” Phys. Rev. E , 063105 (2017). T. Sisan and S. Lichter, “The end of nanochannels,” Microfluid. Nanofluid. , 787–791(2011). D. J. Rankin and D. M. Huang, “The effect of hydrodynamic slip on membrane-basedsalinity-gradient-driven energy harvesting,” Langmuir , 3420–3432 (2016). S. P. Surwade, S. N. Smirnov, I. V. Vlassiouk, R. R. Unocic, G. M. Veith, S. Dai, and S. M.Mahurin, “Water desalination using nanoporous single-layer graphene,” Nat. Nanotechnol. , 459–464 (2015). D. Cohen-Tanugi, L.-C. Lin, and J. C. Grossman, “Multilayer nanoporous graphene mem-branes for water desalination,” Nano Lett. , 1027–1033 (2016). A. Morelos-Gomez, R. Cruz-Silva, H. Muramatsu, J. Ortiz-Medina, T. Araki, T. Fukuyo,S. Tejima, K. Takeuchi, T. Hayashi, M. Terrones, and M. Endo, “Effective NaCl and dyerejection of hybrid graphene oxide/graphene layered membranes,” Nat. Nanotechnol. ,1083–1088 (2017). X. Li, W. Xu, M. Tang, L. Zhou, B. Zhu, S. Zhu, and J. Zhu, “Graphene oxide-basedefficient and scalable solar desalination under one sun with a confined 2D water path,”Proc. Natl. Acad. Sci USA , 13953–13958 (2016). M. I. Walker, K. Ubych, V. Saraswat, E. A. Chalklen, P. Braeuninger-Weimer, S. Caneva,R. S. Weatherup, S. Hofmann, and U. F. Keyser, “Extrinsic cation selectivity of 2Dmembranes,” ACS Nano , 1340–1346 (2017). L. Wang, M. S. H. Boutilier, P. R. Kidambi, D. Jang, N. G. Hadjiconstantinou, andR. Karnik, “Fundamental transport mechanisms, fabrication and potential applications ofnanoporous atomically thin membranes,” Nat. Nanotechnol. , 509 (2017). J. Feng, M. Graf, K. Liu, D. Ovchinnikov, D. Dumcenco, M. Heiranian, V. Nandigana,N. R. Aluru, A. Kis, and A. Radenovic, “Single-layer MoS nanopores as nanopowergenerators,” Nature , 197–200 (2016). W. Li, Y. Yang, J. K. Weber, G. Zhang, and R. Zhou, “Tunable, strain-controllednanoporous MoS filter for water desalination,” ACS Nano , 1829–1835 (2016). R. A. Sampson, “On Stokes’s current function,” Philos. Trans. R. Soc. A , 449–518(1891). 18 H. L. Weissberg, “End correction for slow viscous flow through long tubes,” Phys. Fluids , 1033–1036 (1962). M. Mao, J. Sherwood, and S. Ghosal, “Electro-osmotic flow through a nanopore,” J. FluidMech. , 167–183 (2014). J. E. Hall, “Access resistance of a small circular pore.” J. Gen. Physiol. , 531–532 (1975). C. Lee, L. Joly, A. Siria, A.-L. Biance, R. Fulcrand, and L. Bocquet, “Large apparentelectric size of solid-state nanopores due to spatially extended surface conduction,” NanoLett. , 4037–4044 (2012). A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, and L. Boc-quet, “Giant osmotic energy conversion measured in a single transmembrane boron nitridenanotube,” Nature , 455–458 (2013). J. L. Anderson, “Colloid transport by interfacial forces,” Annu. Rev. Fluid Mech. ,61–99 (1989). J. L. Anderson, M. E. Lowell, and D. C. Prieve, “Motion of a particle generated bychemical gradients Part 1. non-electrolytes,” J. Fluid Mech. , 107–121 (1982). R. F. Probstein,
Physicochemical Hydrodynamics: An Introduction , 2nd ed. (Wiley-Interscience, Hoboken, 1994). P. M. Morse and H. Feshbach, in
Methods of Theoretical Physics , Vol. 2 (McGraw-HillBook Company, 1953) p. 1292. J. Happel and H. Brenner,
Low Reynolds number hydrodynamics (Martinus Nihjoff Pub-lishers, 1983). D. M. Huang, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, “Massive amplification ofsurface-induced transport at superhydrophobic surfaces,” Phys. Rev. Lett. , 064503(2008). S. Marbach, H. Yoshida, and L. Bocquet, “Osmotic and diffusio-osmotic flow generation athigh solute concentration. i. mechanical approaches,” J. Chem. Phys. , 194701 (2017).19 upplementary Material: Entrance effects in concentration-gradient-driven flowthrough an ultrathin porous membrane
Daniel J. Rankin, Lyd´eric Bocquet, and David M. Huang a)1) Department of Chemistry, School of Physical Sciences, The University of Adelaide,Australia Laboratoire de Physique Statistique, CNRS UMR 8550,Ecole Normale Sup´erieure, PSL Research University, Paris,France a) Electronic mail: [email protected] S1 a r X i v : . [ c ond - m a t . s o f t ] J un
1. DERIVATION OF SCALING LAWS FOR DIFFUSIOOSMOSISTHROUGH A CIRCULAR APERTURE IN A 2D PLANAR MEMBRANE
Here we derive scaling laws for the diffusioosmotic flux through a circular aperture in aninfinitesimally thin planar wall in the limit that the solute–membrane interaction range λ ismuch smaller than the aperture radius a , i.e. λ (cid:28) a . To derive analytical expressions, weassume a step-function form for the potential given by Eq. (17) of the main paper. Insertingthis equation into Eq. (12) of the main paper for the diffusioosmotic flux gives Q = − k B T a ∆ cπη ( α − Z d ζ ζ Z ∞ d ν H [ λ − d ( ζ, ν )]1 + ν , (S1)where α , as discussed in the main paper, is a parameter characterizing the solute excessclose the membrane surface, d ≥ H is the Heavisidestep function, which sets the condition that the integrand is only non-zero when 0 ≤ d ≤ λ . A. Solute interaction with entire membrane surface
Using the relationships between the cylindrial ( r, z ) and oblate–spheroidal ( ζ, ν ) coordi-nates, r = a p (1 + ν )(1 − ζ ) and z = aνζ , gives d = p ( a − r ) + z , r ≤ az, r > a , (S2)= a q ν − ζ − p (1 + ν )(1 − ζ ) , p (1 + ν )(1 − ζ ) ≤ aνζ, p (1 + ν )(1 − ζ ) > . (S3)Since the integrand in Eq. (S1) is only non-zero when d ≤ λ , we only need to consider thecase when this condition is satisfied to evaluate Eq. (S1). If d ≤ λ (cid:28) a , then ζ (cid:28) ν (cid:28) p (1 + ν )(1 − ζ ) ≤
1. In this case, using √ ν = 1 + ν / − ν / . . . and p − ζ = 1 − ζ / − ζ / − . . . , Eq. (S3) simplifies to d ’ a ( ν + ζ ) , ν ≤ ζaνζ, ν > ζ . (S4)S2rom this equation, for d (cid:28) a , Z d ζ ζ Z ∞ d ν H [ λ − d ( ζ, ν )]1 + ν ’ Z d ζ ζ (cid:26)Z ζ d ν H [ λ − a ( ν + ζ ) / ν + Z ∞ ζ d ν H ( λ − aνζ )1 + ν (cid:27) , (S5)= Z √ λ/a d ζ ζ Z ζ d ν ν + Z √ λ/a √ λ/a d ζ ζ Z √ λ/a − ζ d ν ν + Z √ λ/a d ζ ζ Z λaζ ζ d ν ν , (S6)= Z √ λ/a d ζ ζ tan − (cid:18) λaζ (cid:19) + Z √ λ/a √ λ/a d ζ ζ tan − (cid:16)p λ/a − ζ (cid:17) , (S7)= − (cid:20) ( a + 3 λ ) π + 2 λa (cid:21) + π (cid:18) a + 2 λa (cid:19) / + 16 (cid:18) λa (cid:19) − (cid:18) a + 2 λa (cid:19) / tan − (cid:18)r aa + 2 λ (cid:19) − (cid:18) λa (cid:19) ln (cid:16) aλ (cid:17) . (S8)Inserting Eq. (S8) into Eq. (S1) yields Q ’ − k B T ∆ cη ( α − ( − [( a + 3 λ ) π + 2 λ ] a π + ( a + 2 λ ) / a / π aλ − π ( a + 2 λ ) / a / tan − (cid:18)r aa + 2 λ (cid:19) − π λ ln (cid:16) aλ (cid:17)(cid:27) , (S9)= − k B T ∆ cη ( α − a ((cid:18)
14 + 1 π (cid:19) (cid:18) λa (cid:19) + (cid:18) − − π −
12 ln a + 12 ln λ π (cid:19) (cid:18) λa (cid:19) + O "(cid:18) λa (cid:19) . (S10)Thus, for λ (cid:28) a , Q ’ − k B T ∆ cη (cid:18)
14 + 1 π (cid:19) ( α − aλ . (S11) B. Solute interaction with pore mouth only
Assuming that the solute interacts only with the part of the membrane surface at the poremouth, the solute–membrane potential U ( d ) will only depend on the distance from the sur-face at radial coordinate r = a . In this case, the distance is d = p ( a − r ) + z , ∀ r, z . Thus,S3ollowing analogous steps to those in the previous section, if d ≤ λ (cid:28) a , d ’ a ( ν + ζ ) and Z d ζ ζ Z ∞ d ν H [ λ − d ( ζ, ν )]1 + ν ’ Z d ζ ζ Z ∞ d ν H [ λ − a ( ν + ζ ) / ν , (S12)= Z √ λ/a d ζ ζ Z √ λ/a − ζ d ν ν , (S13)= Z √ λ/a d ζ ζ tan − (cid:16)p λ/a − ζ (cid:17) (S14)= π "(cid:18) a + 2 λa (cid:19) / − (cid:18) a + 3 λa (cid:19) . (S15)Inserting Eq. (S15) into Eq. (S1) yields Q ’ − k B T ∆ c η ( α − h ( a + 2 λ ) / a / − ( a + 3 λ ) a i , (S16)= − k B T ∆ cη ( α − a ( (cid:18) λa (cid:19) − (cid:18) λa (cid:19) + O "(cid:18) λa (cid:19) . (S17)Thus, for λ (cid:28) a , Q ’ − k B T ∆ c η ( α − aλ . (S18)Thus, the scaling of Q with the pore radius a and solute–membrane interaction range λ andstrength (measured by α ) is the same as in the case of solute interactions with the entiremembrane. S2. THEORY OF CONCENTRATION-GRADIENT-DRIVEN FLOW IN ALONG CYLINDRICAL PORE
Here we derive equations for concentration-gradient-driven flow of a dilute solution in along cylindrical pore (Fig. S1) assuming the same governing equations used in the main textfor flow through a circular aperture in an ultrathin planar membrane. We assume that thepore length L is much larger than the pore radius ( L (cid:29) a ), so that entrance effects canbe ignored. The solute–wall interaction potential U ( r ) in this case is only a function of theradial coordinate r .Since L (cid:29) a , the radial solute flux will be small relative to the axial solute flux (0 ≈ j r (cid:28) j z ), so from Eq. (4) of the main paper, − ∂c∂r − ck B T ∂U∂r = 0 , (S19)S4IG. S1: Cylindrical pore geometry for deriving fluxes through a long pore. The pore hasa radius of a and a length of L . The dashed line indicates the axis of symmetry.which can be integrated to give c = c s ( z ) e − U/k B T , (S20)where c s ( z ) is the solute concentration where U = 0. Similarly, the radial solution velocitywill be negligible compared with the axial solution velocity (0 ≈ u r (cid:28) u z ), so from the r -component of Eq. (1) of the main paper, ∂p∂r + c ∂U∂r ≈ , (S21)which can be integrated to give the solution pressure as p = p ∞ + k B T c s ( z ) (cid:0) e − U ( r ) /k B T − (cid:1) , (S22)where p ∞ , a constant, is the pressure where U = 0. Inserting the z -derivative of Eq. (S22)into the z -component of Eq. (1) in the main paper and integrating twice using the no-slipboundary condition ( u = 0) at the pore surface gives the axial velocity as u z ( r ) = − k B Tη d c s d z Z ar d r r Z r d r r (cid:16) e − U ( r ) /k B T − (cid:17) . (S23)The solution flux at any cross-section of the pore is Q = Z Z S d s u · ˆ n , (S24)into which Eq. (S23) can be inserted (with unit normal vector ˆ n = ˆ z ) to give Q = − πk B T η d c s d z Z a d rr ( a − r ) (cid:0) e − U ( r ) /k B T − (cid:1) . (S25)If ∆ c is the change in c s over the length of the pore, then using the fact that the flux Q isthe same for any cross-section along the pore gives Q = − πk B T η ∆ cL Z a d rr ( a − r ) (cid:0) e − U ( r ) /k B T − (cid:1) . (S26)S5or the step-function interaction potential (Eq. (17) of the main text, with d = r ) and aninteraction range smaller than the pore radius ( λ ≤ a ), from Eq. (S26) the solution flux is Q step = − πk B T η ∆ cL ( α − a − λ ) λ , λ ≤ a. (S27)On the other hand, for an interaction range much larger than the pore radius ( λ (cid:29) a ), e − U ( d ) /k B T can be approximated as a constant, which from Eq. (S26) gives the solution fluxas Q step ’ − πk B T η ∆ cL ( α − a , λ (cid:29) a. (S28)The solute flux density may be written, using Eq. (4) of the main text and Eq. (S20), as j = − De − U/k B T d c s d z ˆ z , (S29)which can be integrated over the pore cross-section using Eq. (14) of the main text to givethe total solute flux as J = − πD ∆ cL Z a d rre − U/k B T , (S30)using the fact that the flux J is conserved for any cross-section along the pore. For thestep-function interaction potential (Eq. (17) of the main text), the solute flux reduces to J step = − πD ∆ cL (cid:2) a + ( α − λ (2 a − λ ) (cid:3) , λ ≤ a (S31)and J step ’ − πD ∆ cL αa , λ (cid:29) a. (S32) S3. FINITE-ELEMENT NUMERICAL SIMULATIONS
The continuum hydrodynamic flow equations (Eqs. (1)–(3) of the main paper) were solvedusing finite-element method (FEM) simulations with COMSOL Multiphysics version 4.3a for a thin planar membrane containing a circular aperture of radius a connecting two largecylindrical fluid reservoirs (Fig. S2). The equations were solved using a fully coupled solver,which is a damped version of Newton’s method. The damping option used to achieveconvergence was “Automatic highly nonlinear (Newton)”. The PARDISO direct solver wasused.The solute–membrane interaction potential was modelled using Eq. (32) of the mainpaper. To avoid discontinuities in the potential, the membrane was given a finite thicknessS6IG. S2: Schematic (not to scale) of the two-dimensional axisymmetric computationaldomain used in the FEM calculations. Solid lines denote solid–liquid boundaries anddashed lines liquid boundaries. The geometry has rotational symmetry about theboundary AH. L and the membrane surface between points D and E in Fig. S2 was given a finite radiusof curvature of L/
2. The membrane thickness was less than 1/5 of the range of the solute–membrane potential in all cases and we verified that the measured solution and solute fluxeswere not sensitive to halving the membrane thickness from the value used for the resultspresented. A boundary layer mesh was used at all solid–liquid interfaces, with 5 boundarylayers, a thickness of the first layer of 1/2 the local domain element height and a boundarylayer stretching factor (mesh element growth rate) of 1.2. The maximum element size at theboundary between points D and E in Fig. S2 was min[ L/ , λ/
20] and that at the boundarybetween points C and D and between E and F was λ/
5, where λ was the solute–membraneinteraction range parameter. The maximum element size in the simulation domain was 28 σ ,where σ was the unit of length. We verified that the measured solution and solute fluxesdid not change significantly with a finer mesh. Table S1 lists the boundary conditions thatwere used to solve the equations and Table S2 lists the parameters used in the simulations.Several supplementary results from the FEM simulations are given in the figures below.Figure S3 shows the simulated solute concentration c and the ratio of the simulated toS7ABLE S1: Boundary conditions used to solve the continuum equations in the FEMsimulations. The vector ˆ n is the surface normal. boundary conditionsAH ˆ n · ∇ c = ˆ n · u = ˆ n · ∇ u = 0AB c = c H , p = p ∞ GH c = c L , p = p ∞ BC and FG ˆ n · j = ˆ n · u = ˆ n · ∇ u = 0CD, DE, and EF ˆ n · j = u = 0 TABLE S2: Parameters used for FEM calculations. The units of length, energy, and timeare σ , k B T , and τ = ησ / ( k B T ), respectively, where T is the temperature and η is thesolution shear viscosity ( σ can be considered the diameter of a fluid molecule). Where arange of values is given, the parameter was fixed at the value in parentheses unlessotherwise indicated. quantity symbol unit valuesolute diffusivity D σ /τ c ∞ σ − − membrane thickness L σ . a σ (cid:15) k B T -3–3 (1 /
10 or 1)solute–membrane interaction range λ σ c σ − − –3 × − reservoir radius/length w σ theoretical solute concentration, c/c theory , where c theory is the solute concentration derivedin the theory (given by Eqs. (5) and (8) of the main paper), near the pore entrance foran intermediate solute concentration difference between the reservoirs (similiar results wereobtained for the range of parameters simulated, with the maximum deviation between c and c theory of around ±
4% observed for (cid:15) = − k B T and ∆ c = 3 × − /σ ). These resultsindicate that the theory accurately predicts the nonequilibrium solute concentration in thisS8ystem geometry, even for a solute–membrane interaction strength at the membrane surfaceseveral times larger than k B T . Figure S4 gives illustrative examples of the simulated axialvelocity u z and axial solute current density j z as a function of the radial coordinate r atthe pore mouth at z = 0, showing in particular the effect of varying the solute–membraneinteraction strength (cid:15) . Figure S5 shows the simulated total solution flux and total soluteflux divided by the concentration difference between the reservoirs as a function of theconcentration difference, demonstrating that both fluxes varied linearly with the appliedconcentration difference for all of simulated parameters. (Results are not shown for varyingpore radius a or varying solute–membrane interaction range λ with the solute–membraneinteraction strength fixed at (cid:15) = k B T /
10, but they are qualitatively similar to those for (cid:15) = k B T .) We have verified that the diffusive solute flux (i.e. not including the contributionto the flux from solute advection) in all the simulations was essentially the same as the totalsolute flux, demonstrating that the simulations were all at low Pecl´et number. REFERENCES a) (b)(c) (d) FIG. S3: 2D plots of the ((a), (b)) solute concentration c (in units of σ − ) from FEMsimulations and ((c), (d)) c/c theory , where c theory is the solute concentration derived in thetheory (given by Eqs. (5) and (8) of the main paper), as a function of the r and z coordinates near the pore entrance for a solute–membrane interaction strength (cid:15) of((a),(c)) k B T /