Diffusion within Ultrathin, Dense Nanoporous Silica Films - American

ARTICLE pubs.acs.org/Langmuir

Diffusion within Ultrathin, Dense Nanoporous Silica Films Thomas C. McDermott, Taslima Akter, J. M. Don MacElroy,* Damian A. Mooney, Michael T. P. McCann, and Denis P. Dowling UCD School of Chemical and Bioprocess Engineering, the SEC Strategic Research Cluster and the Centre for Synthesis and Chemical Biology, the Conway Institute, University College Dublin, Belfield, Dublin 4, Ireland ABSTRACT: In this work the origin of permselectivity in dense silica films which possess a pore structure with pore sizes commensurate with the molecular size of the diffusing gas species is investigated. Much of the recently reported work in this field has involved the development of composite membrane films, and while it is generally assumed that the transport process of the gas species within the selective layer of these films is activated in nature, there are anomalies with this simplified picture. In this paper a new model is developed which, for the first time, explains the permselective behavior of the thin selective coatings ubiquitous to membrane separation processes. The model involves the existence of two primary transport domains within the solid film, one of which rapidly conducts the permeating gas (under non-Fickian conditions), while the second domain involves a slow diffusion mode characterized by normal Fickian transport. To validate the model, molecular dynamics simulations are conducted for diffusion of a number of simple gases (He, N2, and CO2) within silica glasses over a range of solid densities. The silica media employed in these studies are based on a novel approach developed in this work for the construction of three-dimensionally periodic atomistic structures of silica of arbitrary density in which network bond connectivity is ensured. The results obtained from this work are in qualitative agreement with experimental observations and confirm the existence of dual mode transport which is central to the interpretation of the permselectivity in composite membranes systems.

permeances reported in refs 1 and 3 within composite films of silica prepared via chemical vapor deposition of alkoxysilanes decreased exponentially with increasing thickness of the deposited film. This is contrary to Fickian predictions in which the permeance is inversely proportional to the membrane thickness. (2) All gases indicated a finite, nonzero permeance for the thickest SiO2 films deposited. In these studies and in other work reported earlier by Gavalas and co-workers6
9 the primary goal has been to develop the fabrication protocols for composite membranes for use in the separation of H2 from high temperature process gas streams, and it has been demonstrated that high permselectivities with moderately high permeances for H2 can be achieved in such systems. However, the results reported in refs 1 and 3 clearly suggest that if optimally designed membranes are to be fabricated for gas mixture separation for more general applications (and particularly for “difficult” mixtures such as CO2, O2, and N2 in high temperature postcombustion or oxy-fuel carbon capture technologies) then ultrathin permselective coatings with thicknesses less than 25
50 nm will need to be produced. This will require an understanding of the underlying mechanism for the observed exponential dependence of the permeance on membrane

1. INTRODUCTION A major field of ongoing study in gas mixture separation technology involves the development of novel composite ceramic membranes suitable for use across a broad range of process conditions. A common theme in these studies is the use of support membranes composed of a mesoporous/microporous substrate materials (e.g., α-alumina/γ-alumina, Vycor) with a minimum pore size of approximately 5 nm, in some cases supplemented with a further metal or silicon oxide microporous layer with pore dimensions of 1
3 nm, upon which, as illustrated in Figure 1, is deposited a submicrometer thick permselective layer of silica containing pores in the angstrom range (see for example refs 1
4). The latter layer is the selective coating which itself is typically 30
100 nm thick. In order to more fully understand the transport mechanisms within these composite materials and, in particular, how one may develop optimum conditions for both maximum permselectivity and maximum permeation rate, exploratory simulation studies are reported in this work in which equilibrium molecular dynamics simulations are conducted to evaluate the diffusive transport characteristics of simple gases (He, N2, and CO2) within model dense silica media over a range of temperatures. These studies are guided by prior work reported in ref 5 and, most importantly, by two observations which may be drawn from recent work on permeation within composite thin silica membranes:1,3 (1) For a range of simple gases (He, H2, CO, CO2 and CH4) the r 2011 American Chemical Society

Received: October 12, 2011 Revised: November 23, 2011 Published: November 30, 2011 506

dx.doi.org/10.1021/la203994v | Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

0.56 of the bulk silica density demonstrated that the permeance of the gases are well described by the simple expression h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i Pi ¼ χi 1=2πMi RT expðδi =l i Þ expð
L=l i Þ

where χi is the nanopore accessibility coefficient for the distribution of molecules of species i between the bulk gas and the membrane surface (this is distinct from the adsorption equilibrium coefficient) and Mi is the molecular weight of the gas. The parameter δi represents the length scale of the interfacial inhomogeneities which exist at the surface of the membrane.5 Equation 2 can provide a partial explanation for the observations reported in refs 1 and 3 but there are issues which need to be resolved. In particular, it is known that the permeance of the membranes reported in refs 1 and 3 for molecular species such as CO2 and N2 are, by extrapolation, nonzero at infinite thickness, and therefore one must conclude that eq 2 is incomplete. In addition, there were two major assumptions built into the silica models investigated in earlier work:2,5,10 (1) The porous silica media were generated by randomly etching bulk amorphous silica of density 2.2 g/cm3. During the etching process care was exercised vis-a-vis the creation of surface silanols; however, the resulting porous medium was otherwise random. No criterion to ensure bond connectivity within the bonded network of silica tetrahedra was included and it is known that a random structure of this type with isolated clusters of silica can display percolation thresholds for diffusion which are at significantly lower densities than might be expected.12 (2) To simplify the computation, the silica models employed in this earlier work were static. This is a condition which cannot be assumed to hold for small molecule transport within pores which are very close to the size of the gas molecule itself. In this paper, the above assumptions are addressed and molecular dynamics simulations have been performed to recompute the MSDs for three gas species (He, N2, and CO2) diffusing within silica media composed of bonded networks of fully mobile atomistic silica.

