Flow Control through Polymer-Grafted Smart Nanofluidic Channels

NANO LETTERS

Flow Control through Polymer-Grafted Smart Nanofluidic Channels: Molecular Dynamics Simulations

2005 Vol. 5, No. 12 2509-2514

Shashishekar P. Adiga*,† and Donald W. Brenner* Department of Materials Science and Engineering, North Carolina State UniVersity, Raleigh, North Carolina 27695-7907 Received September 14, 2005; Revised Manuscript Received November 7, 2005

ABSTRACT Presented are results of molecular dynamics simulations that demonstrate flow gating through a polymer-grafted nanopore as a function of effective solvent quality. Analysis of density and flow profiles from the simulations show that the difference in drag force exerted on the flowing solvent due to different polymer brush configurations produces the effective fluid gating. Shear-induced permeability changes through these nanopores has also been investigated. These results establish a critical starting point in nanofluidics from which continuum modeling can be developed to design this emerging class of smart nanoporous materials with tailor-made properties.

Polymer chains bonded to the inside of nanoporous structures form an emerging class of smart nanofluidic “valve” in which responsive gating is achieved through collapse transitions of the tethered polymers in response to changes in fluid conditions. It has been demonstrated, for example, that nanoand microporous membranes can be made with the ability to regulate fluid flow in response to changes in environmental variables such as pH, temperature, and solvent quality with respect to the grafted polymers.1-4 In principle, the properties of these smart materials can be tailored to respond to a variety of stimuli, making them useful for applications that include controlled drug delivery, ultrafiltration, molecular sorting, and smart chemical release. Grafting polymers can also reportedly improve fouling resistance of micro- and ultrafiltration membranes.5,6 Recently, molecular dynamics (MD) simulations were used to investigate drift of large molecules through a comb polymer grafted slit nanopore.7 It was demonstrated that solvent-responsive gating can be used to screen oligomers and that the slit nanopore is capable of size sorting. Here, the problem of interest is fluid flow through a polymer-grafted nanopore, in particular, stimuli-responsive control of Poiseuille-type flow and the interaction of flow with grafted polymers. There has been extensive research, both experimental and theoretical, into understanding the nature of fluid flow in nanofluidic geometries that are not modified with grafted polymers. These studies have demonstrated that continuum fluid mechanics equations are capable of accurately describ* Corresponding authors. E-mail: [email protected], [email protected]. † Present address: Materials Science Division, Argonne National Laboratory, 9700 S. Cass Ave., Argonne, IL 60439. 10.1021/nl051843x CCC: $30.25 Published on Web 11/22/2005

© 2005 American Chemical Society

ing Poiseuille flow through channels thicker than about 1020 atomic diameters.8 However, MD simulations have shown that for channels that are only a few atomic diameters wide, solid walls can induce ordering in liquids, which in turn can produce shear viscosities and microscopic stress tensors during flow that do not match predictions of the NavierStokes equations.9,10 Prior analyses of fluid flow in polymer grafted nanopores have used continuum hydrodynamic equations (e.g., the Brinkman equation11) that assume the grafted polymers form a brush layer with solvent permeability related to the monomer density.6,12 While useful for qualitative analysis, the accuracy with which flow characteristics can be obtained within a continuum approximation for these smart nanoscale pores has not been firmly established nor has an expression relating permeability to monomer structure needed for quantitative hydrodynamic calculations been derived. Prior atomistic simulations have focused on flow alongside a polymer layer under good solvent conditions13-16 and have not examined the system reported here. In addition, a majority of the atomistic studies have used an implicit solvent in which the solvent velocity profile inside the brush layer is obtained by solving the hydrodynamic Brinkman equation. Understanding the dynamics of flow through polymergrafted nanopores under both good and poor solvent conditions is a next level of complexity in nanofluidics because it involves the interplay of flow field and polymer brush structure in a confined geometry. As a first step toward a better understanding of fluid flow in these systems, we have carried out extensive MD simulations of flow through idealized polymer chains grafted to the inside of a cylindrical

