6050
Langmuir 1999, 15, 6050-6059
Transport Diffusion of Oxygen-Nitrogen Mixtures in Graphite Pores: A Nonequilibrium Molecular Dynamics (NEMD) Study† Karl P. Travis*,‡ and Keith E. Gubbins Department of Chemical Engineering, Riddick Laboratories, North Carolina State University, Raleigh, North Carolina 27695-7905 Received October 19, 1998. In Final Form: January 10, 1999 The technique of dual control volume grand canonical molecular dynamics (DCV GCMD) is used to study the diffusivity of a mixture of nitrogen and oxygen in a graphite slit pore as a function of pore width. We find evidence in support of combined viscous and diffusive transport through these narrow slit pores. The viscous contribution to the flow becomes weaker as the pore width decreases. The fluid velocity profiles show evidence of microscopic slip but still retain a classical Navier-Stokes parabolic signature. The concentration profiles for each component in the mixture show an approximately linear variation with distance along the pore length, suggesting that cross-coupling effects are weak. We find that the diffusion of oxygen-nitrogen mixtures through a graphite pore displays a complex dependence upon the pore width. Molecular packing appears to play a very significant role in determining the flow of the mixture. Thermodynamic effects are of more importance in these simulations than sieving effects; hence, we do not find a greater diffusivity for oxygen compared to nitrogen.
I. Introduction Of the two categories of diffusion, self-diffusion and transport diffusion, transport diffusion is the most important phenomena from both a practical and theoretical point of view. Transport diffusion is defined as the transport of mass under a gradient in concentration (more strictly, a gradient in chemical potential). Self-diffusion, on the other hand, is the process of mass transport under conditions of macroscopic equilibrium. Transport diffusion is a collective property while self-diffusion is a single particle property. Porous materials are used extensively in the petroleum and chemical process industries as catalysts and adsorbents. Transport of fluid through these porous materials occurs mainly by diffusion, which often controls the overall rate of the process. A detailed understanding of the complexities of diffusional behavior in porous materials is therefore essential for the design of improved catalytic and adsorption processes. The separation of oxygen from a mixture of oxygen and nitrogen by passing a stream of air through carbon molecular sieves (CMS) is a good example of an important industrial process that is diffusion controlled.1 Oxygen selectivities of between 3 and 30 have been reported, even though the kinetic diameters of oxygen and nitrogen differ by less than 0.03 nm. The class of adsorbents known as carbon molecular sieves possess micropores of well-defined structure. The micropores are considered to be slits of approximately 0.5 nm in width between two graphitized carbon layer planes.2 Micropore diffusion of adsorbate molecules occurs by a * To whom correspondence should be addressed. † Presented at the Third International Symposium on Effects of Surface Heterogeneity in Adsorption and Catalysis on Solids, held in Poland, August 9-16, 1998. ‡ Present address: Department of Chemistry, Imperial College, London SW7 2AY, U.K. Tel: 44-(0)171-594-5851. Fax: 44-(0)171594-5804. (1) Ju¨ntgen, H.; Knoblauch, K.; Harder, K. Fuel 1981, 69, 817. (2) Chihara, K.; Suzuki, M.; Kawazoe, K. J. Colloid. Int. Sci. 1978, 64, 584.
process of hopping between energy minima in the adsorbent’s potential energy field.3 This type of process is an activated process and is expected to be temperature dependent. Apart from thermal considerations, many other parameters could affect the transport rates of fluids through CMS, pore size and entropic considerations4 being two such examples. It is important to be able to identify the key parameters in the diffusion process and then to find their optimum values in order to maximize the selectivity of one species over another in a mixture to be separated by kinetic means. Chihara and Suzuki5 attempted to vary the ratio of diffusivities of oxygen and nitrogen in CMS by adjusting the pore width. The pore width was varied in their experiments by adsorption of hydrocarbons followed by heat treatment. They concluded that the absolute diffusivities of oxygen and nitrogen can be decreased by 1 order of magnitude, but their ratio cannot be varied by this method. Seaton and co-workers6 have studied the separation of oxygen and nitrogen in model graphitic pores. They conducted molecular dynamics simulations of self-diffusion in individual pores and found that the diffusivities were strongly dependent on the pore width. Using a randomly etched graphite pore model (REGP), they found that the degree of kinetic separation observed experimentally can be reproduced at the level of individual pores. In a more recent publication the same authors7 looked at transport diffusion of oxygen and nitrogen in the same model pore system, concluding that pore length was a controlling factor in the separation mechanism. The development of the technique of dual control volume grand canonical molecular dynamics (DCV GCMD)8-10 has (3) Ka¨rger, J.; Ruthven, D. M. Diffusion in zeolites and other microporous solids; Wiley: New York, 1992. (4) Singh, A.; Koros, W. J. Ind. Eng. Chem. Res. 1996, 35, 1231. (5) Chihara, K.; Suzuki, M. Carbon 1979, 17, 339. (6) Seaton, N. A.; Friedman, S. P.; MacElroy, J. M. D.; Murphy, B. J. Langmuir 1997, 13, 1199. (7) MacElroy, J. M. D.; Friedman, S. P.; Seaton, N. A. Chem. Eng. Sci., submitted for publication. (8) Heffelfinger, G.; van Swol, F. J. Chem. Phys. 1994, 100, 7548. (9) MacElroy, J. M. D. J. Chem. Phys. 1994, 101, 5274.
10.1021/la981465u CCC: $18.00 © 1999 American Chemical Society Published on Web 03/27/1999
Diffusion of O2-N2 Mixtures in Graphite Pores
enabled the phenomenon of transport diffusion to be simulated directly, with conditions close to those in a real diffusion experiment, in which gradients in the chemical potential drive the flow and pressure gradients can be present. In the DCV GCMD technique, a gradient in chemical potential is established by placing two particle reservoirs at either end of a single pore and maintaining them at fixed, but different, chemical potentials. This is achieved by periodically conducting a series of creations and deletions according to the prescription of grand canonical Monte Carlo (GCMC). By placement of the particle reservoirs within the pore, surface barrier effects can be eliminated. This method has been used to study the transport diffusion of methane in graphite,10 diffusion of various gases in zeolite frameworks,11 diffusion through polymer membranes,12 and diffusion of a mixture of oxygen and nitrogen through a graphite slit-pore.7,13 Previous simulation studies on the transport diffusion of oxygen and nitrogen have been restricted to the pure components. In this work we present results from simulations of transport diffusion of a binary mixture of oxygen and nitrogen through a graphite slit-pore, as a function of pore width at two different temperatures for a single composition. Unlike related published work, we allow for interactions between fluid molecules. This work is the first phase of a longer term study of the molecular nature of the separation mechanism of oxygen and nitrogen in carbon micropores. We ignore pore networking effects and also any pore entry effects, concentrating on the influence of pore width and temperature alone. II. Phenomenological Equations For isothermal mass transport in a simple, binary, bulk fluid the usual nonequilibrium thermodynamic treatment (see for example ref 1414) leads to a single transport law involving a single diffusion coefficient. A treatment of diffusion through porous membranes has been developed by Mason and co-workers.15,16 The reason for the different treatment of porous systems compared to that of bulk systems arises from one of the postulates of nonequilibrium thermodynamics, namely that of local thermodynamic equilibrium. For local thermodynamic equilibrium to hold in a given volume element, the thermodynamic state variables should not vary appreciably with position or time on a microscopic scale. If one ignores time-dependent phenomena, a typical volume element in a microporous system must contain several molecular mean free paths for the postulate of local thermodynamic equilibrium to hold. Such a situation would be realized for a high-density gas in a membrane. On the other hand, if the gas is rarefied, a typical volume element must now include the membrane particles as well as the adsorbate molecules. These two examples may be considered to be limiting cases which give rise to viscous transport and diffusive transport, respectively. Typically, both of these transport mechanisms occur simultaneously, but independently, since, according to Curie’s principle, there are no couplings between viscous phenomena and diffusive phenomena due (10) Cracknell, R. F.; Nicholson, D.; Quirke, N. Phys. Rev. Lett. 1995, 74, 2463. (11) Pohl, P. I.; Heffelfinger, G. S.; Smith, D. M. Mol. Phys. 1996, 89, 1725. (12) Ford, D. M.; Heffelfinger, G. S. Mol. Phys. 1998 94, 673. (13) Travis, K. P.; Gubbins, K. E. Fundamentals of Adsorption, 6th; Meunier, F., Ed.; Elsevier: Paris, 1998; pp 1161-1166. (14) deGroot, S. R.; Mazur, P. Non-equilibrium thermodynamics; Dover: New York, 1953. (15) Mason, E. A.; Viehland, L. A. J. Chem. Phys. 1978, 68, 3562. (16) Mason, E. A.; Malinauskas, A. P. Gas transport in porous media: the dusty gas model. Elsevier: Amsterdam, 1983.
Langmuir, Vol. 15, No. 18, 1999 6051
to their different tensorial character. For binary mixtures in porous membranes the treatment of Mason and coworkers leads to the following expressions for the laboratory frame fluxes of the two fluid components:
J1 ) -L11 ∇µ1 - L12 ∇µ2 -
B0n1 ∇p η
(1)
J2 ) -L21 ∇µ1 - L22 ∇µ2 -
B0n2 ∇p η
(2)
Here the coefficients Lij are the linear phenomenological transport coefficients arising from diffusion flows, ∇µi is the gradient in chemical potential of species i, B0 is a parameter characteristic of the membrane, ni is the number density of species i, η is the viscosity of the fluid, and p is the hydrostatic pressure. Strictly speaking, these flux expressions apply only to fluids composed of structureless molecules. For molecules with internal structure we should allow for the possibility of coupling between gradients in chemical potential and the curl of the angular velocity field. In the simulations described in this paper the average angular velocity is small and has no appreciable positional dependence. These coupling terms are therefore negligible, so that we are justified in using eqs 1 and 2 for diatomic molecules. The third term on the right side of eqs 1 and 2 arises from Poiseuille flow, which occurs whenever a fluid is forced through a membrane by a pressure gradient. The Navier-Stokes equations can be solved for a system of uniaxial molecules undergoing pressure driven flow at low Reynolds number between two infinite parallel plates. The solutions for the center-of-mass velocity and angular velocity are respectively17
ux(z) )
[
2ηr -H2 dp 1 - zj2 + 8η dx (η + ηr)KH
(
cosh(KHzj) -1 cosh(KH) sinh(KHzj) -H zj ωy(z) ) 4η sinh(KH) coth(KH)
[
]
)]
(3) (4)
with zj ) z/(H/2), K ) [ηηr/(η + ηr)(ζ + ζrr)]1/2, where η and ηr are the shear and vortex viscosities, ζ and ζrr are the spin-shear and spin-vortex viscosities, and H is the pore width. If ηr is very small, the classical Navier-Stokes profile is recovered:
ux(z) )
-H2 dp [1 - zj2] 8η dx
(5)
The total average flux in the x-direction is obtained by integrating the product of ux(z) and the fluid number density over half the pore width. Assuming that the density is independent of z and that eq 5 gives the form of the streaming velocity, the total viscous flux is Jv ) (B0n/ η)∇p, with B0 ) H2/12. In reality, the concentration depends significantly upon z in narrow pores, and hence, this is a very simplistic approximation. However we note that in this crude approximation the viscous flux varies as the second power of the pore width, and hence, it is expected to become insignificant in narrow pores such as those we consider in this work. Furthermore, in very narrow pores the validity of the Navier-Stokes equations has been called into question. Simulation studies of planar (17) Eringen, A. C. Contributions to mechanics; Abir, D., Ed.; Pergamon: Oxford, U.K., 1969.
6052 Langmuir, Vol. 15, No. 18, 1999
Travis and Gubbins
Poiseuille flow in pores as narrow as 4 molecular diameters show that the shear viscosity is probably a nonlocal function of position.18,19 Equations 1 and 2 may be simplified further if the fluxes arising from cross-coupling effects are negligible compared to the fluxes -L11∇µ1 and -L22∇µ1. A recent simulation study of transport diffusion of hydrogen-methane mixtures in a graphite pore found that, in that system, crosscoupling effects were immeasurably small.20 Using these approximations, we rewrite eqs 1 and 2 as
J1 ) -Deff 1 ∇c1
(6)
-Deff 2 ∇c2
(7)
J2 )
where in addition to the previously mentioned approximations we have converted the chemical potential gradients into concentration gradients in order to express the flux equations in the more familiar Fickian form (here ci is the concentration of species i in units of mole per cubic meter). The conversion of the gradient term introduces a thermodynamic factor which has been absorbed eff into the effective diffusion coefficients Deff 1 and D2 . Note that c1 and c2 refer to concentrations within the pore and not the concentrations of the two species as they exist in the bulk external gas phase. If these effective diffusion coefficients do not depend appreciably on the local concentration, we would expect to obtain linear concentration profiles for each species if cross coupling effects are indeed negligible. If Poiseuille flow is nonnegligible, the Deff’s will be influenced by the viscous flow. Aside from the effective diffusivities, another quantity useful for describing the transport of material through a membrane is the permeability, Fi:
Jxi Fi ) ∆pi/Lx
(8)
where pi is the partial pressure of species i in the bulk external gas phase and Lx is the length of the pore in the x-direction. The ratio of species permeabilities defines a dynamic separation factor. III. Simulation Method 1. Model Description. We represent the interaction of the adsorbate molecules with the graphite planes in our model by a smooth 10-4-3 potential due to Steele.21 Neglecting the corrugations in the xy graphite planes is expected to be a reasonable approximation for this study but might be a serious omission at lower temperatures. This potential is a function of the z-coordinate only. To represent opposing graphite planes in a slit-pore we employ a superposition of two potential functions
[(
Φic ) 4πicσic2∆nc
)
σic 1 5 (H/2) - z
σic4
10
-
(
(
σic 1 2 (H/2) - z
)
)
4
-
10 σic 1 + 6∆((H/2) - z + 0.61∆)3 5 (H/2) + z 4 σic4 σic 1 (9) 2 (H/2) + z 6∆((H/2) + z + 0.61∆)3
(
)
]
(18) Travis, K. P.; Todd, B. D.; Evans, D. J. Phys. Rev. E 1997, 55, 4288. (19) Travis, K. P.; Gubbins, K. E. Manuscript in preparation.
Figure 1. Schematic diagram of the simulation cell for DCV GCMD. Flow takes place in the x-direction from the source control volume to the sink control volume. Table 1. Lennard-Jones Potential Parameters for Oxygen, Nitrogen, and Carbon Used in This Work atom
σ (nm)
/kB (K)
carbon nitrogen oxygen
0.340 0.3296 0.2940
28.0 60.39 75.49
where ∆ is the inner layer spacing in graphite, which is taken to be 0.335 nm, nc ) 114 nm-3 is the carbon atom number density in graphite, and H is the pore width, defined as the distance between the centers-of-mass of the innermost graphite planes, while ic and σic are the Lennard-Jones parameters appropriate for interactions between a molecular site of species i and a carbon atom. Oxygen and nitrogen are modeled as two-center Lennard-Jones molecules with fixed bond lengths. Interactions between sites on different fluid molecules are modeled with a truncated and shifted Lennard-Jones 12-6 potential:
φij(r) )
{
[( ) ( ) ]
4ij 0
σij r
12
-
σij r
6
- φij(rc)
r e rc r > rc
(10)
Here r is the scalar interatomic distance between a pair of interacting sites, rc is the truncation distance, and φij(rc) is the value of the potential energy at the point of truncation. Lennard-Jones potential parameters for nitrogen, oxygen, and carbon were taken from the literature22 and are given in Table 1. Parameters appropriate for interactions between chemically different species, for example, between an oxygen atom and a nitrogen atom, are given by the Lorentz-Berthelot mixing rules: σij ) 1/ (σ + σ ) and ) x . The bond lengths of the molecules 2 i j ij i j are 0.1097 nm for nitrogen and 0.1169 nm for oxygen. We truncate the Lennard-Jones potential at rc ) 2.5σij. Note that we do not truncate or shift the Steele 10-4-3 potential. We consider only classical dynamics in constructing our model, quantum effects being unimportant for a system such as ours.23 We do not take account of the quadrupole for nitrogen molecules. Justification for this approximation comes from simulation studies of nitrogen adsorption in slit-pores at ambient temperature, which found that the quadrupole had no significant effect on the results.24 (20) MacElroy, J. M. D.; Boyle, M. J. Manuscript in preparation. (21) Steele, W. A. The interaction of gases with solid surfaces; Pergamon: Oxford, U.K., 1974. (22) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230. (23) Beenakker, J. J. M.; Borman, V. D.; Krylov, S. Y. Chem. Phys. Letts. 1995, 232, 379. (24) Kaneko, K.; Cracknell, R. F.; Nicholson, D. Langmuir 1994, 10, 4606.
Diffusion of O2-N2 Mixtures in Graphite Pores
Langmuir, Vol. 15, No. 18, 1999 6053
Control volumes are placed at each end of the slit-pore. Placing the control volumes inside the pore eliminates pore-entrance effects and greatly simplifies the interpretation of our results. The importance of pore-entrance effects on the diffusion mechanism will be dealt with in subsequent publications. We define our coordinate system such that the flow is in the x -direction and the graphite planes are separated along the z -direction. The volume of the control volumes is taken to be the same as that of the flow region between them (see Figure 1). We note that there is no unambiguous definition of volume in a porous membrane. However for simplicity, we define the volume of the flow region and control volumes to be V ) HLxLy, which contains a certain amount of dead space due to the implicit carbon atoms in graphite. The control volumes and the flow region have a length in the x-direction of Lx ) 9.888 nm, while, in the z-direction, the length varied from H ) 9.888 nm to H ) 24.720 nm depending on the pore width. The length in the y-direction varied from Ly ) 9.888 nm to Ly ) 16.48 nm depending on pore width. There is a tradeoff between having too many molecules, which is computationally expensive, and simulation accuracy. The smallest system size we used had a total of 1225 molecules, while the largest had 3642. Larger system sizes were used for the narrower pores. We studied the diffusion at a number of pore widths: {H/∆} ) {5.00, 4.00, 3.00, 2.50}, where ∆ is the graphite layer spacing. 2. Calculation of Hydrodynamic Quantities. In the spirit of the Irving-Kirkwood method25 we can define instantaneous microscopic expressions for the densities of a number of hydrodynamic quantities. The local molar concentration of component γ at position r is given by
cγ(r,t) )
Nγ
1 NA
δ(r - ri(t)) ∑ i)1
(11)
where the summation runs over the members of component γ, δ(r - ri(t)) is the Dirac delta function, ri(t) is the instantaneous center-of-mass location of molecule i, and NA is Avogadro’s constant. Similarly, the molar flux (number of moles per unit area) of component γ can be defined as
cγuγ(r,t) )
Nγ
1
∑ viδ(r - ri(t))
NA i)1
(12)
where vi is the center-of-mass velocity of molecule i. The total fluid concentration and total molar flux follow from eqs 11 and 12 by summing over γ,
c(r,t) ) J(r,t) )
∑γ
coordinate directions. The system we consider is spatially homogeneous in the y-direction. 3. DCV GCMD Algorithm. The DCV GCMD algorithm consists of performing numerous cycles, each of which comprises a molecular dynamics step, in which the trajectories of all fluid atomic sites are incremented, followed by a series of grand canonical Monte Carlo creations and destructions of either species in each of the two control volumes. A creation attempt of a molecule of species γ in control volume c is accepted with a probability, min[1, exp[βµ′γ(c) - β∆νγ + ln(V(c)/(Nγ + 1))]]
(15)
where V(c), Nγ, and µ′γ(c) are the volume, the current number of molecules of species γ in control volume c and the residual chemical potential of species γ in control volume c, respectively, while ∆νγ is the energy change accompanying the insertion of a molecule of species γ into the control volume, and β ) 1/kBT. The residual chemical potential is related to the full chemical potential, µγ, via βµ′γ ) βµγ + lnΛγ-3 where Λγ is the thermal de Broglie wavelength for a molecule of species γ. Destructions of molecules of species γ are accepted with a probability
min[1, exp[-βµ′γ(c) - β∆νγ + ln(Nγ/V(c))]] (16) where ∆νγ is now the energy change accompanying the destruction of a molecule of species γ from control volume c. When molecules are created in either control volume, they are assigned thermal components of both translational velocity and angular velocity selected from a Maxwell-Boltzmann distribution. Furthermore, newly created molecules are given an appropriate initial component of streaming velocity to ensure creations are compatible with mass transport (further details are given below). The system is kept isothermal via the following scheme: The simulation box is divided into 21 bins of equal volume along the flow direction. The local temperature in each bin is then controlled by means of a Nose´Hoover thermostat.26,27 It is important to note that a Nose´Hoover thermostat is to be preferred over a Gaussian thermostat,28 because the latter algorithm requires an explicit form for the substantial time derivative of the component streaming velocities, which may not be zero, even at steady state. Since an explicit form for the local streaming velocities is not generally known, one cannot obtain a expression for its substantial time derivative. The local translational temperature in each bin is defined by Nbin
cγ(r,t)
∑γ cγuγ(r,t)
〈
(13) (14)
In practice we replace the Dirac delta function appearing in eqs 11-14 by a function which takes the value of unity within some finite range of r but is zero otherwise. This procedure is often referred to as the “bins” method. In the simulations, we divide up the simulation box into a number of bins along the x-direction and also along the z-direction. In this way, we can resolve the quantities in both of these (25) Irving, J. H.; Kirkwood, J. G. J. Chem. Phys. 1950, 18, 817.
〈T(xbin)〉 )
∑
[pi - Miuγ(xbin)][pi - Miuγ(xbin)]〉
i∈bin
〈3Nbin - d(Nbin/N)〉
(17)
where uγ(xi) is the local average streaming velocity of species γ in bin i, Nbin is the number of molecules in a particular bin, N is the total number of molecules of that species, and d is the number of parameters that are used in the least-squares fit of the molecular velocities to a (26) Nose, S. J. Chem. Phys. 1984, 81, 511. (27) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (28) Evans, D. J.; Morriss, G. P. Statistical mechanics of nonequilibrium liquids; Academic Press: London, 1990.
6054 Langmuir, Vol. 15, No. 18, 1999
given functional form. In this work we have fitted the streaming velocities to a two-parameter linear equation, and therefore, d ) 2. The factor (Nbin/N) then accounts for the fraction of d alloted to each bin. The presence of an adsorbent interaction field means that the streaming velocities will have a dependence on both the x and z coordinates. However, for simplicity, we average over the z coordinate, leaving only the x dependence. The local average streaming velocities within each bin are calculated using eqs 11 and 12. To ensure that molecular creations are compatible with mass transport, a component of streaming velocity is added to the thermal velocity of newly created molecules. We determine this component of streaming velocity by taking a value for the flux at the center of the flow region and dividing this by the concentration in the appropriate control volume. Since we begin the simulation with zero molecules, initially there will be no measurable flux. As molecules begin to diffuse through the pore, a steadystate flux will gradually develop, which is constant at any yz plane in the simulation cell. To prevent the streaming velocity from diverging, we follow Cracknell et al.10 and introduce a degree of course graining. The method of augmenting the molecular velocities upon creation proceeds as follows: We allow a certain settling time, say 50 000 steps, for obtaining sufficient molecules in the flow region to obtain a flux. After this time, we divide the flux (averaged over the settling time) by the average concentration of that species in the control volume of interest. The streaming velocity thus obtained is then added to the thermal component of velocity of newly created molecules in that control volume. After the settling time, we average the flux over periods of 1000 molecular dynamics steps and use this value in the subsequent calculations of the streaming velocity in the end sections. Because we employ a smooth potential to represent the graphite planes in a real carbon pore, we need to account for the exchange of momentum which would take place between fluid molecules and the carbon atoms in a real adsorbent. We employ the so-called diffuse boundary condition, based on the diffuse boundary conditions used by Cracknell et al.10 but modified here for the case of molecules. Application of our molecular diffuse boundary condition proceeds as follows. After each molecular dynamics time step we check to see if the following two conditions are satisfied: (1) The center-of-mass momentum (in the laboratory frame) of a given molecule in the z-direction has reversed in sign. (2) The center-of-mass of that same molecule is within the repulsive region of the Steele 10-4-3 potential. If, and only if, both of these conditions are satisfied, we reselect the center-of-mass momentum of that molecule in the directions parallel to the confining surface from a Maxwell-Boltzmann distribution at the appropriate temperature. Note that we do not alter the angular velocity of the molecules. This is a reasonable approximation for molecules which are roughly spherical like oxygen and nitrogen, diffusing at moderate to high temperatures. One important control parameter in DCV GCMD is the ratio of stochastic to dynamic steps, nMC/nMD. An optimum value for nMC/nMD was arrived at by finding the smallest value which still yielded the correct concentrations of the two species obtained from conducting separate grand canonical adsorption simulations at the relevant thermodynamic state point for the two control volumes. Using too large a value for this parameter greatly increases the CPU time needed to reach a nonequilibrium steady state. We find that, for our simulations, a value for nMC/nMD of 50 is optimum.
Travis and Gubbins Table 2. Residual Chemical Potentials Used in DCV GCMD Simulationsa t ) 25 °C
t ) 0 °C
species
µ′source (kJ mol-1)
µ′sink (kJ mol-1)
µ′source (kJ mol-1)
µ′sink (kJ mol-1)
nitrogen oxygen
-11.3507 -14.7883
-12.3681 -15.7875
-10.2099 -13.3562
-11.1241 -14.2703
a These values are obtained from separate bulk isobaricisothermal Monte Carlo simulations employing Widom insertions for an 80:20 mixture of nitrogen and oxygen.
Figure 2. Concentration profiles across the pore plotted against the scaled z coordinate. Solid lines represent data at t ) 0 °C, and dashed lines are profiles at t ) 25 °C.
In a molecular dynamics step, a fifth-order Gear algorithm29 was used to integrate Hamilton’s equations of motion, supplemented with a Nose´-Hoover thermostating scheme26,27 and Gaussian constraints28,30 to fix the bond lengths. The integration time step was chosen to be 2.5 fs. Periodic boundary conditions operate in the y -direction. There are no periodic boundary conditions in the x -direction. The ends of the simulation box are dissolving boundaries; that is, if a molecule reaches either of these boundaries, it is removed from the simulation. 4. Operating Conditions. The simulations were performed for a mixture of oxygen and nitrogen at a single composition of 80 mol % nitrogen. The pressures of the bulk gases in equilibrium with the fluid in the control volumes are taken to be 1.5 and 1.0 MPa for the source and sink control volumes, respectively. These pressures were chosen to obtain fluxes with a reasonable signal to noise ratio. The residual chemical potentials of each component at the two pressures were obtained by carrying out separate bulk isothermal-isobaric Monte Carlo simulations employing Widom insertions.31,32 The set of residual chemical potentials appropriate to the two temperatures studied: t ) 0 and 25 °C are given in Table 2. IV. Results and Discussion The fluid concentration profiles as a function of the scaled pore width are shown in Figure 2 for the two temperatures studied: t )25 (dashed lines) and t ) 0 °C (solid lines). Each of the profiles is characterized by the presence of two strong maxima which correspond to the (29) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press: Oxford, U.K., 1987. (30) Daivis, P. J.; Evans, D. J.; Morriss, G. P. J. Chem. Phys. 1992, 97, 616. (31) Widom, B. J. Chem. Phys. 1963, 39, 2802. (32) Frenkel, D.; Smit, B. Understanding molecular simulation; Academic Press: San Diego, CA, 1996.
Diffusion of O2-N2 Mixtures in Graphite Pores
Langmuir, Vol. 15, No. 18, 1999 6055
Table 3. Data Obtained from the t ) 25 °C Simulationsa H/∆
〈c〉 (mmol cm-3)
JN2 (mmol mm-2 s-1)
JO2 (mmol mm-2 s-1)
dcN2/dx (mol mm-4)
dcO2/dx (mol mm-4)
FN2 (10-2 mmol m-1 bar-1 s-1)
FO2 (10-2 mmol m-1 bar-1 s-1)
eff (10-3 DN 2 cm2 s-1)
eff (10-3 DO 2 cm2 s-1)
2.5 3.0 4.0 5.0
9.76 9.49 5.69 4.14
24.33(25) 15.21(37) 12.91(50) 9.37(43)
4.53(19) 3.91(19) 2.36(06) 2.42(12)
-0.0397(10) -0.1279(11) -0.1444(18) -0.1045(13)
-0.0200(8) -0.0426(7) -0.0294(7) -0.0246(7)
6.01(6) 3.76(9) 3.2(1) 2.3(1)
4.48(19) 3.87(19) 2.33(6) 2.39(12)
6.1(2) 1.19(3) 0.89(4) 0.90(4)
2.3(1) 0.92(5) 0.80(3) 0.98(6)
a
Fluxes are evaluated in the laboratory frame. The figures in parentheses are the statistical uncertainties in the last digits. Table 4. Data Obtained from the t ) 0 °C Simulationsa
H/∆
〈c〉 (mmol cm-3)
JN2 (mmol mm-2 s-1)
JO2 (mmol mm-2 s-1)
dcN2/dx (mol mm-4)
dcO2/dx (mol mm-4)
FN2 (10-2 mmol m-1 bar-1 s-1)
FO2 (10-2 mmol m-1 bar-1 s-1)
eff (10-3 DN 2 cm2 s-1)
eff (10-3 DO 2 cm2 s-1)
2.5 3.0 4.0 5.0
10.91 12.32 8.07 5.76
23.34(62) 14.39(50) 15.84(56) 11.66(50)
4.88(19) 4.10(19) 3.48(19) 2.81(19)
-0.017(1) -0.0821(1) -0.146(1) -0.1147(9)
-0.039(1) -0.050(2) -0.046(1) -0.0363(8)
5.77(15) 3.56(12) 3.92(14) 2.88(12)
4.83(19) 4.05(19) 3.44(19) 2.78(19)
13.7(9) 1.75(6) 1.09(4) 1.08(4)
1.25(6) 0.82(5) 0.76(4) 0.77(6)
a
Fluxes are evaluated in the laboratory frame. The figures in parentheses are the statistical uncertainties in the last digits.
adsorbed contact layers. The height of these maxima increases with decreasing temperature, indicating greater adsorption at the lower temperature. This phenomenon can be explained in terms of the increasing dominance of potential effects over kinetic effects. Figure 2 shows evidence of a weak layering of the fluid between the contact layers. At H ) 5∆ a total of 4 fluid layers can be discerned, while at H ) 4∆ three layers are now present. At H ) 3∆, only two layers are in evidence, while at H ) 2.5∆ the concentration profile contains two maxima that are strongly overlapping. For the largest three pore widths, the number of fluid layers seems to follow the rule H/∆-1, while at the lowest pore width, H/∆ is not an integer but is midway between H/∆ ) 3 and H/∆ ) 2, with the result that the maxima overlap but have not completely merged to form a single layer. All of the concentration profiles fall to zero at a point away from the physical location of the innermost graphite planes, indicating dead space between the fluid and the walls. The increased dead space results from overlap of the two wall interaction potentials which produces a larger repulsive region as well as a deeper but narrow attractive well. The fraction of dead space or excluded volume increases as the pore width decreases, such that at the lowest pore width less than half of the pore space is available to the fluid. Column 2 of Tables 3 and 4 gives the mean concentrations for each of the pore widths at the two temperatures. These mean concentrations were obtained by integrating the concentration profiles in Figure 2 over the range [0,1]. From the tabulated values we see that the mean concentration increases with decreasing pore width, with the exception of the lowtemperature data, where the concentration decreases in going from H ) 3∆ to H ) 2.5∆. The optimum pore width for adsorption must be H ) 3∆. Below this width, excluded volume effects ensure that the fluid is virtually onedimensional and so fewer molecules are able to fit into the available pore space. In Figure 3 we show the streaming velocity as a function of the scaled pore width for H ) 2.5∆ and H ) 5∆ at the two temperatures: t ) 0 °C and t ) 25 °C. All the streaming velocity profiles display evidence of microscopic slip; that is, the velocity extrapolates to zero at a point inside the fluid. The point at which they do extrapolate to zero (noslip plane) coincides with the point at which the concentration profile extrapolates to zero. Apart from some structure near the no-slip planes, these profiles have clear quadratic character. We cannot attach too much significance to the peaks at the no-slip planes because the data in this region exhibit very poor statistics (a consequence of the very small number of molecules occupying the bins
Figure 3. Center-of-mass flow velocity plotted against the scaled z coordinate. The solid line represents data from the H ) 5∆ simulation for t ) 25 °C; the broken line is for the same pore width but for t ) 0 °C. The circles represent data for H ) 2.5∆ and t ) 25 °C while the diamonds are for H ) 2.5∆ and t ) 0 °C. The smooth curves are least-squares fits of the data to u(z) ) a + bz2.
centered at these points). The existence of parabolic velocity profiles indicates that we have combined transport taking place, i.e., a combination of diffusion and Poiseuille flow. The fact that the velocity profiles can be fitted by a parabolic equation (dash-dot lines in Figure 3) also confirms that we are in a regime in which the mean angular velocities are effectively zero, i.e., the fluid is conforming to the classical hydrodynamic picture. Despite our running of the simulations for very long times, the velocity profiles are quite noisy, a fact that points to the weakness of the viscous mode of transport in comparison to the diffusive transport. Recent results for methane transport in carbon pores support this conclusion.33 The center-of-mass velocity of the fluid decreases with both decreasing pore width and temperature. The pore width induced lowering of the flow speed is of much greater significance than the temperature induced lowering effect. This can be explained by the fact that the classical NavierStokes equation for the flow velocity varies as H2 whereas the root mean square speed of the molecules varies only as the square root of temperature. Figures 4 and 5 show nitrogen concentration profiles along the flow direction at t ) 25 °C and t ) 0 °C, respectively. The concentration profiles in the flow region (x ) -5 to x ) +5 nm) show a linear variation with distance, (33) Travis, K. P.; Gubbins, K. E. Unpublished work, 1998.
6056 Langmuir, Vol. 15, No. 18, 1999
Figure 4. Nitrogen concentration profiles parallel to the x -axis (flow direction) at t ) 25 °C for 4 different pore widths.
Figure 5. Nitrogen concentration profiles parallel to the x -axis (flow direction) at t ) 0 °C for 4 different pore widths.
in accordance with Fick’s law. Clearly, for nitrogen at these conditions, the coefficient of transport diffusion is only weakly dependent on the local concentration. The gradient of nitrogen concentration at H ) 2.5∆ is very weak in comparison to the gradient at the other pore widths. This weak concentration gradient is a direct result of the lower concentration of nitrogen in the source reservoir. Tables 3 and 4 list the concentration gradients obtained from linear regression of the data in the flow regions. From the tabulated data we see that the nitrogen concentration gradients increase in going from H ) 5∆ to H ) 4∆ but then decrease for the lowest two pore widths. The concentration of nitrogen generally increases as pore width decreasessa consequence of greater adsorption. However, at the lowest temperature, the nitrogen concentration is lower at the lowest pore width than it is for H ) 3∆. This may be a manifestation of different molecular packing in the narrow pore at the lower temperature. The near Fickian behavior displayed by the nitrogen concentration profiles suggests that cross-coupling effects (mass flow of one species caused by a gradient in the concentration of the other species) are negligible in these simulations. Figures 6 and 7 show oxygen concentration profiles along the flow direction at t ) 25 °C and t ) 0 °C, respectively. Despite the large degree of statistical noise present in these profiles, approximately Fickian concentration profiles are again observed. The increased statistical noise in the oxygen concentration profiles simply reflects the smaller numbers of oxygen molecules present in the mixture. The oxygen concentration increases in both
Travis and Gubbins
Figure 6. Oxygen concentration profiles parallel to the x-axis (flow direction) at t ) 25 °C for 4 different pore widths.
Figure 7. Oxygen concentration profiles parallel to the x-axis (flow direction) at t ) 0 °C for 4 different pore widths.
control volumes with decreasing pore width, with the exception of the lowest pore width, at which point it decreases relative to the H ) 3∆ case. The concentration of oxygen in the control volumes is greater at t ) 0 °C than it is at the higher temperature for all of the pore widths studied. The dips in density present at the extremities of the density profiles in Figures 4-7 are a consequence of using dissolving boundary conditions in the x-direction. Molecules created in the region close to the boundary will most likely be destroyed or go to oblivion by passing through the end face, hence the region of low density. The oxygen profiles at H ) 2∆ in Figures 6 and 7 display density peaks at the extremities. This behavior is most likely a consequence of the local fluid structure favoring creation of oxygen molecules in the vicinity of the dissolving end face. We have observed this idiosyncrasy in DCV GCMD simulations using atomic fluids with purely repulsive fluid-fluid interactions. In any case, the fluid behavior at the extremities has no bearing on the fluid properties in the important flow region of the simulation cell. Tables 3 and 4 list the concentration gradients obtained from a linear regression of the data in the flow region. The tabulated data show a steady increase in the concentration gradient with decreasing pore width with the exception of the lowest pore width, at which point the gradient weakens. The gradients are stronger at all pore widths for the lower temperature. The permeability of each species in the mixture is shown in Figure 8 as a function of pore width for the two temperatures studied. Although the permeability is a nonmonotonic function of pore width, there is a general
Diffusion of O2-N2 Mixtures in Graphite Pores
Figure 8. Permeabilities for nitrogen and oxygen as a function of pore width at the two temperatures studied: t ) 0 °C and t ) 25 °C.
Langmuir, Vol. 15, No. 18, 1999 6057
temperature, apart from the highest pore width. At H ) 3∆, the ratio is less than unity, implying a dynamic separation in favor of oxygen whereas, at H ) 2.5∆ and H ) 4∆, nitrogen is the preferred component. The result at H ) 5∆ appears to be temperature dependent, in that the lower temperature favors nitrogen selectivity whereas the higher temperature favors oxygen selectivity. We cannot attach too much significance to this feature because of the large statistical uncertainties in the separtion factors at the H ) 5∆ pore width. The apparent difference between the dynamic separation factors at H ) 3∆ and H ) 5∆ is not statistically significant either. The simulation times would need to be increased by at least a factor of 4 to fully resolve these small differences in dynamic separation factors. It is useful to compare these dynamic separation ratios with the equilibrium selectivities in each of the control volumes. The equilibrium selectivity in a given control volume is defined as eq SN ) 2,O2
Figure 9. Permeability ratios as a function of pore width at the two temperatures studied: t ) 0 °C and t ) 25 °C.
trend toward an increase in permeability with decreasing pore width. Since the partial pressure drops are constant in all simulations, the permeability is merely proportional to the flux of a given species. An increase in permeability with decreasing pore width is expected since the concentration of particles generally increases with decreasing pore width (species fluxes and permeabilities are collected in Tables 3 and 4). However, the increase in flux in going from the H ) 3∆ to the H ) 2.5∆ pore width is much greater than expected from the associated change in concentration. The reason for the large change in permeability could be due to the less tortuous route taken by the diffusing molecules at the narrower pore width. At this pore width, the interaction potential forces the molecules to flow in almost single file. Due to the statistical uncertainties associated with the species permeabilities, the differences in these quantities at the two temperatures have no statistical significance for the lowest two pore widths. The slight increase in oxygen permeability at H ) 5∆ and t ) 25 °C has no statistical significance for the same reason. We define a “dynamic separation factor” as the ratio of nitrogen permeability to that of oxygen. Alternatively, we may define another measure of dynamic separation by the ratio of effective diffusivities. Figure 9 shows the permeability ratio FN2/FO2 plotted against pore width. The permeability ratio is a nonmonotonic function of pore width. As pore width decreases the ratio first increases and then decreases before increasing again at the lowest pore width. The ratio becomes larger at the higher
xN2/xO2 yN2/yO2
(18)
where xi is the mole fraction of species i in the control volume of interest (measured in the DCV GCMD simulations) and yi is the mole fraction of that species in the bulk eq are collected in Table external gas phase. Values of SN 2,O2 5. Figure 10 shows the equilibrium selectivities for the source and sink control volumes (high-pressure and lowpressure regions, respectively) as a function of pore width at t ) 25 °C and t ) 0 °C. The figure shows that nitrogen adsorption is favored over oxygen adsorption for all pore widths except H ) 3∆, where oxygen adsorption is greater (with the exception of the sink reservoir at t ) 25 °C). Nitrogen selectivity is always greater in the sink reservoir than in the source so we conclude that nitrogen adsorption may increase relative to oxygen adsorption as the pressure of the external gas phase is lowered. We are currently conducting a series of adsorption simulations to verify this conclusion. For a given control volume at a given pore width, lowering the temperature lowers the nitrogen selectivity. Comparing Figures 9 and 10, we see qualitatively similar behavior of the two selectivity measures, one dynamic and the other static. Quantitatively, the dynamic separation ratio is not too different from the equilibrium selectivity, and so we conclude that the favorable increase in either species shown in Figure 9 is merely due to thermodynamic effects. The reason the two selectivities oscillate with pore width is presumably a packing effect. We note that while basing separation factors on packing effects is good for simulations, it is inconclusive for CMS materials. For each change in H, there is a change in the fluid structure since one less fluid layer can be accommodated. At H ) 5∆ and H ) 3∆, oxygen molecules must be able to pack in a more energetically favorable arrangement while the opposite is true for H ) 4∆ and H ) 2.5∆. In Figure 11 we show the component effective diffusivities plotted against pore width for the two temperatures studied. Lowering the pore width has only a minor effect on this quantity until H ) 3∆ is reached, at which point the effective diffusivity of nitrogen increases dramatically at both temperatures while a smaller increase is observed for oxygen. For nitrogen, lowering the temperature increases the effective diffusivity at H ) 2.5∆ while the opposite is true for oxygen. While the mobility of a molecule is expected to decrease with decreasing temperature, the Fickian diffusivity may increase with
6058 Langmuir, Vol. 15, No. 18, 1999
Travis and Gubbins
Table 5. Permeability Ratios, Effective Diffusivity Ratios, and Equilibrium Nitrogen Selectivitiesa t ) 25 °C H/∆
FN2/FO2
eff eff /DO DN 2 2
2.5 3.0 4.0 5.0
1.34(6) 0.97(5) 1.37(6) 0.96(6)
2.7(1) 1.29(8) 1.11(6) 0.92(7)
a
t ) 0 °C
eq (source) SN 2,O2
eq SN (sink) 2,O2
1.13 0.985 1.119 1.086
1.235 1.067 1.144 1.105
FN2/FO2
eff eff /DO DN 2 2
eq (source) SN 2,O2
eq SN (sink) 2,O2
1.19(6) 0.88(5) 1.14(7) 1.04(8)
11(1) 2.1(1) 1.43(9) 1.4(1)
0.982 0.865 1.072 1.056
1.123 0.957 1.128 1.092
The figures in parentheses are the statistical uncertainties in the last digits.
Figure 10. Equilibrium nitrogen selectivity in each control volume plotted as a function of pore width at the two temperatures studied: t ) 0 °C and t ) 25 °C.
Figure 12. Effective diffusivity ratios as a function of pore width at the two temperatures studied: t ) 0 °C and t ) 25 °C.
between a nitrogen molecule and an oxygen molecule to be critical. For such small widths we would expect the oxygen diffusivity to be greater than the nitrogen diffusivity (sieving regime). The increasing importance of crosscoupling effects cannot be ruled out; however, we have no way of determining them from this work. Unfortunately, the equilibrium route to the cross-coupling coefficients involving either mean-squared displacements or GreenKubo type methods is extremely inefficient. V. Conclusions
Figure 11. Effective diffusivities for nitrogen and oxygen as a function of pore width at the two temperatures studied: t ) 0 °C and t ) 25 °C.
increasing temperature as observed in our simulations. It should be born in mind that the Fickian diffusivity contains a thermodynamic factor which is related to the inverse slope of an adsorption isotherm. This thermodynamic factor will depend on density as well as temperature. It is likely that the thermodynamic factor is greater at t ) 0 °C than at t ) 25 °C, hence the greater diffusivity at lower temperatures. eff eff /DO is shown The ratio of the effective diffusivities DN 2 2 in Figure 12. This ratio is greater than unity for all but one data point and shows a weak variation with pore width until H < 3∆ whereupon a large increase is observed, the increase being greater at the lower temperature. This particular measure of dynamic separation shows both qualitative and quantitatively different behavior with pore width than the equilibrium selectivity. Clearly, a large difference in the diffusivities of oxygen and nitrogen exists at the lowest pore width, yielding 1 order of magnitude difference at the lowest temperature. This lowest pore width is still not small enough for the difference in size
The technique of dual control volume grand canonical molecular dynamics (DCV GCMD) has been used to study the diffusivity of an 80:20 mixture of nitrogen and oxygen in a graphite slit-pore as a function of pore width at two different temperatures. We find evidence in support of combined viscous and diffusive transport through these narrow slit-pores. The viscous flow is weak compared to the diffusive flow and becomes increasingly weaker as the pore width decreases. The fluid velocity profiles show evidence of microscopic slip, but toward the center of the channel they are approximately parabolic, implying that the fluid behaves much as if it were a structureless continuum with angular velocity playing no significant role in the flow. The concentration profiles for each component in the mixture show an approximately linear variation with distance along the pore length, suggesting that the diffusion flows are decoupled. The ratio of nitrogen to oxygen permeabilities is a complex function of pore width and is generally greater than unity implying a selectivity for nitrogen. However, this ratio is quantitatively similar to the static selectivity factor measured in either control volume and is qualitatively similar in its pore width dependence, suggesting that this apparent nitrogen selectivity is a purely thermodynamic phenomenon. In contrast to this, a 1 order of magnitude difference between nitrogen and oxygen effective diffusivity is observed at t ) 0 °C and H ) 2.5∆.
Diffusion of O2-N2 Mixtures in Graphite Pores
Overall, the diffusion of oxygen-nitrogen mixtures through a graphite pore displays a complex dependence upon the pore width. Molecular packing appears to play a very significant role in determining the flow of the mixture. Thermodynamic effects are of more importance in these simulations than sieving effects, and hence, we do not find a greater diffusivity for oxygen compared to nitrogen. The packing effects seen in these simulations are not expected to be as critical at lower pressures. Our model does suffer from several deficiencies: First, experiments are conducted using molecular sieving carbon which does not have a sharply peaked pore size distribution. Consequently, diffusion through a network of interconnected pores of different sizes is expected to be different to diffusion within an individual pore. Second, our simulations ignore any pore entrance effects, considering only internal resistances to mass transport. Third, the
Langmuir, Vol. 15, No. 18, 1999 6059
use of a smooth potential in place of an explicit atom representation of graphite may be a poor approximation. In a future publication we will investigate the diffusion of this same mixture for smaller pore widths and lower pressures. We will also attempt to quantify the Poiseuille flow contribution to the total pore flux and to estimate the cross-coupling coefficients. Acknowledgment. The authors thank Dr. David Nicholson and Dr. Don MacElroy for useful discussions during the course of this work. We are grateful to the National Science Foundation for the support of this work through research Grant No. CTS-9896195; supercomputer time was provided through an NSF NPACI grant (No. NPA 205). LA981465U