Figure 1. (Left) SEM image of a thin film of SiO2 deposited onto a substrate composed of a secondary 500 nm SiO2 layer (nominal pore size 1 nm), a microporous γ-alumina (nominal pore size 5
7 nm), and a mesoporous α-alumina (nominal pore size in the micrometer range).4 The outermost layer on the right of the SEM image is the primary permselective SiO2 film which is 100 nm in thickness with a nominal pore size of 0.4 nm. (Right) A schematic drawing of the principal elements of a permselective membrane.

thickness as mentioned above in point 1. For example, it is noted that the CO2, CO, and CH4 permeances reported in ref 1 decreased exponentially by 5
6 orders of magnitude as the thickness of the deposited nanomembrane increased from 0 to 30 nm. Very similar results were reported in ref 3. In prior studies2,5,10 of diffusion and mobility within moderately dense silica media the permselective characteristics of CO 2, O2 , and N2 were investigated. In this earlier work it was shown that if the density of the silica medium is high enough then the material becomes nonpercolating for a silica system which is of infinite extent. Here, the molecular dynamics simulation of fully atomistic models of the gas molecules diffusing through the interstices of static amorphous silica structures (also atomistically modeled) was employed, and it was shown that the mean-square displacements of the molecular centers-of-mass of the trace gas species as a function of time, t, are well represented by a scaling theory approximation5,11

2. SIMULATION METHODOLOGY 2.1. Moderate-to-High Density Silica Structures. Starting with a periodically imaged, amorphized silica structure13,14 (containing stoichiometric SiO2) of density 2.2 g/cm3 and dimensions 3.573 nm3, a range of moderate-to-high density nanoporous silica media have been created to serve as model structures for gas transport simulations. The primary issue involved in the creation of these atomistic structures is the problem of maintaining the bond connectivity of the structure within the fundamental cell and across the system periodic boundaries as well as identifying “floating” components of the structure (groups of atoms which have become detached from the main structure). This is solved here using a nonlattice version of the Hoshen-Kopelman cluster counting algorithm.15,16 The creation of lower density structures is carried out in two parts: (a) removal of atoms from the amorphous vitreous SiO2 structure until a specified density is reached (following the etching procedure described in earlier work5 with the exception of further details to be described below), while ensuring that the system remains fully coordinated and fully connected (no floating clusters); (b) checking that if the structure is periodically replicated in the x, y, and z directions there is only one connected molecular network.

0

ÆjΔrðtÞj2 æ ∼

At 2=dw 0 1 þ ðA=6l 2 Þt 2=dw

ð2Þ

ð1Þ

The coefficient d0 w is the fractal dimension of the random walk undertaken by a given particle, the parameter l characterizes the maximum range of the particle displacements under nonpercolating conditions (here referred to as the mean percolating cluster size) and the parameter A is an anomalous scaling term which is a constant for a given species at a given temperature and density and is treated as a simple fitting coefficient. The results obtained demonstrated the existence of percolation thresholds for each of the three gases CO2, O2, and N2 in the density region 0.52 < FSiO2/FBulkSiO2 < 0.56 where FBulkSiO2 = 2.2 g/cm3. A more detailed nonequilibrium molecular dynamics analysis5 of the steady state permeation characteristics of gases diffusing through membranes of finite thickness and densities greater than 507

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

Table 1. Number Fractions of Silicon, Siloxane (Doubly Bonded) Oxygens, Silanol Groups, and the Average Pore Radii of the Silica Structures Employed in the MD Simulations FSiOx (g/cm3)

xSi

xO

xOH

Rp (nm)

1.73

0.307

0.532

0.161

0.128

1.86

0.314

0.571

0.115

0.115

1.97

0.320

0.600

0.0801

0.105

2.03

0.325

0.625

0.0495

0.0997

2.2

0.333

0.667

0.0

0.0935

