NANO LETTERS
Dynamics of the Capillary Rise in Nanocylinders
2002 Vol. 2, No. 11 1281-1285
Lev D. Gelb* and Alicia C. Hopkins Department of Chemistry, Florida State UniVersity, Tallahassee, Florida 32306-4390 Received August 13, 2002; Revised Manuscript Received September 5, 2002
ABSTRACT Molecular dynamics simulations are used to investigate the dynamics of fluid flow into empty cylindrical pores with diameters between 1.5 and 5 nm. Model system parameters are chosen to represent xenon in silica pores. The principal data collected are the amount of fluid in the cylinders and the density profile of the fluid. The effects of both pore diameter and fluid temperature are considered. The temperature is found to have little effect on the penetration of the fluid into the pore until temperatures close to the critical point are reached. The cylinder diameter has a substantial effect, and for diameters g2 nm, the rate of fluid flow into the pore is proportional to its cross-sectional area. Last, at subcritical temperatures in all but the smallest (1.5 nm) cylinders, the total amount of fluid in the cylinder is observed to grow as the square root of the time.
1. Introduction. Fluid flow in porous materials and capillaries has been studied extensively and is relevant in many industrial and laboratory processes.1-5 With recent advances in nanotechnology, it is possible to observe fluid flow in single nanoscale capillaries, especially various types of nanotubes.6,7 Hence the applicability of conventional fluid dynamics in nanoscale systems, as well as theories of the capillary rise, may be tested directly. Such nanoscale capillaries are important components in future “lab-on-achip” designs,8-10 in nanofluidic devices,6-10 and for the delivery of quantities of fluid as small as a few hundred molecules.11 We have performed a series of molecular simulations of the uptake of a simple fluid into capillaries a few nanometers wide. We have identified and visualized the flow behavior that may be observed, estimated the rate of fluid uptake, and tried to better understand the effects of different system parameters on this rate. Molecular simulation is an ideal tool for such exploratory studies because the length-scales and time-scales of interest can be easily accessed, system geometries are precisely controlled, and microscopic quantities such as the number of atoms in the pore can be measured directly. 2. Model and Simulation Details. A simplified model for xenon adsorption on silica was used throughout this study. The pore material was modeled with an amorphous configuration of particles representing the oxygen atoms in silica. Silicon atoms are omitted, following other studies.12,13 Silicon atoms are small, of low polarizability, and are generally located below the surface; therefore, they would have only * Corresponding author. Present address: Department of Chemistry, Washington University, Campus Box 1134, One Brookings Drive, St. Louis, MO 63130. 10.1021/nl025748p CCC: $22.00 Published on Web 09/28/2002
© 2002 American Chemical Society
very weak interactions with a nonpolar adsorbing fluid, and their omission does not roughen the surface. The material model was generated by simulating a Lennard-Jones liquid at the density corresponding to that of oxygen atoms in silica glass, and then filling the few remaining holes in this configuration with additional Lennard-Jones particles. To prepare pore models, cylindrical holes were then “drilled” out of the material, as explained below. No surface relaxation was applied, and pore-wall atoms were kept frozen throughout the fluid uptake simulations. The pore models used in this study have surfaces that are rough on an atomic scale, although quite smooth on larger scales. Such models have also been used in previous studies.14,15 The xenon fluid was also modeled with a Lennard-Jones potential. The potential parameters used for xenon were σ ) 0.391 nm, /kB ) 267.5 K, and for the pore particles, σ ) 0.27 nm, and /kB ) 271 K. The potential interactions between all particles were truncated at 2.5σXe. Because in molecular dynamics simulations only forces contribute to the behavior of the system, this corresponds to a “cut-andshifted” potential for a Monte Carlo simulation. These values were chosen so that the critical temperature of the model fluid (Tc ) 1.079 in reduced units, with this cutoff16) would correctly describe xenon, with a critical temperature of 289 K. The value of pore was chosen to maintain the same Xe/pore ratio as in earlier work,17 and the values of σ are those used in the same studies. Interspecies parameters were generated from the Lorentz-Berthelot mixing rules. In ref 17, the values of Xe and pore used made for poor quantitative comparison between simulated and experimental isotherms, as the model fluid had different critical properties than the experimental one. Good qualitative agreement was achieved, however, and the parameters given above, which do have
Figure 1. Schematic of system geometry and preparation. A liquid in equilibrium with its own vapor is placed in contact with a completely evacuated cylindrical pore mouth, via a three-step process, and is subsequently drawn into the pore.
Figure 2. Density profiles of fluid uptake in a 3 nm cylinder at 217 K. The four profiles shown were taken at times of 2, 10, 18, and 26 ns, in descending order. There is considerable noise in the density histogram in the center of the pore, due to the very small volumes being sampled in the radial histogram. The pore material is not shown.
the correct critical properties, should be a substantial improvement. The system geometry is shown in Figure 1. The simulation cell is periodic in all three directions. The configuration at the start of the fluid uptake simulations is prepared as follows. Initially, the cell contains only a block of pore material (Figure 1a). A quantity of liquid is excised from a separate simulation cell, placed in contact with this material, and molecular dynamics are run until the total energy stabilizes (Figure 1b). During this process, some of the liquid evaporates to fill the void space with vapor, and some condensation of liquid occurs on the other side of the block of solid material. This adsorbed layer is present in the rest of the simulation but does not have any effect on the porefilling process as it is too far away. Finally, a cylindrical channel (Figure 1c) is excised from the pore material in a single time step, without any manipulation of the fluid atoms. The clock is then reset to t ) 0; as the molecular dynamics simulation continues from this point, fluid will flow into the cylinder, drawn by the affinity of the fluid atoms for the pore surface. This geometry can be considered as a first approximation to an experimental system where a material containing cylindrical nanopores is brought into contact with a liquid surface. Direct simulation of such a procedure is also possible, but introduces the velocity profile of the moving surface as an additional complication. A similar approach was previously used by Spohr et al.18 to study the adsorption of water into attractive and repulsive slit pores, though with a focus on equilibrium properties rather than dynamics. This process was repeated for each fluid temperature and pore diameter. Pore diameters of 1.5, 2.0, 3.0, 4.0, and 5.0 nm were considered, and fluid temperatures of T ) 190, 217, 244, and 271 K were run for each diameter. The total number of fluid particles in the system ranged between 10 000 and 20 000; more fluid was used for the deepest and widest pore systems. The x and y dimensions of the periodic cell were both 8.1 nm. In almost all cases, the depth of the cylindrical
channel was 31 nm, with a few simulations done in 47 nm deep cylinders for comparison. A “Gaussian constraint” constant-temperature thermostat was used throughout these simulations, implemented within a fifth-order Gear predictor-corrector integration algorithm.19 The time-step is 10.3 fs (0.005 reduced time units.) Thermostatting is required to remove the heat generated as the fluid contacts the pore interior, or the temperature would rise very quickly in the simulation. Given that the motion of the fluid is quite slow on a molecular time scale and that, in an experimental system, the porous material would serve as an effective heat bath, the constant-temperature approach seems reasonable. Simulations were typically run for a total of 4 million time steps, or just over 40 ns. In all cases the liquid layer outside the pore remained thick, ensuring that the fluid in the pore was in contact with a bulk system at liquid-vapor equilibrium; for the single-component fluid used here, this condition and the temperature completely specify the thermodynamic state point. 3. Results and Discussion. For visualization purposes, we have collected two-dimensional histograms of the fluid density in the radial and axial directions over 30 ps time slices, which can be viewed as a series of images. Figure 2 shows a series of these density profiles showing the progress of fluid into the 3 nm cylinder at 217 K. In the first image, a well-defined hemispherical meniscus is clearly visible in the cylinder. Fluid coats the first 40% of the cylinder depth, but there is essentially no density further inside. The liquid is seen to be strongly layered near the surface, as has been observed in many other simulation studies.15,20,21 As time progresses, the hemispherical meniscus moves smoothly into the cylinder, and the cylinder walls become coated with fluid at greater depth. In the third image, only a small volume of vapor remains in the cylinder, and in the last image it is completely filled. Outside the cylinder, the thickness of the liquid layer is reduced. In these simulations the rate of motion of the liquid meniscus is approximately 1 nm/ns. A second sequence of density snapshots, taken in the 4 nm diameter
1282
Nano Lett., Vol. 2, No. 11, 2002
Figure 3. Density profiles of fluid uptake in a 4 nm cylinder at 271 K. The four profiles shown were taken at times of 2, 10, 14, and 22 ns, in descending order. While the initial stages are similar to those shown in Figure 2, the later behavior is quite different, with large density fluctuations visible in the cylinder.
Figure 5. Number of Xe atoms in the cylinder vs time, n(t), for all systems studied. Graphs (a-e) each show data for a particular pore diameter at four different temperatures. The legend in (a) applies to graphs (a-e). In (d), the two dot-dashed lines are data taken in the extended-length cylinder at 190 and 217 K. Graph (f) shows data for different cylinder diameters at 190 K, for comparison.
Figure 4. Density profiles of fluid uptake in a 1.5 nm cylinder at 217 K. The four profiles shown were taken at times of 2.6, 7.7, 15.5, and 23.2 ns, in descending order. The fluid is seen to be very highly structured in the pore, and no well-defined meniscus is visible.
system at 271 K, is shown in Figure 3, and a third series, taken in the 1.5 nm system at 217 K, is shown in Figure 4. Because the accumulation time for these density histograms is very short (only 3000 time steps), some artifacts are observed. The first is that a noisy white line is present along the axis of the cylinder. This occurs because the volume elements in the (r, z) histogram are extremely small near the axis, and therefore are often empty for the entire time of accumulation. The second artifact is that the bulk fluid outside the pore appears to show density variations; because of the very short time of averaging, atoms do not move very far from their original position, leading to the illusion of inhomogeneity even far from the material surface. Both of these effects would be alleviated by taking longer time averages, but then the liquid-vapor menisci in the cylinder would become blurred due to their motion during the time of accumulation, obscuring important detail. In Figure 5 are shown the principal data collected in this study, the number of particles in the cylinder as a function of time, n(t), for all systems at all temperatures. Each set of Nano Lett., Vol. 2, No. 11, 2002
curves is associated with a single pore diameter; the temperature of the fluid has only a minor influence on the rate at which fluid enters the pores. In all cases, the n(t) curves flatten out after approximately 25 ns, because the fluid has contacted the back end of the cylinder. This is visible in the last density profile in Figure 2, for example. Two additional runs were made using a 4.0 nm pore system that was more than twice as long; these data are shown in Figure 5d. In the longer system, no flattening out of n(t) was observed, and the divergence of data from the long systems and short ones occurs precisely at the time when the fluid visibly contacts the back end of the cylinder. In the discussion that follows, only data taken before this contact time are considered, because only these data are free of effects introduced by the use of a finite-depth cylinder. The temperature dependence of the data is as follows. For the 3.0, 4.0, and 5.0 nm diameter cylinders (Figure 5c, d, and e), data taken at all temperatures superimpose well except for those at the highest temperature, 271 K. For the 1.5 and 2.0 nm cylinders (Figure 5a and b), the data superimpose reasonably well except for those at the lowest temperatures, 190 and 217 K. That is, for larger cylinders, the flow rate is independent of temperature at low temperatures, and for smaller cylinders, the flow rate is independent of temperature at higher temperatures. In the three larger cylinders, the deviation of 271 K data from that taken at lower temperatures becomes more pronounced as the diameter is increased. This 1283
temperature is quite close to the critical temperature, and we observe in density histograms that the fluid in the pore exhibits substantial density fluctuations in the axial direction and that the liquid-vapor meniscus is much broader. An additional pore-filling mechanism seems to be at work in these systems, as fluid that has moved along the walls can collect in the center of the cylinder in adVance of the hemispherical meniscus, forming small bubbles which are later filled in. Such an event is visible in the end of the sequence shown in Figure 3. One would expect that the amount of fluid contained in a cylinder under pressure-driven flow conditions should be proportional to the cross-sectional area of the cylinder. If this is the case here, then dividing each data set by its crosssectional area should “collapse” all the data onto the same curve. As has previously been observed when discussing nanoporous materials, the area in question will depend on the choice of molecular surface used to define the diameter of the system and is to some extent arbitrary.14 The radius used in constructing the pore systems is not an appropriate choice, because it defines the distance from the center of the cylinder to the centers of the atoms on the pore surface; the actual volume accessible to the fluid is smaller than this. Instead, we have chosen to use an effective radius reff ) r - σfw/2; this takes into account the volume excluded by the radius of the pore atoms, σfw/2, and corresponds to the molecular surface definition commonly used in simulations of biochemical systems.14,22,23 By scaling n(t) data by effective diameter, we found that in the 1.5 nm cylinder the amount of fluid is not proportional to the cross-sectional area. As mentioned above, the 1.5 nm cylinder exhibits different behavior at 190 and 217 K than at the higher temperatures. Density plots for fluid flow in this cylinder at 217 K are shown in Figure 4. The cylinder is so narrow that exactly four liquid layers fill it entirely, and a hemispherical meniscus such as those seen in larger cylinders is not visible. The flow in this system appears to be principally along the cylinder walls. The fluid is highly structured in this environment, and we believe that this substantially retards the rate of fluid penetration into the pore. This effect would be accentuated by low temperature, because at lower temperatures simple liquids such as xenon are denser and more highly structured. A similar deviation from temperature independence is observed in the 2 nm pore at 190 K. We have attempted to collapse a large subset of the scaled n(t) data onto a single curve. This is only possible within the temperature limits described above, so that large-pore data at high temperatures and small-pore data at low temperatures are not used. Data from the 1.5 nm system are not used, as different flow mechanisms are clearly present in the smallest cylinder. Every remaining data set is divided by the appropriate effective area as calculated from reff above. Finally, the combined data were fit to a simple power-law form n(t)/Aeff ) BtC, with B and C as the fit parameters. The results are shown in Figure 6. The superposition of data is excellent, especially given that each simulation was run only once and that the n(t) curves are not perfectly smooth. The 1284
Figure 6. Number of particles in the cylinder as a function of time, n(t), data for all systems excluding 3, 4, and 5 nm cylinders at 271 K, all 1.5 nm systems, and the 2.0 nm system at 190 K, scaled by effective area. The solid curve is the accompanying power-law fit.
fit to the exponent gave C ) 0.497 ( 0.005, with R2 ) 0.992. This can be interpreted with the following simple model. The principal thermodynamic driving force for fluid flow into the cylinder is the heat of adsorption of the fluid against the cylinder wall, which is proportional to the distance that the fluid front has penetrated, l(t). A reasonable assertion is that the work required to move the fluid front one unit distance should be proportional to the amount of fluid that needs to be moved, itself proportional to l(t). Then, one obtains dl(t)/dt ∝ 1/l(t), which integrates to l2(t) ∝ t, or equivalently, l(t) ∝ xt. This is exactly what is observed; a power-law growth in the quantity of fluid in the pore, with an exponent of 1/2, to within numerical error. This model assumes that the motion of the fluid front can be simply described with a scalar function l(t). In fact, the density visualizations shown earlier suggest that the fluid may move along the pore walls and through the pore center at different rates, which may indicate that this simple picture applies only during the initial stages of the pore filling process, and that different mechanisms and rate laws may appear at later times. The rates of fluid flow are such that, for an experimental observation of systems comparable to those described here, nanosecond time resolution would be required. In longer cylinders lower time resolution would be required, provided that this flow model holds at longer times. If such capillaries were to be used for the delivery of small amounts of fluid, the fluid would have to be removed from the capillary, which could be accomplished by heating. We have seen that substantial deviations from simple continuum-like fluid behavior occur only in the smallest (1.5 nm) cylinder, in which the cylinder diameter is only about four times the particle diameter. While it is not at all surprising that fluid in such a small cylinder has different flow properties than in the bulk, it is curious that the slightly larger 2.0 nm cylinder (five particle diameters) does follow the effective-area proportionality described above, because Nano Lett., Vol. 2, No. 11, 2002
the fluid is highly structured in that system as well. To better understand this, closer analysis of the fluid structure and flow will be required in additional studies. In the future, this work will also be extended to carbon nanotubes, which have substantially smoother (and more strongly adsorbing) walls than the cylinders used here and may show more dramatic effects due to enhanced liquid structure. References (1) Ford, D. M.; Heffelfinger, G. S. Massively parallel dual control volume grand canonical molecular dynamics with LADERA II. Gradient driven diffusion through polymers. Mol. Phys. 1998, 94(4), 673-683. (2) Kainourgiakis, M. E.; Kikkinides, E. S.; Stubos, A. K.; Kanellopoulos, N. K. Simulation of self-diffusion of point-like and finite-size tracers in stochastically reconstructed Vycor porous glasses. J. Chem. Phys. 1999, 111(6), 2735-2743. (3) MacElroy, J. M. D. Nonequilibrium molecular dynamics simulation of diffusion and flow in thin microporous membranes. J. Chem. Phys. 1994, 101, 5274. (4) Tamai, Y.; Tanaka, H.; Nakanishi, K. Molecular simulation of permeation of small penetrants through membranes. 1. Diffusion coefficients. Macromolecules 1994, 27, 4498-4508. (5) Travis, K. P.; Gubbins, K. E. Combined diffusive and viscous transport of methane in a carbon slit pore. Mol. Sim. 2000, 25(3-4), 209-227. (6) Gogotsi, Y. Libera, J. A.; Guvenc-Yazicioglu, A.; Megaridis, C. M. In situ multiphase fluid experiments in hydrothermal carbon nanotubes. Appl. Phys. Lett. 2001, 79(7), 1021-1023. (7) Megaridis, C. M.; Yazicioglu, A. G.; Libera, J. A.; Gogotsi, Y. Attoliter fluid experiments in individual close-end carbon nanotubes: liquid film and fluid interface dynamics. Phys. Fluids 2002, 14(2), L5-L8. (8) Karlsson, R.; Karlsson, M.; Karlsson, A.; Cans, A. S.; Bergenholtz, J.; Akerman, B.; Ewing, A. G.; Voinova, M.; Orwar, O. Movingwall-driven flows in nanofluidic systems. Langmuir 2002, 18(11), 4186-4190.
Nano Lett., Vol. 2, No. 11, 2002
(9) Pepin, A.; Youinou, P.; Studer, V.; Lebib, A.; Chen, Y. Nanoimprint lithography for the fabrication of DNA electrophoresis chips. Microelec. Eng. 2002, 61-2(927-932). (10) Cao, H.; Yu, Z. N.; Wang, J.; Tegenfeldt, J. O.; Austin, R. H.; Chen, E.; Wu, W.; Chou, S. Y. Fabrication of 10 nm enclosed nanofluidic channels. Appl. Phys. Lett. 2002, 81(1), 174-176. (11) Piner, R. D.; Zhu, J.; Xu, F.; Hong, S.; Mirkin, C. A. “Dip-pen” nanolithography. Science 1999, 283(5402), 661-663. (12) Brodka, A.; Zerda, T. W. Molecular Dynamics of SF6 in porous silica. J. Chem. Phys. 1991, 95(5), 3710-3718. (13) Brodka, A.; Zerda, T. W. Properties of liquid acetone in silica pores: Molecular dynamics simulation. J. Chem. Phys. 1996, 104(16), 6319-6326. (14) Gelb, L. D.; Gubbins, K. E. Characterization of porous glasses: Simulation models, adsorption isotherms, and the BET analysis method. Langmuir 1998, 14, 2097-2111. (15) Steele, W. A.; Bojan, M. J. Simulation studies of sorption in model cylindrical micropores. AdV. Colloid Interface Sci. 1998, 76-77, 153-178. (16) Shi, W.; Johnson, J. K. Histogram reweighting and finite-size scaling study of the Lennard-Jones fluids. Fluid Phase Equilib. 2001, 187188, 171-191. (17) Gelb, L. D.; Gubbins, K. E. Simulations of capillary condensation in porous glasses. In AIChE Symposium Series #325, volume 97, pages 292-295. AIChE, 2001. (18) Spohr, E.; Trokhymchuk, A.; Henderson, D. Adsorption of water molecules in slit pores. J. Electroanal. Chem. 1998, 450, 281-287. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (20) Heffelfinger, G. S.; van Swol, F.; Gubbins, K. E. Liquid-vapor coexistence in a cylindrical pore. Mol. Phys. 1987, 61(6), 13811390. (21) Heffelfinger, G. S.; Tan, Z.; Gubbins, K. E.; Marini Bettolo Marconi, U.; van Swol, F. Fluid mixtures in narrow cylindrical pores: Computer simulation and theory. Int. J. Thermophys. 1988, 9(6), 1051-1060. (22) Connolly, M. L. Analytical molecular surface calculation. J. Appl. Crystallogr. 1983, 16, 548-558. (23) Connolly, M. L. Computation of molecular volume. J. Am. Chem. Soc. 1985, 107, 1118-1124.
NL025748P
1285