nanopore under both good solvent conditions where the chains expand and inhibit flow and poor solvent conditions in which the chains collapse to open the pore. Polymer and solvent densities as well as flow fields as a function of distance from the center of the pore have been calculated, with the latter compared to flow fields for a nanopore without grafted polymers. For the fully extended chains, the fluid permeability, calculated from the simulations by integrating over the flow field, is reduced by a factor of 10 relative to the pore without the grafted polymers. At the fullest state of collapse, the grafted chains form multichain micelles, and the fluid permeability increases by a factor of 3 relative to the fully extended chains. The simulations clearly demonstrate that the difference in drag force exerted on the moving solvent by the extended and collapsed chains is sufficient to result in the magnitude of gating observed experimentally for this class of smart system. MD simulations of different flow rates for the case of a good solvent were also carried out. These studies demonstrate that when extended, grafted chains may be deformed along the flow direction at high shear rates, which results in shearinduced thinning of the grafted layer that increases the pore permeability. This shear-induced permeability change provides an additional mechanism for flow control through these nanopores, as has been experimentally demonstrated by Castro et al.6 The MD simulation system consists of polymer chains of length N end-attached to the inside of a cylindrical pore with radius R that is filled with a pool of explicit solvent atoms. The polymer chains, which are grafted at a density of one chain per an area of s, are represented by a bead-spring model. The grafting surface is represented by stationary wall atoms of the same size as the polymer beads placed on a cylindrical surface. Nonbonded interactions between monomer-monomer (mm), monomer-solvent (ms), solventsolvent (ss), monomer-wall (mw) atom, and solvent-wall (sw) atom pairs are modeled using a Lennard-Jones shifted force potential of the form USF(rij) )

{

ULJ(rij) - ULJ(rc) - (rij - rc)

(

)

dULJ(rij) drij

rij)rc

0

rij erc

(1)

rij > rc

where ULJ(rij) is the Lennard-Jones potential given by ULJ(rij) ) 4ij

{( ) ( ) } σij rij

12

-

σij rij

6

(2)

and rc is the cutoff distance beyond which the interaction is zero. The quantities ij and σij are energy and distance parameters, respectively, corresponding to a particular pair of particles i and j. Values of ij ) 0.5 kJ/mol and σij ) 2.7 Å are used for all nonbonded pair interactions unless stated otherwise. A mass of 14 amu is assigned to each particle. Neighboring monomers along a chain interact via a harmonic bond potential with the value of the spring constant kb ) 2510

13.72 kJ/(mol Å2) and equilibrium distance ro ) σmm. The cylindrical nanopore, which has radius of R ) 30σmm and a length of L ) 120σmm, is grafted with polymer chains with N ) 100, s ) 100σmm2. The pore is filled with 142 191 solvent atoms of the same size and mass as the polymer beads. The total number of particles in the system is 187 753 at a reduced number density of 0.5533. A periodic boundary condition is applied along the z-axis, which is parallel to the axis of the cylinder, and a single-body external force Fz on each solvent atom generates flow in the z-direction. A constant temperature of T ) 200 K is maintained by coupling the system to a Nose´-Hoover thermostat. The MD simulations, which used a time step of 5 fs, were performed using the DL_POLY_2.12 simulation code.17 The structure was equilibrated for 2 × 105 steps before the external field was switched on. Once the external field was switched on, the simulations were carried out for a sufficiently long time such that the flow was fully developed and the grafted polymer layer responded to the flow. Typically, sampling was begun 5 × 105 time steps after switching on the flow field. The solvent velocity profile was generated by calculating the time average of the z component of the velocity of the solvent atoms that are at a distance r from the center of the pore, averaged over 5 × 105 time steps. To analyze the system’s ability to control flow, MD simulations of flow under five different solvent conditions were carried out. For the good solvent case the truncation distance rc in eq 2 for the mm interactions was chosen as 21/6σmm, which truncates the potential such that only repulsive interactions are included (system A). For the other four solvent conditions, a cutoff of 2.5σmm (which includes the attractive part of the potential) and an interaction potential mm ) 0.2, 0.3, 0.5, and 1.0 kJ/mol (systems B, C, D, and E, respectively) were used for the mm pair. In all five systems sw, mw, ss, and ms interactions were truncated so that they are purely repulsive. This choice of parameters was made so that the mm interaction was increasingly favored over the ms interaction and the effective solvent quality becomes progressively poorer as the system goes from A to E. The standard way of varying the solvent quality is by either changing the temperature or changing the well depth of the van der Waals interaction between explicit solvent particles and polymer beads. However, for the current study, badness in solvent quality is induced by making the polymerpolymer interaction progressively attractive. The rationale behind this is as follows. The purpose of this study is to investigate drag on flowing particles as a function of conformation of the grafted chains alone. Varying either the temperature of the system or the solvent-polymer interaction to induce conformational changes would also imply changing viscosity of the flowing solvent. There are experimental systems where a change in interaction between monomers alone induces conformational transition. For example, Ito et al. have demonstrated pH-responsive control of water flow through a polypeptide-grafted nanoporous membrane.1 Essentially the Coulombic repulsion between monomers upon Nano Lett., Vol. 5, No. 12, 2005