The full version of the nonlattice Hoshen-Kopelman algorithm is not required in step (a) to check the connectivity of the molecular network—only the nodes (corresponding to the atoms in the present case) in the original algorithm are used with the links omitted and use is made of a simplified form of the algorithm. Hence the labeling algorithm required for the etching process only requires the atom connectivity table (and the number of atoms) as input. To conduct step (b) the list of bonds which cross the periodic boundaries of the system are easily obtained by considering if the components of the separation vector between two adjacent atoms rij are greater than half the simulation box length L. In order to determine if a molecular cluster spans the simulation box in a particular dimension, the periodic boundary conditions in that dimension can be switched off (keeping the periodic boundary condition in the other dimensions) with a subsequent check if the atoms in the boundary crossing bond still have the same cluster number. To check structures of lower, but not very low, density (down to 70% say) it is usually not necessary to carry out a detailed analysis to check the spanning in each dimension (turning off one periodic boundary condition at a time), and instead, all the periodic boundary conditions may be switched off and one can check for a cluster that still spans all three dimensions. The above algorithm has been applied to generate silica structures with densities 1.73 g/cm3, 1.86 g/cm3, 1.97 g/cm3, and 2.03 g/cm3. Initial etching runs were conducted (with full bond connectivity) to determine the appropriate level of silicon atom removal, with subsequent condensation of a prescribed set of dangling |
Si
O bonds (modeled as silanols), which would lead to the desired target density. In all cases the prescription applied for the level of “silanol” condensation was that all dangling oxygens created as a result of the raw etching process which were within 0.28 nm of each other (which is fractionally larger than the oxygen nearest neighbor distance in bulk silica) were “condensed” with the removal of one of the oxygens in every given pair. For each pair of dangling bonds condensed, 105 Metropolis Monte Carlo relaxation steps were undertaken for the silica system locally within a range of 1 nm of the site of a given condensation reaction. Selected average properties for each of these structures, determined from samples of 25 configurations each for the four system densities, are provided in Table 1. The pore size distributions within these media (see Figure 2) and the average pore radii were determined using the technique described by Bhattacharya and Gubbins.16 2.2. Equilibrium Molecular Dynamics Method and Potential Functions. For each of the 25 silica configurations at the four different densities summarized in Table 1, 50 gas particles were inserted at random locations to achieve reasonable statistical sampling of the starting points of the individual trajectories. The

Figure 2. Normalized probability density distributions of pore sizes within moderate-to-high density nanoporous silica media: 9, 1.73 g/cm3; 0, 1.86 g/cm3; 2, 1.97 g/cm3; Δ, 2.03 g/cm3.

integration technique for the molecular dynamics method employed in this work for the fully atomistic structures was the velocity Verlet algorithm18 with a time step of 0.5 fs. An initial period (100 ps) for equilibration was employed with velocity rescaling and thereafter the trajectories evolved within the microcanonical (NVE) ensemble. All of the trajectories were of 2.5 ns duration and the mean-squared displacements of the centers-of-mass of the particles were evaluated over 0.5 ns at time intervals of 0.05 ps with equal weight given to each of the 104 data points. In the computations the gas species are treated as ideal gases, i.e., there are no explicit interactions between the gas molecules within any given silica configuration, and the individual atom
atom pair interactions are considered to be composed of both intermolecular nonbonded terms and intramolecular contributions. Note that, to ensure no distortions arise locally within the solid, the initial insertions of the gas particles are conducted in such a way that no two gas particles are within a molecular diameter of one another. The nonbonded interactions are described by Lennard-Jones ULJ(rij) and Coulombic Uel(rij) (partial charge) potential terms Uðrij Þ ¼ ULJ ðrij Þ þ Uel ðrij Þ The Lennard-Jones interaction is 0 !12 !6 1 σ σ ij ij A ULJ ðrij Þ ¼ 4εij @
rij rij

ð3Þ

ð4Þ

where rij is the distance between atoms i and j, εij represents the magnitude of the minimum in well-depth, and the distance σij characterizes the separation between the two atoms at which repulsion dominates the interaction. A cutoff radius for the Lennard-Jones interaction of 1.2 nm was employed in the simulations. In the evaluation of the interactions between the partial charges on the gas molecules and the charges within the silica system, the Coulombic potential in its screened form19,20 is employed in this work. In these computations the following simplified form which is known to provide reliable estimates for 508

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

the electrostatic interactions21 was utilized ! ( ) qi qj 1 1 1 Uel ðrij Þ ¼
þ ðrij
Rc Þ Rc2 4πε0 rij Rc

between bonds i and j, and θij0 is the corresponding equilibrium value. kib and kijθ are the bond stretching and bond bending force constants, respectively. The parameters of the Keating potential for the silica system and for the intramolecular modes of N2 and CO2 (also evaluated using the functional form in eq 7) are provided below in Table 3. In the case of CO2 the mathematical form of the bond angle bending potential for bonds i and j appearing in eq 7 is replaced by27

ðrij e Rc Þ ð5Þ

where rij is the distance between two partial atomic charges qi and qj on particles i and j, respectively, ε0 is the electrical permittivity of space and the cutoff Rc was set equal to 1 nm. The Lennard-Jones and electrostatic (partial charge) parameters on each atom used in the model systems are reported in Table 2. In this work the OH groups are treated as united atom structures with the charge on the OH group defined by combining the individual charges on the oxygen and hydrogen atoms. To assign the charge on the surface silicon atoms, qSi, a simple formula has been employed which is given by qSi ¼
qOH NOH
qOðSiloxaneÞ NOðSiloxaneÞ

Uij, Bend ¼

3. MODELING The underlying premise of the analysis leading to eqs 1 and 2 is believed to be fundamentally correct with one exception: in order to satisfy observations that normal Fickian diffusion (no matter how slow) of simple gas species can take place within silica films of the type shown in Figure 1 an additional term (or terms) must be included within the analysis. A model which extends the prior analysis leading to eq 2 is described in the following. (We note that in view of the assumption of gas phase ideality as stated at the beginning of Section 2.2, the model developed in this work and the simulation results reported below apply strictly to transport and permselectivity in the free particle Henry’s Law limit for adsorption/ partitioning. The more complex situation of competitive adsorption and permeation which would prevail, in particular, at low-to-moderate temperatures requires development in future work.) The structure depicted in Figure 3 contains two primary domains, the gray area which represents a poorly conducting region and the white “channels” which have a conductance much larger than the gray domain. As described by Bunde and Havlin28 random walkers within the white “channels” very rapidly sample the configurational space of any given cluster within this domain. Their transfer to the poorly conducting domain is extremely slow and will only occur during a time scale of the order of tX or greater as indicated in Figure 3b. At short times the walker moves on the finite (white) clusters and the MSD is considered to be described be eq 1 during this time range (note that at very short times the MSD will be proportional to t2 in the ballistic regime). For tξ < t < tX The MSD saturates at the value 6l 2 in three dimensions. At intermediate-to-long times the random walker transfers into the gray domain and slowly diffuses through the medium. A simple extension of eq 1 is

