Ind. Eng. Chem. Res. 1999, 38, 723-730
723
Molecular Dynamics Simulations of Diffusion in Mesoporous Glass Neil E. Fernandes and George R. Gavalas* Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125
The effect of gas-solid intrapore potential and surface roughness on diffusion in a single capillary was investigated by molecular dynamics simulations. Calculations were carried out for nitrogen and isobutane under free molecular flow conditions in pores of diameter 4-14 Å at temperatures of 200-800 K. The gases were treated as Lennard-Jones atoms and the pore surface was taken as cylindrical, exerting a 9-3 potential. No energy transfer was considered between the gas and solid, but interaction with the roughened pore wall provided the scattering required for diffusive transport. Two effects of the gas-solid potential were examined in some detail. One is the enhanced intrapore gas concentration which increases the flux, and the other is the bending of the molecular trajectories which decreases the flux. In pores of radius 20 Å, both effects were significant for temperatures as high as 500 K and were enhanced as the temperature decreased. For nitrogen, the two effects partially canceled each other over the temperature range examined, resulting in a temperature dependence similar to that of Knudsen diffusion. For isobutane, the partitioning effect dominated the path curvature effect at temperatures as high as 500 K. 1. Introduction Gas transport in porous solids is encountered in many technologies, including adsorption, catalysis, membrane separations, and gas or petroleum production. The pore size can vary widely between low surface area catalyst supports such as R-Al2O3 with pores of around 1 µm and activated or molecular sieve carbons with pores of 1 nm or less. While some materials such as microfiltration or ultrafiltration media have a unimodal and relatively narrow pore size distribution, other materials such as coal chars have a broad and multimodal distribution spanning a range of a few angstroms to several microns. The mathematical description of gas transport in porous solids has been the subject of numerous research topics and will continue to occupy researchers from several disciplines. In view of the multifaceted nature of the problem, progress has been incremental and different investigators have emphasized different aspects of the problem. One can roughly distinguish between three such aspects: the geometry of the individual pore elements, the interconnection of the individual pores (network effects), and, of concern in this paper, the role of gas-solid interactions. Considerable progress has been made using tools from statistical physics to address network effects,1 while computer molecular dynamics simulations have been used to explore the effect of the gas-solid interactions. Less attention has been devoted to the geometry of individual pore elements, perhaps because it is typically a random function of position and cannot be conveniently characterized. We investigate gas transport in pores of diameter from 10 to 40 nm, i.e., from the upper microporous to the lower mesoporous range. Materials whose structure partly or wholly falls in this range of size include sorbents and catalyst supports such as silica gels, activated carbons and aluminas, and polymeric and ceramic membranes. The pores here are such that their diameters are sufficiently larger than the gas molecules * To whom correspondence should be addressed. Phone: (626)395-4179. Fax: (626)568-8743.
so that (i) transport is not activated, but small enough that (ii) transport is influenced by the intermolecular forces between gas and solid. The analysis will be limited to the case of a sufficiently dilute gas such that the interaction between gas molecules is negligible and hence (iii) only the interactions between gas molecules and the solid need be accounted for. Because the gassolid potential affects a significant part of their cross section, small pores have a higher gas concentration resulting in enhanced fluxes. To concentrate on the effect of the intermolecular geometry in as simple a setting as possible, gas transport is considered in an infinitely long cylindrical capillary with molecular scale roughness. For the low gas density considered here, the molar flow rate in the capillary is given by
N)-
K ∆p DporeπRp2 RT ∆l
(1)
where ∆p is the pressure difference across the capillary and ∆l is the length of the capillary. The length considered is usually much larger than the capillary radius Rp so that the diffusivity Dpore calculated for an infinite capillary is applicable. Rp is the nominal capillary radius on which some roughness will be incorporated. K is a partitioning coefficient which accounts for the difference between internal and external concentrations. In addition to increasing the internal gas concentration, the potential field inside the capillary has significant influence on the diffusion coefficient Dpore, as will be discussed below. The diffusion coefficient can be roughly expressed as
Dpore ∝ λuav
(2)
where uav is the mean molecular speed and λ is the mean step size along the capillary axis. While the speed uav depends on the temperature and molecular mass, the mean step size depends on the capillary diameter and the gas-solid interaction. Frequent reversals of direction reduce the step size, whereas persistence of
10.1021/ie9801150 CCC: $18.00 © 1999 American Chemical Society Published on Web 01/09/1999
724
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
the velocity component along the capillary axis after a collision results in longer steps between reversals. A fundamental calculation of the mean step size must, therefore, consider the momentum and energy exchange between gas molecules and the solid. Numerous experimental and theoretical studies have been published examining gas-solid interactions,2 but these studies have not been applied to gas transport in porous solids. Molecular beam studies on single crystal faces show that the angular distribution of the reflected molecules depends on the angle and the energy of the incoming beam and generally deviates from the cosine law of reflection. Likewise, the degree of thermal equilibration (the accommodation coefficient) also depends on the energy and angle of the incoming molecules. In some theoretical studies, the motion of the gas molecule has been simulated using coupling with a few atoms or ions of the solid which are, in turn, coupled to the vibrational modes of the rest of the lattice.3 Such calculations require drastic approximations and have been performed only for smooth crystal faces. A much simpler approach is taken in the customary derivation of the Knudsen diffusion coefficient,4,5 where the molecules are treated as hard spheres executing straight-line trajectories reflected diffusely from the pore walls (cosine law) and with complete thermal equilibration. Using these assumptions, one can trace the molecular trajectories through complicated pore geometries. For example, Monte Carlo simulations were used to compute the diffusion coefficients in multidimensional bundles of capillaries, in fibrous structures,6-8 and in irregularly shaped pores.9 In the aforementioned calculations using ballistic trajectories, the gas-solid interactions are not taken into account except by postulating that each collision results in a diffuse reflection so that all gas-solid pairs behave identically. A more refined approach has been taken in a series of studies by Nicholson, Petropoulos, and collaborators.10-12 These authors developed a molecular path tracing method (MPT) using a passive intrapore potential to calculate the molecular trajectories but retained the condition of diffuse reflection with the pore wall. The term “passive” here indicates that the gas-solid interaction is not treated as a coupled motion but that the effect of the solid is simply to impose a steady position-dependent potential on the gas molecules. The MPT analysis first utilized a linear potential for a two-dimensional slit-shaped pore10 and was subsequently extended to a cylindrical pore with a more realistic integrated Lennard-Jones (LJ) potential11 and to series and parallel combinations of capillaries.12 These studies showed that the gas-solid potential influences transport not only by increasing the intrapore gas concentration but also by bending the molecular trajectories. Trajectories between collisions with the wall can either lengthen or shorten the mean step, as compared with ballistic trajectories, depending on the molecule velocity after reflection. It was also emphasized that the customary separation of the flux into a surface component and a gas-phase component becomes increasingly ambiguous as the pore diameter decreases.10,11 In the MPT method, the pore wall is treated as a geometrical surface on which hard-sphere collisions with diffuse reflection take place. However, in the case of real (soft) potentials, defining a collision becomes increasingly ambiguous with decreasing pore diameter. Instead of assuming diffuse hard-sphere reflection, one can
employ a soft potential to represent the effect of the solid. Integrated potentials have been used extensively in Monte Carlo calculations of adsorption isotherms, but such smooth potentials do not allow for reversals of the axial velocity which are necessary for the diffusive transport. To allow for velocity reversals, it is necessary that the potential incorporate some form of heterogeneity or roughness. For example, the potential obtained by summation rather than integration over individual atoms in a geometrically perfect crystal possesses a certain degree of heterogeneity. Alternatively, and computationally more simple, one can use a smooth integrated potential with superimposed potentials of discrete clusters, or scatterers, to introduce a mechanism for velocity reversal, with the size and number density of the scattering clusters determining the frequency of the velocity reversals. The purpose of this paper is to explore the use of such composite potentials to describe the effects of path bending, partitioning, and path reversal on gas transport. The trajectories obtained using these composite potentials are averaged to calculate the diffusion coefficient according to Einstein’s equation
2tD ) 〈z2〉
(3)
where z is a displacement along the capillary axis. As already mentioned, the composite potential used in the calculations is passive in that the total energy of the individual molecules is preserved. However, the equilibrium distribution of initial positions and velocities chosen for the sampled molecules is maintained throughout the duration of the simulation providing the correct thermal averaging. A separate Monte Carlo calculation using the same gas-solid potential is used to obtain the second factor K in eq 1. To allow comparison with prior experimental data, the potential parameters are selected to correspond to those of silica glass. 2. Construction of the Pore Surface and Intrapore Potential There are many ways of simulating the solid potential. One approach utilized randomly distributed nonoverlapping LJ carbon atoms in a cylindrical volume and constructed a pore by hollowing out a cylinder through this structure.13 A more refined procedure for simulating a mesopore in silica was developed in a study of SF6 diffusion in porous glass. The pore was constructed similarly to that in ref 13, but the solid substrate in this case was atomistically modeled by using the BornMayer-Huggins two-body potential between the silicon and oxygen atoms. When the system was relaxed from high temperature to low temperature via a computational annealing method, a silica structure was found that agreed quite well with many of the experimentally observed parameters of silica. Determining molecular trajectories through such pore constructs is computationally intensive because it requires knowing the Cartesian coordinates of all of the solid substrate atoms and summing over all of the interactions between these stationary atoms and the diffusing molecules. As mentioned earlier, the potential used in our simulations consists of two contributions. The first term is the potential due to atoms centered at and below a perfect cylindrical surface of radius Rp, and the second term is the sum of potentials of the scattering centers introduced as surface roughness. To construct the first
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 725
term, an approach similar to that of ref 13 is used, drawing a cylindrical surface of radius Rp within a randomly placed collection of oxygen atoms. The potential energy, averaged over thousands of different axial positions, is calculated as a function of radial position by summing over all of the pairwise LJ potential gasoxygen interactions. In this calculation the LorentzBerthelot mixing rule is used to calculate the parameters of interaction between the gas molecule and oxygen atom in the solid from the gas-gas and oxygenoxygen parameters. After ref 14, the interactions with silicon atoms have been neglected because of the low polarizability of silicon and shielding by the oxygen atoms. The Kirkwood-Mu¨ller formula was used to estimate oo/k and σoo as 64 K and 3 Å, respectively.15 For convenience in the molecular dynamics calculations, values of the aforementioned potential were tabulated at radii ri and parametrized by the analytical expression
[(
σgs3
σgs9
)
3. Molecular Trajectory Calculations
2 2 Φ(r) ) nπgsσgs3 + 3 15 (R + r)9 (R - r)9 p p
(
σgs3
(Rp + r)
+ 3
)]
σgs9
(Rp - r)3
(4)
derived by Powl and Everett for the potential between two semi-infinite slabs.16 Here n is the volumetric number density of atoms in the solid and σgs and gs are the LJ parameters for the individual gas-solid atom interactions. These authors also derived the potential within an infinite cylindrical pore but in the form of a slowly converging series inconvenient for calculations. Expression (4) increases to infinity at a pore radius of Rp, and except for very small values of the pore radius, it is very close to the potential in a cylindrical pore. Equation 4 was used as an interpolation formula for the numerically tabulated potential. The LJ parameters, σgs and gs, were used as adjustable parameters to minimize the sum of errors
MSE )
[ ( ) ( )]
∑i exp -
Φ(ri) kT
- exp -
Φi
kT
To complete the pore assembly, the clusters are arranged at a cylindrical surface of radius Rp -3.2 Å, with 1.6 Å being the approximate Si-O bond length. The arrangement algorithm places the clusters in a random, nonoverlapping configuration by allowing the replacement of a cluster a preset number of times (∼1000) if after an attempted placement an intersection with other surface clusters occurs. The scattering center concentration is taken to be equal to the concentration of silanols on the surface of Vycor glass (with vicinal and geminal pairs being counted as one). From previous experimental measurements, 2.3 nm-2 is a reasonable value for the silanol number.17 To avoid correlation effects, pore lengths of 1000 Å were used and periodic boundary conditions were imposed. The interactions between gas molecules and surface scatterers were assumed to be negligible for intermolecular distances beyond 6σgc.
2
(5)
where Φi is the averaged potential calculated by summation at ri and Φ(ri) is the potential at ri calculated from the analytical expression (4). The second term of the intrapore potential comprises the contributions of the scattering centers, simulating the surface roughness. These centers are taken to have the size and chemical composition of regular silica tetrahedra, but their potential is approximated by a spherically symmetrical function of distance from the center of the tetrahedron as follows. Oxygen atoms are placed at the corners of a tetrahedron consistent with the dimensions of a silica tetrahedron, and the potential is averaged over randomly chosen points on the surface of a sphere centered at the tetrahedron origin. The averaged potential is then parametrized with a LJ potential with the parameters σgc and gc estimated by nonlinear regression. Comparisons between the silica cluster and silica cylinder potentials obtained by summation, and the potentials fitted by the LJ analytical expressions were very favorable, especially near the region of the energy well. Differences in the potential well depth were less than 4%.
The molecular trajectories calculated using the potential defined in the previous section can be called soft sphere-soft wall (SSSW) trajectories. In addition to these SSSW simulations, we have also calculated trajectories for hard-sphere gas molecules and hard wall and scattering centers (HSHW). The latter trajectories consist of straight line segments between collisions and can be rapidly constructed by solving a sequence of quadratic equations.18 HSHW simulations were used to investigate surface roughness effects that would be too time-consuming with the SSSW simulations. The hard pore wall was located at a radius Rp, and the scattering centers were located at a radius Rp - 3.2 Å. The gas molecules were treated as hard spheres with diameters equal to the LJ diameters σgg obtained from the literature, while the diameter of the scattering centers σcc was calculated from the Lorentz-Berthelot mixing rule
σgc )
σcc + σgg 2
(6)
where σgc was the same parameter as that estimated in the SSSW simulations. The HSHW trajectories were obtained as straight line segments with either specular or diffuse reflections on the hard surface or the hard scattering centers. With the diffuse reflection condition, the assumptions crucial for the derivation of the Knudsen diffusion coefficient are fulfilled, and so it is expected that the diffusivities obtained from these simulations will be similar to the values calculated from the Knudsen diffusion formula:
2 DK ) Rpuav 3
(7)
This comparison provides a valuable check of the simulation algorithms. The initial velocities were sampled from the equilibrium distribution, while the initial positions were chosen from a uniform distribution reflecting in this case the absence of an intrapore potential. The SSSW trajectories can be calculated using the classical equations of motion
dqi pi dpi ) , ) -∇Φ dt m dt
(8)
where Φ is the intrapore potential, and qi and pi are
726
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 Table 1. Intermolecular Potential Parameters hydrogen nitrogen isobutane
(gg/k)/K
(gc/k)/K
(gs/k)/K
σgg/Å
σgc/Å
σgs/Å
29 80 313
111 338
65 165
2.78 3.75 5.34
4.07 4.53 5.19
3.03 3.53 4.17
tor algorithms. Typically 300-500 atoms were necessary to obtain reasonable statistics, and simulation times were chosen so that each molecule executed between 100 and 200 wall collisions. The HSHW simulations were fast enough to be run on a Pentium 100-MHz processor, whereas the SSSW code was run on a Cray T3D computer comprising 256 DEC Alpha (21064) microprocessors. In both cases run times of approximately 1 h were required. 4. Partition Coefficient Once a capillary has been specified, as described in section 2, the partition coefficient K, appearing in eq 1 is defined by Figure 1. Distribution of velocity components obtained from the SSSW simulation of nitrogen at 200 K for scatterer density ) 0.5 nm-2.
the coordinates and momenta of the molecules, respectively. Under these conditions the total energy of the molecule remains constant. Taking a large number N of molecules with an equilibrium distribution of initial positions and velocities is equivalent to sampling from the microcanonical (NVE) ensemble. Improved statistics are obtained by sampling from the canonical (NVT) ensemble. This is accomplished by modifying the equations of motion according to the Nose´-Hoover formalisn:19
dqi pi dpi ) , ) -∇Φ - ζpi dt m dt Q dζ 2 dt
)
|pi|2
3
∑i 2m - 2NkT
(9)
(10)
It is interesting to note that the Nose´-Hooverian term ζ in eq 9, in conjunction with the scattering clusters, is sufficient to produce an equilibrium velocity distribution even when the initial velocities are nonequilibrium. Figure 1 shows the velocity distributions in the three Cartesian directions obtained from the simulation of nitrogen in a model rough-walled pore. The initial positions were chosen at random, but the initial velocities in each of the Cartesian directions were set at (π/ 8)1/2(RT/M)1/2 for all molecules. Despite this nonequilibrium initial condition, it is clear that MaxwellBoltzmann velocity distributions were obtained as a result of the Nose´-Hooverian energy exchange and the orientational randomness engendered by the scattering clusters. The numerical results presented in this paper were obtained using eqs 9 and 10, but the results of NVE simulations were very similar as long as the initial conditions (positions and velocities) for the trajectories were sampled from the equilibrium distribution. The equations of motion for all SSSW simulations were integrated using a fourth-order Runge-Kutta algorithm with adaptive step size control. This algorithm is slow and is usually shunned by investigators in this field, but in our simulations it was found to be particularly robust when compared to the classical predictor correc-
K)
cpore ) cbulk
∫Ve-Φ/kT dV V
(11)
where V is the volume of the cylinder of radius Rp and a suitably large length L. The integral is approximated by
K=
1
N
∑
Ni)1
exp
( ) -Φ(xi) kT
(12)
where a large number N of points xi are randomly selected inside the cylinder, making no exception for the presence of scattering centers within that volume. A Henry’s law constant KH can also be defined by a similar integral or summation except that if at a randomly selected point the potential energy is larger than 100kT, the point is discarded and a new point chosen. If we write
K ) K HK V
(13)
the factor KV < 1 represents the volume fraction remaining after subtraction of the volume occupied by the scattering clusters. The flux enhancement due to the factor KH is often referred to as “surface” flow. 5. Results and Discussion 5.1. Molecular Interaction Parameters. Using the procedures outlined in section 2, the parameters appearing in Table 1 were calculated. The resulting interaction parameters were found to be very weakly influenced by temperature and only slightly dependent on the pore radius, with deviations in the parameters for the latter case being less than 10% over the range of pore sizes investigated. In light of this weak dependence and the qualitative nature of the simulation, the parameters in Table 1, used for all of the simulations, were calculated at 393 K for pores of radius Rp ) 10 Å for hydrogen, isobutane, and nitrogen. The gas-gas interaction parameters were obtained from ref 20. 5.2. HSHW Simulation Results. The HSHW simulations were used to explore the effects of the surface roughness separately from the effects of the intrapore potential. Figure 2 illustrates the results for simulations
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 727
Figure 2. HSHW diffusion results for isobutane, nitrogen, and hydrogen in pores of Rp ) 25 Å at 393 K.
Figure 3. Intrapore potential for nitrogen in silica pores of various radii. The potential has been normalized with kT at 393 K.
in pores of radius Rp ) 25 Å. The pore diffusivities of hydrogen, nitrogen, and isobutane, normalized with their mean molecular speeds, are plotted as functions of the scatterer surface density. Here, the filled markers represent the results of simulations utilizing specular reflection and the empty markers represent the results for diffuse (or cosine) reflection. From eq 7 the normalized diffusivity for diffuse reflections is expected to be ∼12 × 10-10 m, which is in accord with the simulations. The curves for specular reflection exhibit a minimum near a scatterer density of ∼1.7 nm-2, at which the diffusivities are approximately twice those obtained using diffuse reflection. This minimum is not unexpected in that at both very low and high cluster densities the roughness of the surface diminishes. At low cluster densities molecules mostly collide with the smooth wall, resulting in large effective step sizes in the axial direction. At high densities the clusters begin to overlap, effectively producing a smooth pore wall, resulting once more in a larger step size. It is also observed that the effect of surface roughness depends on the diffusing molecule. At high scatterer densities the normalized diffusivity clearly is an increasing function of molecular size; small molecules find the surface rougher than large molecules. As one might expect, the diffusivities obtained under diffuse reflection conditions are only slightly dependent on the scatterer density. Because the scattering engendered by the cosine law overwhelms the additional scattering caused by the clusters, the value of the diffusivity obtained using specular reflection can be regarded as an upper bound for a given degree of surface roughness. If the intrapore potential field can be ignored, as for example at high temperatures or for weakly adsorbing gases such as hydrogen, one can expect the true pore diffusivity to lie between the results obtained for specular reflection and diffuse reflection. Obtaining more quantitative results would require incorporating into the simulation other scattering mechanisms such as the phonon-molecule coupling and accounting for the shape and internal degrees of freedom of the gas molecules. 5.3. SSSW Simulations. Because of the computational expense of these simulations, only isobutane and nitrogen were studied and only in pores of cluster
Figure 4. Intrapore potential for isobutane in silica pores of various radii. The potential has been normalized with kT at 393 K.
density 2.3 nm-2. From prior experimental results and from their boiling points, it is known that isobutane and nitrogen show strong and moderate adsorbing effects, respectively. The potential energy curves are shown for both gases in Figures 3 and 4. Here the potential energy, averaged over the axial and tangential coordinates, has been normalized against kT at 393 K to give an indication of the adsorbing potential of the pores at moderate temperatures. Because of the ambiguity associated with the pore radius in the presence of the scatterers, an effective pore radius Ri has been defined as
1 Ri ) Rp - 3.2 - σic 2
(14)
where the subscript i refers to the gas molecule, σic is the LJ diameter for the gas-scatterer interaction, and all dimensions are in angstroms. Equation 14 also accounts for the finite sizes of the molecules, which
728
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
Figure 5. Temperature dependence of the partitioning constant K and its factors KH and KV, for nitrogen in silica pores of radius Rp ) 20 Å.
Figure 6. Variation with pore radius of the partitioning constant K and its factors KH and KV, for nitrogen at 393 K.
become more important in small pores. As the effective pore size decreases, the average potential energy becomes more negative, suggesting that partitioning effects will be enhanced. However, it is also noticeable that the well depth ∆E, as defined in Figure 4, also decreases. Figures 5 and 6 respectively demonstrate the effects of temperature and pore radius on the partitioning of nitrogen. As expected the volume exclusion factor KV decreases with decreasing pore radius and is essentially independent of temperature because of the very steep repulsive branch of the potential. The Henry’s law constant KH decreases with increasing temperature and with increasing radius. The overall partition constant K ends up decreasing with temperature but increasing with radius. The effect of temperature on the pore diffusivity is shown in Figure 7 for nitrogen and isobutane in pores of Rp ) 20 Å. The pore diffusivities have been normalized by the average molecular speed to eliminate the
Figure 7. Temperature dependence of the nitrogen and isobutane diffusivities in silica pores of radius Rp ) 20 Å.
Figure 8. Temperature dependence of nitrogen and isobutane normalized permeabilities in silica pores of radius Rp ) 20 Å.
effect of mass and the effect of temperature as prescribed by the Knudsen diffusion expression. Hence, the normalized diffusivity is expected to be constant for molecules strictly behaving with a Knudsen temperature dependence, as indicated by the solid guide line in Figures 7 and 8. This expectation is borne out by the results for nitrogen at high temperatures, but the reduction in the diffusivity with decreasing temperature is dramatic for both gases. The decreased diffusivity is attributed to the path curvature effect which is expected to be a strong function of exp(-∆E). Figure 8 presents the normalized product of diffusion coefficient and partition coefficient versus temperature. The product KD is proportional to the molar flux, which is the experimentally significant quantity. The factor K approximately cancels the path curvature effect in the case of nitrogen and overwhelms it in the case of isobutane. The results for isobutane are in accord with measurements of isobutane permeation in untreated Vycor at 450 K.17 The results for nitrogen are also in accord with the experimental results in (17) but not with previous
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 729
Figure 9. Variation with the pore radius of the nitrogen and isobutane pore diffusivities, normalized with the Knudsen diffusivity, at 393 K.
Figure 10. Variation with the pore radius of the nitrogen and isobutane effective pore diffusivities, normalized with the Knudsen diffusivity, at 393 K.
measurements21 showing “surface flows” of nitrogen in untreated Vycor at temperatures as low as 270 K. This discrepancy may be due to a too weak intrapore potential adopted in the simulations or to geometric effects. The pores of Vycor, although approximately cylindrical, may contain converging and diverging segments involving length scales larger than the size of the scattering centers used in the present simulations. As the intrapore potential becomes more negative with decreasing pore radius, it is reasonable to expect stronger partitioning and path curvature effects. The diffusivities of the two gases normalized with the Knudsen diffusion coefficient (as calculated at the effective pore radius Ri) are shown as a function of pore size in Figure 9. As the pore radius decreases, the normalized diffusivities of both gases increase toward the Knudsen limit. Here again the dramatic effect of path curvature is noticeable. In the HSHW simulations (Figure 2) specular reflection typically gave pore diffusivities that were at least twice the Knudsen diffusivities, but in Figure 9 at pore radii of around 20 Å the pore diffusivity is significantly below the Knudsen value. The decrease of the path curvature effect with decreasing pore radius might have been anticipated after the discussion concerning the size of the potential energy wells in Figures 4 and 5. For isobutane in pores of radius ∼4 Å the potential well is negligible. Figure 10 presents the diffusivities normalized to incorporate the effect of partitioning. Once more it is clear that for isobutane partitioning has the dominant influence on diffusivity. The diffusivity of nitrogen, on the other hand, slightly decreases with decreasing radius because of the volume exclusion factor included in the partitioning constant K. However, with decreasing pore radius and increasing intrapore gas concentration, the assumption of negligible interactions among the gas molecules eventually breaks down and the diffusivity obtained by this simulation procedure ceases to be physically meaningful.
single capillary. By neglect of additional momentum and energy exchange due to coupling of the gas dynamics with solid vibrations, the simulations probably overestimate the diffusion coefficients by as much as a factor of 2. Two important effects of the intrapore potential were explored: partitioning, i.e., enhanced intrapore gas concentration, and path curvature. For isobutane by far the most important effect was the partitioning on account of which diffusivities in excess of the Knudsen values were calculated at temperatures as high as 500 K. For nitrogen the simulations suggest that partitioning and path curvature effects cancel and result in a temperature dependence close to that of Knudsen diffusion over the temperature range of 200-900 K. This result is in some variance with experimental results showing that surface flows of nitrogen are significant at temperatures as high as 270 K.
6. Conclusions The potential between gas molecules and a rough solid surface was employed in molecular dynamics simulations to calculate the diffusion coefficient in a
Acknowledgment This research was supported by DOE University Coal Research Program (Grant DE-FG22-92PC92525). Notation Symbols not presented here are defined in the text. Roman Symbols c ) gas concentration, mol m-3 Dpore ) pore diffusivity, m2 s-1 DKnudsen ) Knudsen diffusivity K ) partitioning constant KH ) Henry’s constant KV ) volume exclusion constant L ) length of the macroscopic sample, m m ) molecular mass, kg N ) molar flow rate, mol s-1 n ) number of attracting centers per unit volume ∆p ) pressure difference across macroscopic porous samples, N m-2 R ) gas constant, J mol-1 K-1 Rp ) pore radius as defined by the 9-3 potential, m r ) radial distance in the pore T ) temperature, K uav ) mean molecular speed, ms-1
730
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
V ) pore volume, m3 z ) axial distance in the pore, m Greek Symbols R ) void fraction ) measure of the depth of the potential well in the pore and also the Lennard-Jones interaction energy, J λ ) mean free path, m σ ) Lennard-Jones molecular diameter, m φ ) potential energy, J kg-1 Subscripts c ) scatterer g ) gas p ) pore s ) surface
Literature Cited (1) Burganos, V. N.; Sotirchos, S. V. Diffusion in Pore Networks; Effective Medium Theory and Smooth Field Approximation. AIChE J. 1987, 33 (10), 1678-1689. (2) Steele, W. A. The Interaction of Gases with Solid Surfaces; Pergamon Press: New York, 1974; Chapter 6. (3) Wang, L.; Ge, Q.; Billing, G. D. Molecular Dynamics Study of Diffusion on a Cu(111) Surface. Surf. Sci. 1994, 301, 353-363. (4) Present, D. R. Kinetic Theory of Gases; McGraw-Hill Book Company: New York, 1958. (5) Jackson, R. Transport in Porous Catalysts; Elsevier Scientific Publishing Company: Amsterdam, The Netherlands, 1977. (6) Burganos, V. N.; Sotirchos, S. V. Knudsen Diffusivities in Parallel Multi-dimensional or Randomly Oriented Capillary Structures. Chem. Eng. Sci. 1989, 44, 2541-2562. (7) Tomadakis, M. M.; Sotirchos, S. V. Effective Knudsen Diffusivities in Structures of Overlapping Fibers. AIChE. J. 1991, 37, 74-86. (8) Tomadakis, M. M.; Sotirchos, S. V. Knudsen Diffusivities and Properties of Structures of Unidirectional Fibers. AIChE. J. 1991, 37, 1175-1186. (9) Burganos, V. N.; Payatakes, A. C. Knudsen Diffusion in Random and Correlated Networks of Constricted Pores. Chem. Eng. Sci. 1992, 47 (6), 1383-1400.
(10) Nicholson, D.; Petrou, J.; Petropoulos, J. H. Calculation of the Surface Flow of a Dilute Gas in Model Pores from First Principles: 1. Calculation of Free Molecule Flow in an Adsorbent Force Field by Two Methods. J. Colloid Interface Sci. 1979, 71 (3), 570-579. (11) Nicholson, D.; Petropoulos, J. H. Calculation of the Surface Flow of a Dilute Gas in Model Pores from First Principles: 2. Molecular Gas Flow in Model Pores as a Function of Gas-Solid Interaction and Pore Shape. J. Colloid Interface Sci. 1981, 83 (2), 420-427. (12) Nicholson, D.; Petropoulos, J. H. Calculation of the Surface Flow of a Dilute Gas in Model Pores from First Principles: 3. Molecular Gas Flow in Single Pores and Simple Model Porous Media. J. Colloid Interface Sci. 1985, 106 (2), 538-546. (13) Bojan, M. J.; Vernov, A. V.; Steele, W. A. Simulation Studies of Adsorption in Rough-Walled Cylindrical Pores. Langmuir 1992, 8, 901-908. (14) Brodka, A.; Zerda, T. W. Molecular Dynamics of SF6 in Porous Silica. J. Chem. Phys. 1991, 95 (5), 3710-3718. (15) Steele, W. A. The Interaction of Gases with Solid Surfaces; Pergamon Press: New York, 1974; p 30. (16) Powl, J. C.; Everett, D. H. Adsorption in Slit-like and Cylindrical Micropores in the Henry’s Law Region. J. Chem. Soc., Faraday Trans. 1 1976, 72 (61), 619-635. (17) Fernandes, N. E.; Gavalas, G. R. Gas Transport in Porous Vycor Glass subjected to Gradual Pore Narrowing. Chem. Eng. Sci. 1998, 53 (5), 1049-1058. (18) Haile, J. M. Molecular Dynamics SimulationsElementary Methods; John Wiley & Sons Inc.: New York, 1992. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1993. (20) Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; John Wiley & Sons Inc.: New York, 1964. (21) Hwang, S. T.; Kammermeyer, K. Surface Diffusion in Microporous Media. Can. J. Chem. Eng. 1966, 44, 82-89.
Received for review February 18, 1998 Revised manuscript received September 4, 1998 Accepted October 19, 1998 IE9801150