Figure 1. Snapshots from the MD simulations of systems A (left) and E. Fluid atoms are not shown for better visualization.

deprotonantion at high pH drives the helix-coil transition in grafted chains. Snapshots from MD simulations for systems A and E are illustrated in Figure 1. In the case of system A, which corresponds to a good solvent, the chains are in an extended conformation and are stretched away from the surface. As the attractive part of the mm interaction is increased, the chains in systems B and C start to collapse. The further increase in mm interaction results in the grafted chains aggregating to form pinned micelles with globular cores and extended legs connecting the core with the grafting point (Figure 1). The formation of pinned micelles has been proposed by Williams18 and Zhulina et al.19 in the case of sufficiently dense brushes exposed to a poor solvent. In system E, which had the largest degree of collapse, there were on average nine chains per micelle. Plotted in the top panel of Figure 2 is the monomer density profile for systems A, B, and E. The monomer density profile is diffuse and extends to the center of the cylindrical pore for the system A, and drops off within a short distance of the pore wall in the case of system E, leaving a central polymer-free region in the pore. The density profile for system B is between those of A and E. The solvent density profile plot (center panel Figure 2) illustrates that in each case solvent atoms are slightly expelled from the polymer layer, with the degree of expulsion from the sides of the pore increasing with decreasing solvent quality. However, there remains a significant amount of solvent throughout the pore in each case. Plotted in the bottom panel of Figure 2 is the average velocity of the solvent atoms as a function of distance from the pore center for an applied force Fz ) 1fe for each of the three systems where fe ) 1.6605 × 10-3 pN. For comparison the flow profile in a polymer-free pore is also plotted, where a parabolic profile characteristic of circular Poiseuille flow is obtained. In contrast to the solvent density profiles, where the density was relatively unchanged throughout the pore, the structure of the polymer due to the solvent quality significantly influences the velocity profile. For the good solvent case (A) the extended chains produce significant drag on the solvent atoms throughout the pore, resulting in a flow Nano Lett., Vol. 5, No. 12, 2005

Figure 2. Density and flow profiles from the simulations. Monomer density (top), solvent density (center), and solvent velocity (bottom) profiles for systems A (4), B (0), and E (2). Also plotted in bottom graph is the velocity profile inside a polymer free pore ([). Fz ) 1fe.

velocity that is greatly reduced across the pore, with the only significant flow near the pore center. As the solvent quality decreases, and hence the chains begin to collapse, the velocity of the solvent atoms in the pore center increases and the distance from the pore wall at which solvent moves decreases. Both of these effects result in a higher overall 2511