ð6Þ

The intramolecular bonded interactions in the silica media are computed using the Keating Si
O bond stretching and O
Si
O, Si
O
Si bond angle bending interaction potential.25,26 This potential represents the energy of the solid as a function of the nearest neighbor positions and is given (in its simplified form25) by 2 1 b 1 θ ki ½ri
ri0 2 þ kij cos θij
cos θij0 E¼ i ∈ bonds 2 ij ∈ angles 2

∑

∑

ð7Þ In eq 7, ri is the distance between two bonded atoms, ri0 is the (average) equilibrium length of bond i, θij represents the bond angle Table 2. Lennard-Jones and Coulombic Partial Charge Parameters (from Schumacher et al.22,23 and Chakravarty24)a ε/kB (K)

σ (nm)

O(silica) Si

185.0 0.0

0.2708 0.0


0.64025 +1.2805

O(OH)

185.0

0.3


0.533

H(OH)

0.0

0.0

+0.206

Atom

Charge, qi(eo)

He

10.22

0.228

N(qb)

34.897

0.33211


0.5475 (+1.095)b

0.0

O (CO2)

82.997

0.3064


0.33225

C

29.999

0.2785

+0.6645

1 θ k ðθij
θij0 Þ2 2 ij

a

See also MacElroy and Raghavan.13,14 b The charge +1.095 is the charge on the center-of-mass required to reproduce the quadrupole moment on N2

0

ÆjΔrðtÞj2 æ ¼

At 2=dw 0 þ 6Dt 1 þ ðA=6l 2 Þt 2=dw

ð8Þ

Table 3. Bonded Intramolecular Potential Parametersa Atomic interaction Si
O

Bond stretching constant (MJ/mol/nm2)

Bond bending constant (kJ/mol)

260.6

Equilibrium distance (nm) 0.161

Si
O
Si O
Si
O

θSi
O
Si,0 = 144 θO
Si
O,0 = 109.471

193 417

N
N

292.9

0.10464

CdO

585.8

0.1161

OdCdO Si
OH

Angle (o)

θOdCdO,0 = 180

418 292.9

0.1487

a

The pure silica parameters (bond length Si
O, and bond angles Si
O
Si, O
Si
O) are taken from von Alfthan et al.26 The bond stretching and bond bending constants for N2, CO2, and Si
OH are from Mayo et al.27 and the N
N, CdO, and Si
OH equilibrium bond lengths are from Schumacher et al.22,23 509

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

Figure 3. (a) “Superconducting” clusters (white)/poor conductor (gray) mixture. (b) Schematic drawing of the MSD of a random walker diffusing within the medium illustrated in (a).

where D is the long-term diffusivity of the gas species within the poorly conducting domain. Noting that the interrelationship between the autocorrelation function of the center-of-mass velocity (VACF) of the gas molecules and the center-of-mass MSD is given by 1 dÆjΔrðtÞj2 æ ¼ 2 dt

Z t 0

jðτÞ dτ

Figure 4. Mean square displacement of helium atoms within moderateto-high density nanoporous silica media: 9, 1.73 g/cm3; 0, 1.86 g/cm3; 2, 1.97 g/cm3; Δ, 2.03 g/cm3. (a) 473 K and (b) 873 K. (Note that, for clarity, only 12 points, equally spaced in time, of the total 104 points for the entire MSDs (solid lines) are shown as symbols in each case).

ð9Þ

one may then consider two contributions to the VACF representing both of the additive terms appearing in eq 8, i.e. jðtÞ ¼ jðξÞ ðtÞ þ jðXÞ ðtÞ

within the downstream chamber of the membrane separation unit (see Figure 1). Omitting the influence of the mesoporous/ microporous support substrate (this can be readily included as needed to refine the analysis), the downstream concentration c(B)(L) is simply related to the upstream selective barrier layer concentration c(0) (and hence c(B)(I)) through an exponential dependence on the barrier layer thickness L5,11 and for conditions in which L > l the net “superconducting” flux across the membrane is given by

For one-dimensional diffusion through membranes of the kind illustrated in Figure 1, one has Z Jðz, tÞ ¼

t


0

Z

t

¼
0

jðt
τÞ

∂cðz, τÞ dτ ∂z

∂cðz, τÞ jðξÞ ðt
τÞ dτ
∂z

Z

t 0

∂cðz, τÞ jðXÞ ðt
τÞ dτ ∂z

ð10Þ

h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i   J ðξÞ ¼ χ RT=2πM expðδ=l Þ expð
L=l Þ cðBÞ ðIÞ
cðBÞ ðIIÞ

The first term on the right corresponds to the “superconducting” flux and, following on from the analysis reported in refs 5 and 11 under stationary conditions, this term is simply given by net effusive flux from the downstream side of the membrane  1  ð11Þ J ðξÞ ¼ v̅ cðBÞ ðLÞ
cðBÞ ðIIÞ 4

ð12aÞ In the stationary limit the second term on the right of eq 10 simply provides the Fickian term  DK  ðBÞ c ðIÞ
cðBÞ ðIIÞ ð12bÞ J ðXÞ ¼ L

where v is the molecular mean speed of the gas molecules, c(B)(L) is bulk gas concentration on the downstream interface of the membrane and c(B)(II) is the bulk gas concentration

where K is the equilibrium adsorption coefficient of the gas ((moles of adsorbed gas/m3)/(moles of bulk gas/m3)). Combining gives 510

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

Figure 5. Mean square displacement of nitrogen molecules within moderate-to-high density nanoporous silica media at 473 K (a) and 873 K (b): 9, 1.73 g/cm3; 0, 1.86 g/cm3; 2, 1.97 g/cm3; Δ, 2.03 g/cm3. The full MSDs (104 points) are shown as the long dashed curves and the nonlinear fits using eq 8 are shown as solid lines in each case.

Figure 6. As in Figure 5 but for CO2.