flow rate and hence a larger effective pore size. Hence the changes in permeability are not due to a literal change in pore size reflected in the solvent density, but rather for this system the change in pore size is due to the change in velocity profile resulting from changes in the drag forces acting on the solvent atoms from the different degrees of collapse of the polymer chains. To further quantify the simulated smart nanopore’s ability to regulate flow in response to effective solvent quality, a dimensionless fluid permeability K of the pore is defined as K ) Q/Q0. Here Q and Q0 are flow rates through the pore with and without a grafted polymer, respectively, under the same flow conditions. From standard fluid mechanics equations, the flow rate through the polymer-free nanopore is given byQ0 ) (πR4FzF)/8µ, where F and µ are the density and viscosity of the solvent, respectively. The viscosity of the Lennard-Jones fluid used in the MD simulations was determined by carrying out MD simulations of Poiseuille flow of the fluid through a polymer-free cylindrical nanopore at Fz ) 1fe. The velocity profile generated was fit to a parabolic profile, and the viscosity was estimated from the equation for Q0. The value of viscosity thus extracted is used in calculating Q0 at different applied flow fields. Flow rates through the polymer-grafted nanopore are extracted from the MD velocity profiles as Q)

∫0R 2πrVz(r) dr

(3)

where νz(r) are the flow distributions plotted in Figure 2. To characterize the effect of solvent quality on layer thickness, we define the height of the grafted polymer layer as the square root of the second moment of monomer density (Fm) distribution

H)

(

)

∫0R (R - r)2Fm dr 1/2 ∫0R Fm dr

(4)

Plotted in Figure 3 are K (top) and H/HA (bottom) as a function of mm. HA is the thickness of the grafted layer in a good solvent (system A). Overall the flow rate changes by a factor of 10 between the polymer-free pore and the fully extended polymer and by a factor of about 3 between the fully extended and fully collapsed states. Changes in the fluid permeability also correlate well with changes in the normalized height of the brush structure. When the chains are collapsed, there is more room in the pore within which the solvent can flow without significant drag, while the extended structure produces increased drag throughout the pore. Although an idealized polymer model, our results are consistent with experimental studies. Park, Toshihiro, and Imanishi characterized fluid permeability as a function of pH for water flowing through a nanoporous material into which a self-assembled thiol brush of poly(glutamic acid) was deposited.1 In their system, pore radii ranged from about 56 to 160 nm, with a majority of pores having radii of about 80 nm. The water permeability decreased by a factor of about 2512

Figure 3. Fluid permeability (top) and normalized grafted layer height H/HA (bottom) as a function of mm. HA is the grafted layer height in system A. The broken line in the top panel corresponds to the permeability of system A. Fz ) 1 fe.

2.3 as the water went from an acidic (pH of 3) to a neutral pH due to a configuration change within the polypeptide brush. In related experiments Mika, Childs, and Dickson3 demonstrated water gating with changes in pH using poly(4-vinylpyridine) anchored within a microporous membrane. Pores with a diameter of 820 nm showed relative permeabilities of about 1600, while smaller pores of diameter 192 nm showed relative water permeabilities of about 85. Extrapolating their relative permeabilities versus pore size data using an exponential function yields a relative permeability of 2.7 for a pore diameter of 16 nm. While the relative permeabilities will depend strongly on properties such as the molecular weight and grafting density of the polymer, the pore size distribution, and the change in configuration of the polymer when extended and collapsed, the favorable comparison between our simulations and the experimental studies suggests that the idealized grafted polymer brush captures the essential physics of this system. Reported in the rest of this Letter are results from an investigation of the effect of solvent flow rate on permeability for a polymer-grafted nanopore under good solvent conditions. The shear response of grafted polymers under good solvent conditions has been a subject of considerable debate.20,21 Using a surface force apparatus, Klein et al. demonstrated that a pair of opposing brushes slide past each other, in a good solvent, experiencing a repulsive normal force.22 This repulsion between sheared brushes was interpreted as shear-induced swelling of brushes. Following these experiments, theoretical predictions based on polymer scaling arguments23 and hydrodynamic models24 have mostly supported shear-induced swelling. On the basis of scaling arguments, Sevick and Williams proposed a novel design for a “smart” constant discharge flow-control valve made Nano Lett., Vol. 5, No. 12, 2005