4. RESULTS AND DISCUSSION The results for the MSDs for each of the gases He, N2, and CO2 at the lowest (473 K) and highest (873 K) temperatures simulated in this work are reported in Figures 4 to 6. These MSDs were computed from 2.5 ns molecular trajectories and a nonlinear regression algorithm was employed to determine the four parameters, A, l , d0 w, and D appearing in eq 8. The nonlinear regression was carried out using a version of the SOLVOPT program29 which implements the Shor algorithm.30 Versions of this program are available;31 however, the program used here was a Mathematica implementation developed in-house. This allowed the results of the nonlinear fitting to be validated by comparison with the built-in nonlinear solvers used in the Mathematica function NonlinearModelFit. It was found necessary to use the numerical minimization NMinimize in order to obtain converged results; otherwise, the Mathematica results were similar to those obtained from SOLVOPT. However, the results from SOLVOPT usually produced the tightest fits and were more robust with respect to the changes in the parameters caused by eliminating points from the initial part of the curve from the region to be fitted. In all of the work reported below (with the exception of the analysis of the helium

the general expression for a single ideal gas species J ¼ PΔp where Δp is the pressure difference across the membrane and P is the permeance " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi #   1 L DK expðδ=l Þ exp
þ P¼ χ 2πMRT l RTL

ð13Þ

In Section 4 below the expression presented in eq 8 is used to correlate MSD data generated as outlined earlier based on EMD simulation of the diffusion of the three gases He, N2, and CO2 in silica media at bulk densities of 1.73 g/cm3, 1.86 g/cm3, 1.97 g/cm3, and 2.03 g/cm3. In this section of the paper attention is solely focused on eq 8. Equation 13 and its implications for the analysis and assessment of permselectivity of silica membranes will be considered briefly in the concluding paragraphs of Section 4. 511

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

Figure 8. Diffusivity of helium as a function of temperature within silica nanoporous structures at densities of 9, 1.73 g/cm3; 0, 1.86 g/cm3; 2, 1.97 g/cm3; Δ, 2.03 g/cm3.

to be less than (20% in all cases based on averages over five bins of five configurations each). In the analysis of the helium MSD results it was very soon clear that the first term in eq 8 was only of marginal importance due to the significantly larger diffusivity for the gas in this case. The log
log plots for the He MSDs shown in Figure 4 clearly indicate anomalous behavior in most cases up to approximately 50 ps; however, for most of the sampled time range the transport process is Fickian. This behavior is not unexpected in view of the small kinetic diameter of the helium atom. To more clearly demonstrate the comparatively small effect the anomalous contribution plays in the case of helium, results for this contribution, obtained by subtracting the Fickian contribution from the overall MSD in Figure 4, are reported in Figure 7 at the two highest densities. These results also demonstrate the scatter associated with this term for Helium and hence the difficulty in estimating the parameters A and l for this gas. The results provided in Figure 7 do, however, demonstrate that the plateau values, 6l 2, for the anomalous component of the MSD are such that at these densities the percolation cluster sizes within the “superconducting” domains for helium are ∼0.29 nm at 2.03 g/cm3 and ∼0.33 nm at 1.97 g/cm3. As will be shown below these are consistent with the results computed for N2 and CO2. The diffusivities for helium obtained from the long-time slopes indicated in Figure 4 (and those √ not shown for 673 K) are plotted in Figure 8 in the form D/ T where it is observed that at all densities the diffusivity displays activated behavior, albeit comparatively weak at the lower densities. A comparison of the model results (eq 8) for both N2 and CO2 are shown in Figures 5 and 6. In all cases the quality of the fit is excellent and the deconvolution of the anomalous and Fickian terms, as shown for example for CO2 in Figure 9 under two extreme conditions, demonstrates that both terms are of similar magnitude over much of the sampling time. The diffusivities computed for N2 and CO2 are plotted as functions of density and temperature in Figures 10 and 11, respectively. At all temperatures D would appear to drop rapidly as the density increases above 2.0 g/cm3 indicating significant confinement effects in the quasi-static silica medium. At 873 K, the highest temperature investigated here and hence the temperature at which one would anticipate the least resistance to

Figure 7. Anomalous component of the mean square displacement of helium atoms within nanoporous silica media at (a) 1.97 g/cm3 and (b) 2.03 g/cm3: - - -, 473 K;


, 673 K; —, 873 K.

MSDs as discussed below) the nonlinear fit was conducted for all data within the time range 10
500 ps. Initial studies employing this regression analysis for all four parameters resulted in significant scatter in the parameters themselves. For this reason the regression procedure was limited in the following respect. The fractal dimension dw0 for continuum percolation is known to be dependent on the nature of the distribution of bond conductivities within the three-dimensional system.28,32 If this distribution is assumed to be uniform, then d0 w is bounded by 3.479 < d0 w < 3.768.32 From the full four parameter regression on the N2 and CO2 data, the parameter d0 w was estimated, on average, to be d0 w ∼ 3.98 ( 0.8 which overlaps with the bounds of the values cited above, particularly the upper bound. On this basis and in view of the scatter observed in the full nonlinear analysis, it was deemed best to fix the value of d0 w at 3.768. The results reported in this work are based on this assumption with the three coefficients A, l , and D in eq 8 now serving as the regression parameters for the fit to all of the MSD data. For the time range 10
500 ps, this approach provided reliable results for both N2 and CO2 as illustrated in Figures 5 and 6 (we note that the standard errors on each of these parameters were determined 512

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

Figure 9. Examples of the deconvolution of the MSD based on eq 8 for CO2 at (a) T = 473 K and FSiOz = 1.73 g/cm3 and (b) T = 873 K and FSiO2 = 2.03 g/cm3. The long dash is the Fickian contribution, the short dash is the anomalous term (the first on the right of eq 8), and the solid line is the combination of both these terms. The filled circles are selected equispaced original MD data.

Figure 10. Fickian diffusivity as a function of density: 9, 473 K; 0, 673 K; b, 873 K. (a) N2 and (b) CO2.

significant density dependence of this parameter (see Figure 12) would suggest that at densities higher than those investigated here the cavities within which N2 or CO2 molecules are trapped are very tight with the cluster pathways for wormlike movement of the N2 and CO2 molecules rapidly approaching the cavity size itself. Transport under these conditions would involve diffusive jumps of the gas particles determined by rare events resulting from the formation of apertures between trapping sites in the solid. The formation of apertures sufficiently wide for transport to take place would require the transformation of significant levels of thermal energy into potential energy, and it is anticipated that the slopes of the lines shown in Figure 11 would increase significantly as the density approaches 2.2 g/cm3. The density dependence of the percolation cluster size, l , as shown in Figure 12a is also of particular interest when compared with the earlier results for static media.2 As was observed for static media there is a distinct discontinuity in the semilogarithmic plot, which is believed to be a signature of the percolation transition for the superconducting channels depicted in Figure 3a. As the “superconducting” channels become fully percolating at lower densities, the percolation cluster size, l , extends to infinity. For the fully mobile, connected silica network media examined here, the influence of the growth in the percolation length scales for

diffusion due to the greater frequency of appearance of low energy pathways (through “free volume” mobility) arising from the intrinsic vibrations in the silica structure, results reported by Kajihara et al. 33 for N2 in vitreous silica provide D ∼ 6.6  10
16 m2/s, which is 4 orders of magnitude smaller than the lower bound of Figure 10a. This suggests that as the density increases toward 2.2 g/cm 3 the trend lines indicated in Figure 10 are indeed correct with a very sharp drop in diffusivity in all cases as the density approaches 2.2 g/cm3 . This may also be inferred from the temperature dependence indicated in Figure 11. At low densities the activation energies for diffusion are comparatively high, reflecting a behavior that might be more akin to diffusion in relatively “wide” pores (see Table 1 and the pore size distributions provided in Figure 2) with hopping between surface sites. At the highest density, however, the gas molecules are tightly confined and the diffusion depends much more on the intramolecular dynamics of the silica tetrahedra and their vibrations relative to neighboring tetrahedra. The rather weak dependence of l on temperature in contrast to the 513

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

Figure 12. (a) Percolation cluster size, l , as a function of density: 9 CO2 and b N2 at 873 K; Δ CO2 and 2 N2 at 673 K; ( CO2 and ) N2 at 473 K. (b) Percolation cluster size as a function of temperature: 9 CO2 and 2 N2 for 1.97 g/cm3; 0 CO2 and Δ N2 for 2.03 g/cm3.

Figure 11. Fickian diffusivity as a function of temperature: 9, 1.73 g/ cm3; 0, 1.86 g/cm3; 2, 1.97 g/cm3; Δ, 2.03 g/cm3. (a) N2 and (b) CO2.

CO2 is clearly observed in the range 79
84% of the density of bulk silica in contrast to the much lower density of 56% of the bulk value for the static, random systems reported by Cuffe et al.2 For N2 the percolation transition occurs below 80% in the present case in comparison with the upper limit of 52% for the structures reported in ref 2. These observations were to be anticipated based on the prior work on random vs connected overlapping and nonoverlapping spheres models of microporous media12 in conjunction with the greater freedom of movement of the gas molecules within the “superconducting” channels in the mobile silica structures as implied by the temperature dependence of the percolation cluster size illustrated in Figure 12b. Before concluding, an approximate analysis for permeation across hypothetical silica films is conducted using eq 13. The five parameters required to estimate the permeance for a given gas species using this equation are the Henry’s law adsorption equilibrium coefficient, the diffusivity of the species in the infinite medium, the percolation cluster size, l , the inhomogeneity factor, δ, and the nanopore accessibility coefficient, χ. Based on nonequilibrium molecular dynamics simulations reported by Moloney34 it may be inferred that as a reasonable approximation the group of terms χ exp(δ/l ) is ∼1.0 for each gas species. The Henry’s law K values for

each of the gases in the individual etched structures are computed from the canonical average of the Boltzmann factor as estimated via the Widom insertion algorithm35 * !+ 0 1Z Uðr Þ 0 Ki ¼ exp
dr ð14Þ V V kB T Ω

where V is the volume of the fundamental cubic cell and the angular bracket Æ...æΩ refers to the canonical average over orientations of the inserted molecule, represented by the solid angle Ω. A total of 108 trial insertions were performed for each of the three etched media at densities 1.73, 1.86, and 1.97 g/cm3 with canonical averaging over orientations of the inserted species and averaging of the 25 nanoporous silica configurations in each case. The results obtained at 873 K were as follows: 1:73 g=cm3 : KðN2 Þ ¼ 0:22

and

KðCO2 Þ ¼ 0:67

1:86 g=cm3 : KðN2 Þ ¼ 0:15

and

KðCO2 Þ ¼ 0:43

1:97 g=cm3 : KðN2 Þ ¼ 0:097 514

and

KðCO2 Þ ¼ 0:26

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

largely determined by the DK/RTL term in eq 13, the permselectivity α for the CO2/N2 systems can be very large for thin films under certain circumstances. This raises the possibility of fabricating composite membranes with both high permeance for one of the major species and high selectivity. The results for helium shown in Figure 13a at 1.97 g/cm3 also clearly demonstrate the influence of the anomalous first term in eq 13. This term is less dominant for the small helium atom, and the results at this density exhibit the very large selectivities (∼1000) frequently associated with the lightest gases such as helium for transport in composite silica membranes relative to other molecular species such as CO2 and N2. Central to future work on this topic will be the determination of the process parameters which will ultimately enable precise control of the parameter, l , which itself largely determines both the magnitude of the peak in Figures 13b and peak position relative to the membrane thickness. Density is clearly an important control variable; however, as shown in a recent study of the deposition of alkoxysilane precursors via chemical vapor deposition36 care should also be undertaken to ascertain the influence of the precursor, CVD reactor temperature, and precursor partial pressure. Process conditions clearly influence the magnitude of the percolation cluster size, l , as may be inferred from the results reported by Araki et al.3 Another element of control relates to the Fickian term involving D. Lowering the diffusivity significantly enhances the relative importance of the non-Fickian contribution and this may be achieved by controlling film density and/or by introducing varying levels of other ionic species into the silica film (see for example Altemose37). Ultimately, however, care will need to be taken (a) to ensure that the permeance is increased as much as possible to ensure separation capacity while achieving high selectivity and (b) to assess the stability of the ultrathin layer in the typically harsh environments in which oxide membranes of these kinds will find preferred applications. Figure 13. (a) Predicted permeances as functions of membrane thickness for CO2 (open symbols) and N2 (filled symbols) for etched model silica densities of 1.73 g/cm3 (squares), 1.86 g/cm3 (circles), and 1.97 g/ cm3 (triangles); the solid line without symbols is the prediction for helium at 1.97 g/cm3. (b) The permselectivities for CO2/N2 computed from the results shown in (a) for 1.73 (9), 1.86 (0), and 1.97 g/cm3 (b).

5. CONCLUSIONS In this work, the diffusion processes of the light gases He, N2, and CO2 within dense amorphous silica systems have been investigated theoretically and a novel approach has been proposed which describes the permeation characteristics of fluids through membranes in a unique way. The model allows one to determine key parameters from a simulation of homogeneous silica systems which may then be adapted to apply to a relatively complete description of the permeance of fluids in the barrier layers of composite membrane systems. The core of the theory is the superposition of two distinct modes of transport within disordered media: (a) anomalous diffusion which exists over time scales less than a few hundred picoseconds in length and which takes place on finite ‘superconducting’ clusters of channels within the disordered structure; (b) normal Fickian diffusion which takes place within the underlying continuous background structure of the amorphous medium. The existence of both these modes of transport is demonstrated by a deconvolution of the mean-square displacement of the center-of-mass of the diffusing particles which ultimately results in a novel expression for the permeance of gases within the disordered system. The theory can explain experimental observations in a very simple way and implies that, with appropriate control of experimental and process variables to achieve specific density and morphological conditions within the barrier films of

The permeances for N2 and CO2 for the three etched materials at these densities are shown in Figure 13a as functions of the membrane thickness using eq 13. Estimates of the permeance of Helium within the higher density structures at 1.97 g/cm3 were also obtained using the approximate value l ∼ 0.33 nm as cited above with reference to Figure 7. The Henry’s Law K value for helium at 873 K for these structures was determined to be 0.189 and, combined with the diffusivity computed from the Fickian limit shown in Figure 4, 2.461  10
8 m2/s, these results provided the predictions shown by the solid line in Figure 13a. The (ideal) permselectivities for CO2 relative to N2 defined by α¼

PCO2 PN2

ð15Þ

for each of the solids are reported in Figure 13b. The two key points which the results in Figure 13 demonstrate are as follows: (i) the permeances decrease rapidly with increasing silica film thickness in agreement with experimental observations,1,3 and (ii) although the selectivity for films greater than 30 nm is 515

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516

Langmuir

ARTICLE

composite membranes, highly selective membranes for difficult mixtures such as CO2/N2 may be fabricated.

(31) http://www.uni-graz.at/imawww/kuntsevich/solvopt/ (32) Havlin, S.; Ben-Avraham, D. Adv. Phys. 2002, 51, 187. (33) Kajihara, K.; Hirano, M.; Takimoto, Y.; Skuja, L.; Hosono, H. Appl. Phys. Lett. 2007, 91, 071904. (34) Moloney, J. Generation and Characterisation of Nano-porous Inorganic Membranes for High Temperature Gas Separation using Molecular Modelling, PhD Thesis, University College Dublin, Ireland, 2011. (35) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, United Kingdom, 1987. (36) Akter, T.; McDermott, T. C.; MacElroy, J. M. D.; Mooney, D. A.; Dowling, D. P. Langmuir, 2011, DOI:10.1021/la2031329. (37) Altemose, V. O. J. Appl. Phys. 1961, 32, 1309.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT T.C.M. and M.T.P.M. gratefully acknowledge the financial support provided by the EU Sixth Framework Programme (Grant No.: NMP3-CT-2005-456014032) and T.A. acknowledges the support of Science Foundation Ireland grant 05RFP/ ENG0044. ’ REFERENCES (1) Lee, D; Zhang, L.; Oyama, S. T.; Niu, S.; Saraf, R. F. J. Membr. Sci. 2004, 231, 117. (2) Cuffe, L.; MacElroy, J. M. D.; Tacke, M.; Kozachok, M.; Mooney, D. A. J. Membr. Sci. 2006, 272, 6. (3) Araki, S.; Mohri, N.; Yoshimitsu, Y.; Miyake, Y. J. Membr. Sci. 2007, 290, 138. (4) McCann, M. T. P. Plasma Deposition of Nanoporous Inorganic membranes for Gas Separation, PhD Thesis, University College Dublin, Ireland, 2010. (5) MacElroy, J. M. D. Mol. Phys. 2002, 100, 2369. (6) Tsapatsis, M.; Gavalas, G. R. J. Membr. Sci. 1994, 87, 281. (7) Kim, S.; Gavalas, G. R. Ind. Eng. Chem. Res. 1995, 34, 168. (8) Tsapatsis, M.; Gavalas, G. R. AIChE J. 1997, 43, 1849. (9) Fernandes, N. E.; Gavalas, G. R. Chem. Eng. Sci. 1998, 53, 1049. (10) Mooney, D. A.; MacElroy, J. M. D.; Kozachok, M.; Cuffe, L.; Tacke, M. 16th International Congress of Chemical and Process Engineering; CHISA: Prague, the Czech Republic, 2004. (11) MacElroy, J. M. D.; Friedman, S. P.; Seaton, N. A. Chem. Eng. Sci. 1999, 54, 1015. (12) Park, I.-A.; MacElroy, J. M. D. Molec. Sim. 1989, 2, 105. (13) MacElroy, J. M. D.; Raghavan, K. J. Chem. Phys. 1990, 93, 2068. (14) MacElroy, J. M. D.; Raghavan, K. J. Chem. Soc., Faraday Trans. 1991, 87, 1971. (15) Hoshen, J.; Kopelman, R. Phys. Rev. B 1976, 14, 3438. (16) McDermott, T. C. Molecular Modelling of Nanopore Structures, PhD Thesis, University College Dublin, Ireland, 2011. (17) Bhattacharya, S.; Gubbins, K. E. Langmuir 2006, 22, 7726. (18) Verlet, L. Phys. Rev. 1967, 159, 98. (19) Wolf, D.; P. Keblinski, P.; S. R. Phillpot, S. R.; Eggebrecht, J. J. Chem. Phys. 1999, 110, 8255. (20) Fennell, C. J.; Gezelter, J. D. J. Chem. Phys. 2006, 124, 234104. (21) Carre, A.; Berthier, L.; Horbach, J.; Ispas, S.; Kob, W. J. Chem. Phys. 2007, 127, 114512. (22) Schumacher, C.; Gonzalez, J.; Perez-Mendoza, M.; Wright, P. A.; Seaton, N. A. Ind. Eng. Chem. Res. 2006a, 45, 5586. (23) Schumacher, C.; Gonzalez, J.; Wright, P. A.; Seaton, N. A. J. Phys. Chem. B 2006b, 110, 319. (24) Chakravarty, C. J. Phys. Chem. B 1997, 101, 1878. (25) Keating, P. N. Phys. Rev. 1966, 145, 637. (26) von Alfthan, S.; Kuronen, A.; Kaski, K. Phys. Rev. B 2003, 68, 073203. (27) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (28) Bunde, H.; Havlin, S.Fractals and Disordered Systems; SpringerVerlag: Berlin, 1991; Chapters 2-3. (29) Kappel, F.; Kuntsevich, A. V. Comput. Optim. Appl. 2000, 15, 193. (30) Shor, N. Z. Minimization Methods for Non-Differentiable Functions; Springer Series in Computational Mathematics, Vol. 3; SpringerVerlag: Berlin, 1985. 516

dx.doi.org/10.1021/la203994v |Langmuir 2012, 28, 506–516