from two opposing polymer brushes.25 The idea behind the proposed design is that as the flow velocity through the polymer grafted valve increases, the brush layer swells, in effect causing a reduction in overall flow and thus maintaining a constant flow rate over a range of pressure gradients. These models have simplified the solvent flow with a shear force applied to the free ends of the grafted chains, essentially ignoring the details of solvent flow inside the brush layer. In addition, these models consider a step-function monomer density profile in the grafted layer as opposed to a parabolic density profile predicted by self-consistent field theories.26 More accurate modeling of the shear response of polymer brushes has been performed using Monte Carlo,13 Brownian dynamics,14,15 and MD simulations.16 A majority of these studies have predicted little or no increase in brush thickness in contrast to the analytic models. In fact, they have predicted a reduction in brush height at high shear rates owing to deformation of grafted chains in the direction of flow. In subsequent experimental work, Ivkov et al. performed X-ray and neutron reflectivity studies on the effects of solvent flow on the height and density profile of a polymer brush and observed no change in either quantity for a range of shear rates.20 While the simulations have eliminated the need for assumptions regarding monomer density profile and chain deformation in response to a flow field, they still use a continuum model for the solvent. These simulations treat the brush layer as a porous medium and solve the Brinkman hydrodynamic equation to obtain the solvent velocity profile inside the brush layer. Also, these investigations have focused on a pair of opposing polymer-brush-bearing surfaces sliding past each other. Our simulations, on the other hand, involve flow of explicit solvent atoms and are aimed at investigating changes in permeability due to flow through a polymergrafted nanopore. More recently, Castro et al. have demonstrated an increase in hydraulic permeability of a polymer-grafted microporous silica membrane at high flow rates.6 They attribute this shearinduced increase in permeability to a reduction in polymer layer thickness (up to 47%) caused by significant chain deformation. Shear-induced thinning is expected in the case of polymer chains in coil conformations when subjected to a shear flow that is stronger than the entropic restoring force. When the rate of deformation (applied strain rate dγ/dt) is faster than the time required for the chain to relax to the equilibrium configuration (characterized by relaxation time of the polymer chain, τR), one expects stretching of the chain in the flow direction. The strength of the flow field is characterized by the dimensionless Weissenberg number Wi ) τR dγ/dt.27 To determine the effect of flow rate on the structure of a grafted layer and in turn the permeability of the nanopore, a series of simulations of system A (good solvent) was conducted under flow fields of Fz ) 1, 2, 3, 4, 5, 10, 15, and 20 fe. Depicted in Figure 4 are monomer density profiles for Fz ) 1, 10, 15, and 20 fe. For comparison the monomer density profile for Fz ) 0 (no flow field) is also shown. It is clear from the density profiles that flow alters the brush structure, and the higher the shear rate the more compact is Nano Lett., Vol. 5, No. 12, 2005

Figure 4. Monomer density profile inside the pore for Fz ) 0 (×), 1 (0), 5 (4), 10 (9), and 20 (O) fe.

Figure 5. Hydraulic permeability (top) and normalized grafted layer height H/H0 (bottom) as a function of Weissenberg number. H0 is the grafted layer height at Fz ) 0. The points in the plot correspond to values of flow field Fz ) 1, 2, 3, 4, 5, 10, 15, and 20 fe.

the density profile. The density profiles for Fz ) 0 and 1 (which was used in the simulations reported above) are effectively the same indicating that in the latter case flow is not strong enough to cause significant chain deformation. To characterize the effect of shear on hydraulic permeability of the pore the permeability K was determined as a function of Wi. The flow rates required for calculating K are obtained using velocity profiles extracted from the simulations (eq 3). A value for Wi for each value of flow field is calculated as follows. The strain rate is taken to be dγ/dt ) Vz(r)0)/R. To determine the polymer relaxation time, an MD simulation of an isolated free polymer chain of length 100 in a pool of solvent atoms was carried out. A value for τR was calculated to be approximately 125 ps by fitting to a single exponential the decay of the autocorrelation function of the deviation of the radius of gyration from the mean value. Plotted in Figure 5 are permeability K and normalized brush thickness H/H0 as a function of Wi, where H0 is the brush height under no flow. There is an increase in K of about 2.5 times in going from Wi ) 0.07 to 3.04. The corresponding decrease in brush height is about 34%. In conclusion, reported in this Letter are the results of an investigation of simulated solvent flow through polymergrafted cylindrical nanopores. Using MD simulations it was demonstrated that the system is capable of changing the permeability drastically upon changing the solvent quality 2513

from good to poor. It was also established that for a good solvent shear can deform the brush layer altering the permeability of the system and thus providing yet another mechanism of flow control. The onset of shear-induced permeability is related to the dimensionless Weissenburg number that characterizes the flow strength with respect to polymer relaxation time. One could imagine a smart nanovalve that releases solvent in a container by adjusting its permeability in response to pressure inside. The MD simulation results presented here establish a critical starting point from which continuum mechanics based modeling can be developed to design smart nanoporous materials with tailor-made properties. The need for continuum modeling arises as MD simulations are computationally very expensive to use in designing specific polymersolvent systems. We will report in a more detailed publication a continuum modeling approach to designing polymer-grafted flow control valves. The simulation results presented here form an independent means to evaluate the validity of continuum models. Acknowledgment. This work was supported by NASAAmes through a NASA-North Carolina State University cooperative research agreement. The MD calculations were performed through a grant of computing time from the North Carolina Supercomputing Center. D. Srivastava and M. Meyyapan are thanked for their support of this work. References (1) Park, Y. S.; Toshihiro, I.; Imanishi, Y. Langmuir 2000, 16, 5376. (2) Akerman S.; et al. Int. J. Pharm. 1998, 164, 29. (3) Mika, A. M.; Childs, R. F.; Dickson, J. M. J. Membr. Sci. 1999, 153, 45.

2514

(4) Iwata, H.; Hirata, I.; Ikada, Y. Macromolecules 1998, 31, 3671. (5) Castro, R. P.; Cohen, Y.; Monbouquette, H. G. J. Membr. Sci. 1996, 115, 179. (6) Castro, R. P.; Monbouquette H. G.; Cohen, Y. J. Membr. Sci. 2000, 179, 207. (7) Adiga, S. P.; Brenner, D. W. Nano Lett. 2002, 2, 567. (8) Koplik, J.; Banavar, J. R.; Willemsen, J. F. Phys. Fluids A 1989, 1, 781. (9) Travis, K. P.; Todd, B. D.; Evans, D. J. Phys. ReV. E 1997, 55, 4288. (10) Travis, K. P.; Gubbins, K. E. J. Chem. Phys. 2000, 112, 1984. (11) Brinkman, H. C. Appl. Sci. Res. 1947, A1, 27. (12) Milner, S. T. Macromolecules 1991, 24, 3704. (13) Lai, P. Y.; Binder, K. J. Chem. Phys. 1993, 98, 2366. (14) Doyle, P. S.; Shaqfeh, E. S. G.; Gast, A. P. Phys. ReV. Lett. 1997, 78, 1182. (15) Neelov, I. M.; Borisov, O. V.; Binder, K. Macromol. Theory Simul. 1998, 7, 141. (16) Grest, G. S. Phys. ReV. Lett. 1996, 76, 4979. (17) DL_POLY is a package of molecular simulation routine written by W. Smith and T. R. Forester, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington, 1996. For details see http://www.dl.ac.uk/ TCSC/Software/DL_POLY/main.html (18) Williams, D. R. M. J. Phys. II 1993, 3, 1313. (19) Zhulina E. B.; et al. Macromolecules 1995, 28, 8612. (20) Ivkov, R.; Butler, P. D.; Satija, S. K.; Fetters, L. J. Langmuir 2001, 17, 2999. (21) Doyle, P. S., Shaqfeh, E. S. G.; Gast A. P. Macromolecules 1998, 31, 5474. (22) Klien, J.; Perahia, D.; Warburg, S. Nature 1991, 353, 143. (23) Barrat, J.-L. Macromolecules 1992, 25, 832. (24) Kumaran, V. Macromolecules 1993, 26, 2464. (25) Sevick, E. M.; and Williams, D. R. M. Macromolecules 1994, 27, 5285. (26) Milner, S. T.; Witten, T. A.; Cates, M. E. Macromolecules 1998, 21, 2610. (27) Doyle, P. S.; Ladoux, B.; Viovy, J. Phys. ReV. Lett. 2000, 84, 4769.

NL051843X

Nano Lett., Vol. 5, No. 12, 2005