Two Level Model of Hydrogen Diffusion in Nanoporous Solids with

Article pubs.acs.org/JPCC

Two Level Model of Hydrogen Diffusion in Nanoporous Solids with Strongly Bound Adsorption States Carlos Rodríguez-Castellanos*,†,‡ and Javier Q. Toledo-Marín†,§ †

Facultad de Física and Instituto de Ciencia y Tecnología de Materiales, Universidad de La Habana, San Lázaro y L. Vedado, La Habana 10400, Cuba ‡ Centro de Investigación en Ciencia Aplicada y Tecnología Avanzada (CICATA−Legaria), Instituto Politécnico Nacional, México DF 11500, México § Instituto de Física, Universidad Nacional Autónoma de México, México DF 04510, México S Supporting Information *

ABSTRACT: A detailed understanding of hydrogen adsorption and diffusion in nanoporous solids is essential to design high-density storage materials. A two level lattice model with a dynamics of thermally activated molecular jumps, inspired by experimental evidence on cubic nitroprussides, is proposed for solids with strongly and weakly bound states coexisting in nanopores. The model provides a simple, unified theoretical framework to study both the thermodynamics and kinetics of hydrogen adsorption−desorption and to calculate, with minimum computational effort, quantities which can be directly compared to experimental results. In the linear and continuous limit, retarded diffusion is described by an integro-differential equation and a frequencydependent diffusion coefficient. The transport diffusion coefficient shows a very strong loading dependence, with a high coverage behavior determined by intermolecular interactions. Desorption kinetics at constant external pressure and temperature is obtained from the analytical solution of the diffusion equation and compared to the results obtained from the numerical solution of nonlinear balance equations and from Monte Carlo simulations. The relevance of nonlinear, grain-size, and fluctuations effects is discussed. Fluctuations inside pores lead to fast, subdiffusive growth of desorption times with grain size in the nanometric range.

I. INTRODUCTION

experimental evidence on hydrogen adsorption in cubic nitroprussides.9−12 Cubic nitroprussides of the family M[Fe(CN)5NO], where M = Ni, Co, show a three-dimensional network of large interconnected cavities, each with six partially naked metal centers on its surface.9,10 Their strong electrostatic interaction with guest molecules can lead to high density adsorption of other species.11 An extensive study of H2 adsorption at 77 K and pressures up to 10 bar in the series Ni1−xCox [Fe (CN)5NO] has been presented.12 According to the experimental evidence, it is expected that six hydrogen molecules can be strongly bound to the cations, and still up to six molecules could be adsorbed in the same cavity, in states of higher energy. Cubic nitroprussides have gravimetric densities that are too large and adsorption energies that are too low for practical applications in hydrogen storage systems, but controlled variation of the proportion of Ni and Co systematically modifies the strength of the electrostatic interaction and provides an excellent experimental model to evaluate its role in hydrogen adsorption. They could also be a model system to study desorption kinetics, because diffusion of hydrogen

Adsorption in light microporous solids, where hydrogen molecules are confined in very small pores or bound to unsaturated metal centers, is one of the most promising alternatives for high-density hydrogen storage. High surface areas and large adsorption enthalpies are required.1−4 Under appropriate conditions, the combined effect of dispersion, orbital, and electrostatic interactions could be both long-range and strong enough to overcome intermolecular repulsion, binding two or more molecules per adsorption center at room temperature.5,6 At low pressures some molecules could be adsorbed in strongly bound, highly localized states. At higher pressures, a large pore volume could still be available to adsorb weakly bound and less confined molecules.7,8 However, strong attraction by charge centers leads to large activation energies for intraparticle diffusion, which could become very slow or irreversible, with undesirable high desorption temperatures. In this probable scenario, to assess the potential for hydrogen storage of a given material, diffusion and adsorption should be studied together and interpreted within the same theoretical framework. The present paper focuses on hydrogen diffusion in a hypothetical nanoporous material with strongly and weakly bound adsorption states coexisting in nanopores. The discussion is based on a two level lattice model inspired by © 2015 American Chemical Society

Received: December 30, 2014 Revised: February 4, 2015 Published: March 5, 2015 7234

DOI: 10.1021/jp513000g J. Phys. Chem. C 2015, 119, 7234−7242

Article

The Journal of Physical Chemistry C

frequency dependent diffusion coefficient, whose low frequency limit gives the transport diffusion coefficient, was introduced. The relative amount of adsorbed hydrogen during desorption of a grain at constant pressure and temperature was calculated. The present work extends those studies by proposing a lattice model of adsorption and diffusion that allows us to discuss the role of intermolecular interactions, nonlinearity, grain size, and fluctuations. The integro-differential diffusion equation is derived in the linear and continuous limit. The frequency-dependent diffusion coefficient is related to the information which can be obtained from experiments of desorption kinetics. Hopping rates are calculated with transition state theory and used to discuss the loading dependence of diffusion coefficients. Desorption kinetics is studied by comparing the analytical solutions to the diffusion equation, the numerical solution of nonlinear balance equations, and the results of Monte Carlo simulations. The remaining of the paper is organized as follows. Theory is presented in section II. Numerical results for the adsorption isotherms, loading dependence of diffusion coefficients, and desorption kinetics obtained by three different methods are presented and discussed in section III. Conclusions are summarized in section IV.

molecules in cubic nitroprussides at 77 K could have a similar behavior to that shown at room temperature by materials with proportionally higher adsorption energies. Hydrogen diffusion in porous materials has been extensively investigated, both theoretically and experimentally.13,14 As the adsorbent is usually a powder, transport across macropores and interphases is very fast and adsorption−desorption rates are determined by intraparticle diffusion.15 This process can be considered isothermal, because heat propagates very fast in the framework while molecules diffuse slowly through pores, channels, or interstitials, and is usually described by a classical diffusion equation. Loading dependent diffusion coefficients can be determined from transport measurements, pulsed field gradient nuclear magnetic resonance, or quasi-elastic neutron scattering (QENS).15−18 They are also calculated with the aid of molecular dynamics simulations (MDS), in most cases considering a rigid framework and using the interaction potentials which correctly describe equilibrium adsorption.19−28 In the case of zeolites and activated carbons, with weak dispersion forces and pore diameters much larger than hydrogen molecular size, the adsorption−desorption cycle is fast and reversible, adsorption coefficients at 300 K range from 10−7 to 10−8 m2 s−1, and the activation energies are about 1−4 kJ/mol.20 As for hydrogen diffusion in metal organic frameworks (MOFs), carbon nanotubes (CNTs), and other nanoporous materials, progress in simulations is impressive,21−27 but experimental results are not so abounding.15−18 In MOFs MIL-53(Cr) and MIL-47(V), QENS experiments at 77 K, supported by MDS, show self-diffusivities over 10−7 m2 s−1, which decrease as loading increases.16 This behavior was explained by considering jumps in one and three dimensions, respectively, with activation energies of 0.6 and 1.6 kJ/mol. Many simulations show that self-diffusivities decrease, corrected diffusivities are nearly constant, and transport diffusivities increase with increased loading.25 Deviations from this behavior have also been obtained. In MOFs with wide pores, MDS lead to self-diffusion coefficients which first increase with loading and then decrease at higher pressures due to the presence of different adsorption sites.26 Hysteresis has also been found in MOFs with narrow one−dimensional channels and small pores17,18 and is generally attributed to framework flexibility allowing the opening of small windows at high pressures. However, MC simulations show that hysteresis is also possible in a rigid framework as a result of very low probability jumps from sites located at very small pores.28 Despite the impressive progress in simulations, analytical lattice model approaches to adsorption and diffusion in nanoporous solids29−32 are very useful to describe different behaviors with the same physical model. A few free parameters, with a well-defined physical meaning, can be calculated with minimum computational effort and compared to the values obtained by fitting experimental data. These methods complement each other. A coarse-grained lattice description is particularly useful to describe diffusion in the presence of deep traps, where MDS face their limits and infrequent events techniques become necessary.14 Hydrogen desorption kinetics in a porous solid with strongly and weakly bound adsorption states coexisting in the cavities was studied in a previous paper.33 A linear integro-differential equation, which takes into account the retardation of the diffusion current with respect to the concentration gradient, a consequence of molecular trapping at the strongly bound adsorption states, was obtained and solved. A complex,

II. THEORETICAL METHODS Consider a crystalline porous material containing a periodic network of identical interconnected cavities. Inside a cavity there are i = 1, ..., N identical adsorption centers, each binding niα = 0, 1 hydrogen molecules in two different states α = 1, 2 with energies εα. The energy of a configuration {niα} will be defined as N

H=

∑ (ε1ni1 + ε2ni2 + Vni1ni2) + i=1

U 2

N

∑

ni2nj2

i≠j=1

(1)

In eq 1 V and U represent, respectively, the interaction energies between two molecules adsorbed at different states in the same site and at state α = 2 in different sites. State α = 1 will be considered strongly bound to the adsorption center, and the interaction between molecules adsorbed at this state and those adsorbed at other sites is neglected. The total number of molecules in state α can take values nα = 0, ..., Nα. Although inspired in experimental evidence on cubic nitroprussides,12 this model Hamiltonian does not pretend to give a detailed description of these or other materials. Similar, alternative models can be obtained with slight modifications to account for specific information about adsorption sites and energies, obtained from computer simulations. For example, when states of type 2 cannot be ascribed only to one adsorption center, it would be easy to include in the Hamiltonian the interaction of a molecule adsorbed at state type 2 with several molecules adsorbed at states of type 1. After partial summation over occupation numbers, the grand potential per cavity for the adsorbate in equilibrium with a gaseous phase at temperature T and fugacity f can be easily expressed as a sum over the total number of molecules adsorbed in a cavity at the weakly bound states α = 2 N2

Ω = −ln

⎛ N ⎞ n −ε n − Un (n − 1)/2 ⎜ ⎟f 2 e 2 2 2 2 (1 + f e−ε1− V )n2 ⎝ n2 ⎠ =0

∑ n2

(1 + f e−ε1)N − n2 7235

(2) DOI: 10.1021/jp513000g J. Phys. Chem. C 2015, 119, 7234−7242

Article

The Journal of Physical Chemistry C All energies have been expressed in units of kBT. The equilibrium averages of adsorbed molecules are given by nα 0(f ) =

2

∂Ω ∂εα

n 0 (f ) =

∑ nα0(f ) = −f α=1

∂Ω ∂f

(3)

Explicit expressions for nα0, n0 derived from eqs 2 and 3 and used for their numerical evaluation are given in the Supporting Information. Generalization to an arbitrary number of adsorption states per site is straightforward. In a recent paper,34 a mean field procedure to take into account the influence of neighboring cavities on equilibrium adsorption has been proposed, but these corrections are not considered here. To describe diffusion we start from nonequilibrium balance equations. Let niα(r,⃗ t) be the average number of molecules adsorbed at time t in site i and state α, inside a cavity centered at r .⃗ Under well-known assumptions (see Supporting Information), the rate of change of niα(r,⃗ t) can be expressed as a balance of uncorrelated jumps between states in the same and neighboring adsorption sites ∂niα( r ⃗ , t ) = ∂t

Figure 1. Schematic two-dimensional representation of strongly bound α = 1(●), weakly bound α = 2 (○) and transient 2* (+) states for hydrogen molecules in three neighboring nanopores, and the molecular jumps (↔) contributing to diffusion.

∂n2( r ⃗ , s) ∂n ( r ⃗ , s) z2 =− 1 + ∂s ∂s N Γ21

∑ {Γjβ ,iα( r ⃗ − a ⃗ ; a ⃗)njβ( r ⃗ − a ⃗ , t )

(4)

J ⃗ ( r ⃗ , s) =

For an open system in contact with a reservoir of heat and particles, Γiα,jβ(r;⃗ a⃗) is the thermally activated over-barrier hopping rate from state type α at site i in a cavity centered at r ⃗ to the state type β at site j in a cavity centered at r ⃗ + a.⃗ In equilibrium, Γiα,jβ(r;⃗ a)⃗ is independent of r ⃗ and detailed balance conditions lead to

z2 2N Γ21v

⎧

∑ a ⃗⎨Γ22( r ⃗ − a ⃗ ; a ⃗)n2( r ⃗ − a ⃗ , s) a⃗

⎩

⎡ n 2( r ⃗ , t ) ⎤ ⎢1 − ⎥ + Γ22( r ⃗ ; a ⃗)n2( r ⃗ , s) ⎣ N ⎦ ⎡ n 2 ( r ⃗ + a ⃗ , s ) ⎤⎫ ⎢1 − ⎥⎬ ⎣ ⎦⎭ N

(8)

The equilibrium hopping rates Γ22(a) and Γ21 can be obtained from conventional transition state theory (TST).35

(5)

The description is simplified by averaging fluctuations inside a cavity n (r ⃗ , t) niα( r ⃗ , t ) = α N

⎩

A dimensionless time variable s = tΓ21 has been introduced in eqs 7a and 7b. Designating by v the volume occupied by a cavity, the current density is given by

[1 − niα( r ⃗ , t )] − Γiα , jβ( r ⃗ ; a ⃗)niα( r ⃗ , t )

Γjβ ,iα(a ⃗)nβ 0[N − nα 0] − Γiα , jβ(a ⃗)nα 0[N − nβ 0] = 0

a⃗≠0

⎡ n (r ⃗ , t) ⎤ ( r ⃗ − a ⃗ , t )⎢ 1 − 2 ⎥ − Γ22( r ⃗ ; a ⃗)n2 ⎣ N ⎦ ⎡ n ( r ⃗ + a ⃗ , t ) ⎤⎫ ( r ⃗ , t )⎢ 1 − 2 ⎥⎬ ⎣ ⎦⎭ N (7b)

j,β ,a⃗

[1 − njβ ( r ⃗ + a ⃗ , t )]}

⎧

∑ ⎨Γ22( r ⃗ − a ⃗ ; a ⃗)n2

Γ22(a) =

u P(2*) a P(2)

Γ21 =

u b

(9)

In eq 9 the hopping distances are a, b and the thermal velocity is u. It is assumed that no potential barrier needs to be crossed in transitions between states at the same site. In the expression for the transition rate Γ22(a) between weakly bound states in neighboring cavities, P(2), P(2*) are, respectively, the equilibrium occupation probabilities of the initial and transient states 2, 2*. The last consists of a molecule removed from state 2 at site j inside a cavity centered at r ⃗ and placed at the top of the potential barrier in the channel connecting it with the nearest site in a neighbor cavity centered at r ⃗ + a ⃗ . This approximation overestimates transition rates, because all molecules reaching the transient state with nonzero velocity are assumed to jump successfully to the final state. Different ways to improve TST by including a dynamical correction factor have been proposed,36−38 but they will not be considered here. The equilibrium occupation probabilities of the initial and transition states are

(6)

Since molecules at states α = 1 are strongly bound, only jumps to weakly bound states α = 2 in the same site with transition rate Γ12 will be considered. In contrast, molecules adsorbed at states α = 2 are allowed to jump to void states α = 1 in the same site with rate Γ21 and also to other void states α = 2 in the nearest sites of the same and neighboring cavities with rate Γ22(r;⃗ a⃗), as illustrated in Figure 1. Let z be the number of different cavities to which a molecule can jump from state α = 2 at a given site. Then, for a cubic array of identical cavities, where z = 3 and Γ22(r;⃗ a⃗) = Γ22(r;⃗ a) ⎡ n20(N − n10) ∂n1( r ⃗ , s) n ( r ⃗ , s) ⎤ = n 2 ( r ⃗ , s ) ⎢1 − 1 ⎥− ⎣ ∂s N ⎦ n10(N − n20) ⎡ n ( r ⃗ , s) ⎤ n1( r ⃗ , s)⎢1 − 2 ⎥ ⎣ N ⎦ (7a)

P(2) = 7236

n20 N

P(2*) = ⟨e−ΔH ⟩

(10)


Article

The Journal of Physical Chemistry C where ⟨...⟩ is the grand canonical average over states with nj2 = 1 and ΔH = H[{niα − δi , jδα ,2}] − H[{niα}]

D̃ (ω) = D2τ0

(11)

(12)

Then, eqs 7a, 7b, and 8 become ∂ρ1 ∂s ∂ρ2 ∂s

=

N − n10 n ρ2 − 20 ρ1 N − n20 n10

=−

∂ρ1 ∂s

(13a)

δn(s) = + D2∇2 ρ2

(13b)

J ⃗ ( r ⃗ , s) = −D2∇ρ2

(14)

n ⎞ ∂D ⎛ D2 = Ds + n20⎜1 − 20 ⎟ s ⎝ N ⎠ ∂n20

(15)

z 2 Γ22 2 a N Γ21

n(s) − n(∞) n(0) − n(∞)

I (̃ ω) = −D̃ (ω)

s

ds′e−s ′∇ρ( r ⃗ , s − τ0s′)

N − n10 n τ0−1 = + 20 N − n20 n10

(23)

∬ dS ⃗·∇ρ ̃( r ⃗ , ω)

(24)

To evaluate the surface integral, one has to solve the diffusion equation in Fourier space. Boundary conditions correspond to equilibrium between hydrogen adsorbed at the grain surface ∑ and the external reservoir

(17)

D̃ (ω)∇2 ρ ̃( r ⃗ , ω) = iωρ ̃( r ⃗ , ω)

(18)

ρ ̃( r ⃗ , ω)|Σ = −

where ρ = ρ1 + ρ2 is the density deviation from equilibrium and τ0 is the relaxation time inside a cavity. For times much larger than τ0 or for a time independent concentration gradient, the current response is governed by a transport diffusion coefficient Dt always smaller than D2: ⎡ n N − n10 ⎤ Dt = D2⎢1 + 10 ⎥ n20 N − n20 ⎦ ⎣

(22)

Now, I(̃ ω) is calculated by integrating eq 18 over the grain boundary:

(16)

∫0

(21)

I (̃ ω) = 1 − 2πiωδn(̃ ω) ̃ I (0)

By combining eqs 13a and 14, a retarded current response to a density gradient is obtained. J ⃗ ( r ⃗ , s) = −Dt

iωτ0 + 1

Let I(̃ ω) and δñ(ω) be, respectively, the Fourier transforms of I(s) and δn(s). It is easy to show that

where Ds is the self-diffusion coefficient Ds =

n20 n10

Obviously, retardation becomes important at low pressures, when n20 ≪ n10. The transport diffusion coefficient Dt is the low frequency limit of D̃ (ω). According to eq 19, it depends on loading through D2 and the concentration dependence of retardation inside a cavity. In the general case, a complex diffusion coefficient can be obtained from experiments of adsorption or desorption kinetics. Consider, for example, desorption at constant external pressure of a grain. Let n(s) be the number of molecules which remain adsorbed after time s and I(s) = −ṅ(s) the flux rate. The relative amount δn(s) of hydrogen adsorbed after time s is given by

Detailed expressions for P(2*) to be used in numerical calculations are given in the Supporting Information. A linear diffusion equation can be obtained for systems not far from equilibrium by passing to a linear and continuous description in the vicinity of the equilibrium distribution. Let ρα be the deviation from equilibrium density at states type α nα( r ⃗ , s) = nα 0 + vρα ( r ⃗ , s)

iω +

(25)

Δρ 2πi(ω − iϵ)

(26)

This problem can be easily solved for simple grain geometries. In the case of spherical grains of radius R it can be shown that (see Supporting Information) I (̃ ω) 3 = 2 [1 − γ cot γ ] ̃ I (0) γ

−1

(19)

(27)

⎡ ωR2 ⎤1/2 γ=⎢ ⎥ ⎣ iD(ω) ⎦

By combining eqs 2−3 and 9−11 the loading dependence of the diffusion coefficients Ds and Dt defined by eqs 16 and 19 can be easily obtained (see Supporting Information). A “corrected” diffusion coefficient can also be calculated dividing Dt by the corresponding thermodynamic factor14 derived from the adsorption isotherms. Diffusion is retarded because the model assumes two different states for storage, but only one for diffusion. Molecules trapped at the strongly bound states α = 1 must jump first to the weakly bound states α = 2 and only after that can diffuse to other cavities. Retardation gives rise to a complex, frequencydependent diffusion coefficient D̃ (ω) and Fick’s law becomes a linear relation between the Fourier transforms of current density and density gradient → ∼⎯ J ( r ⃗ , ω) = −D̃ (ω)∇ρ ̃( r ⃗ , ω) (20)

(28)

In the low frequency or small grain limit ωR ≪ |D(ω)| one has 2

R2 D̃ (ω) = [δn(̃ ω)]−1 30π

(29)

The transport diffusion coefficient can be derived from experimental data by using eq 29 and the results compared to the theoretical value given by eq 21. To describe adsorption or desorption kinetics, the time dependence of directly measurable quantities like δn(s) and I(s) can be derived by solving eqs 13a and 13b. The set can be reduced to a single integro-differential diffusion equation,33 with the integral term describing retardation (see Supporting Information) 7237


Article

The Journal of Physical Chemistry C n N − n10 ∂ρ( r ⃗ , s) − 20 n10 N − n20 ∂s

∫0

s

to void states j, 2 in a different site of the same cavity. Transition probabilities are taken to be one when Hi > Hf, and w = eHi−Hf when Hi < Hf. Next, transitions from occupied states i, 2 to nearest void states j, 2 in a neighboring cavity have probability w = eHi−H*i, where H*i is the energy of the transient state. Finally, when the initial state is in a nanopore located at the boundary and the final state is inside the grain, a molecule is created with probability w = eHi−H*i+μ. In the opposite case, a

ds′ρ( r ⃗ , s − s′)e−n20 / n10s ′

= D2∇2 ρ( r ⃗ , s)

(30)

Consider a grain charged at f 0, which for t > 0 is connected to a reservoir of fugacity f at the same temperature T. The solution of eq 30 with the corresponding initial and boundary conditions is ρ( r ⃗ , s) = Δρ[∑ ⟨uk⟩uk( r ⃗)Gk (s) − 1]

molecule is annihilated with probability w = eHi−H*i−μ, where μ is the chemical potential of the gaseous phase. This approach preserves the initially proposed model dynamics, described by eq 4 without neglecting fluctuations inside a cavity, as implied by eq 6.

(31)

k

where the function Gk(s) is given by G k (s ) =

Dt λk ⎧⎡ D2λk − bk ⎤ −aks ⎡ ak − D2λk ⎤ ⎨⎢ +⎢ ⎥e ⎥ (ak − bk ) ⎩⎣ ak bk ⎦ ⎣ ⎦ ⎪

⎪

⎫ e−bks⎬ ⎭

III. RESULTS AND DISCUSSION In this section the results of numerical calculations for a model system with N = 6 sites, and maximum occupation numbers N1 = 6 and N2 = 5 per cavity are presented and discussed. For illustrative purposes, binding energies will be taken ε1 = −6 and ε2 = −2 and parameters U and V will take different values to account for attractive or repulsive intermolecular interactions. Adsorption isotherms are shown in Figures 2 and 3. The equilibrium averages of the number of adsorbed molecules as a

⎪

⎪

(32)

⎛ ak ⎞ 1 ⎧ −1 ⎜ ⎟ = ⎨(τ0 + D2λk ) 2⎩ ⎝ bk ⎠ ⎪

⎪

±

(τ0−1 + D2λk )2 − 4D2λk

n20 ⎫ ⎬ n10 ⎭ ⎪

⎪

(33)

In eqs 31−33, λk and uk(r ⃗ ) are, respectively, the eigenvalues and normalized eigenfunctions of a boundary value problem (see Supporting Information). The eigenvalues λk only depend on the shape and size of a grain. ∇2 uk + λk uk = 0

⟨uk⟩ =

uk( r ⃗ , s)|Σ = 0

∭ uk*( r ⃗) d3r ⃗

(34) (35) Figure 2. Adsorption isotherms. Average number of molecules adsorbed as a function of fugacity in a cavity with 6 sites and maximum occupation numbers N1 = 6, N2 = 5, for ε1 = −6, ε2 = −2 and different values of the interaction energies U, V. All energies are in units of kBT.

From eq 31, the relative amount δn(s) of hydrogen which remains adsorbed in a grain of volume V after time s and the flux rate I(s) is given by δn(s) =

1 V

∑ |⟨uk⟩|2 Gk(s) k

I(s) = −Δρ ∑ |⟨uk⟩|2 Ġ k (s) k

(36)

Equations 30−36 are equivalent to those obtained in ref 33, except for the change of self-diffusion coefficient by the loading dependent coefficient D2 defined by eq 15. Nonlinear and small grain size effects, not included in eq 30, can be considered by using the numerical solution of eqs 7a and 7b for a grain containing a finite array of interconnected cavities. A comparison between desorption kinetics obtained from eqs 30 and 7a−7b can be used to consider the role of nonlinear effects and the continuous limit in the case of very small grains. Desorption kinetics can also be described by combining Monte Carlo simulations for a finite array of interconnected cavities with detailed balance and TST relations between transition rates. Starting with saturated cavities, consider random molecular jumps between occupied and void states in the same or neighboring cavities. Let Hi, Hf be, respectively, the energies of two successive configurations. A molecule in an occupied state i, α can jump to a void state i, β in the same cavity and site. Molecules at occupied states i, 2 can also jump

Figure 3. Adsorption in weakly bound states. Average number of molecules adsorbed at the weakly bound states as a function of fugacity in a cavity with 6 sites, maximum occupation numbers N1 = 6, N2 = 5, for ε1 = −6, ε2 = −2 and different values of the interaction energies U, V.

function of fugacity have been calculated from eqs 2−4. All isotherms show type I behavior. Strongly bound states α = 1 are fully saturated at low pressures and intermolecular interactions have little effect on their occupation. Weakly bound states α = 1 are occupied at higher pressures, especially if interactions are repulsive (U, V > 0). Other values of the parameters could be obtained by fitting experimental data. To compare with excess 7238


Article

The Journal of Physical Chemistry C adsorption isotherms,12 fugacity must be translated into pressure values and a correction must be subtracted to account for the amount of gas in the cavity volume at given temperature and pressure. The loading dependence of diffusion coefficients Ds and Dt was calculated by combining eqs 2 and 3 with 15, 16, and 19. Figure 4 shows the dependence of the self-diffusion coefficient

33, where the loading dependence of the self-diffusion coefficient and the spatial variation of transition rates under nonequilibrium conditions were not considered. Desorption kinetics has been studied by three different methods. Consider desorption at 298 K from a cubic grain with 103 cavities. The relative amounts δn(s) of hydrogen which remains adsorbed after time s, calculated by solving the linear integro-differential diffusion eq 22 and nonlinear balance eqs 7a, 7b, are shown in Figure 6. To obtain almost complete

Figure 4. Loading dependence of the self-diffusion coefficient. Calculated at 298 K for cavities with 6 sites, N1 = 6, N2 = 5, ε1 = −6, ε2 = −2 and different values of the interaction energies U, V.

Figure 6. Desorption kinetics from diffusion and balance equations. The relative amount of hydrogen adsorbed after time s during desorption to a final equilibrium state at T = 298 K and μ = −12 in a cubic grain with 1000 cavities, each with 6 sites, and maximum occupation numbers N1 = 6, N2 = 5, has been calculated for ε1 = −6, ε2 = −2 from the solutions of the linear diffusion equation (DE) and nonlinear balance equations (BE). The chemical potentials corresponding to the initial conditions and the calculated desorption times are shown in the inset.

Ds with the average number n0 of molecules adsorbed in a cavity. The behavior is governed by the character of the interactions between adsorbed molecules. In the absence of interactions (U = V= 0), Ds is independent of n0. When all interactions are repulsive (U, V > 0), the presence of other molecules makes it easier for a molecule to jump to a neighbor cavity, leading to an increase of Ds. The effect of attractions (U, V < 0) is the opposite. When repulsion and attraction coexist, the dependence is non-monotonic: the low coverage (n0 < 6) behavior is controlled by the character of the interaction between molecules adsorbed at the same site (V). At higher loadings, when states α = 2 are populated, the attractive or repulsive character of the interaction between weakly bound molecules (U) becomes dominant. In any case, the selfdiffusion coefficient Ds has the same order of magnitude for all values of coverage. The loading dependence of the transport diffusion coefficient Dt is stronger, as shown in Figure 5. Notice that there is a

evacuation, the chemical potential of the gas reservoir has been taken μout = −12. No dependence of the interaction parameters U, V is found, because all transition probabilities have been evaluated for the final equilibrium state, where loading is very small and intermolecular interactions become irrelevant. The upper curve in Figure 6 corresponds to the solution of the linear integro-differential diffusion equation. In the linear approximation the relative quantity δn(s) is independent of the initial loading. The other curves were obtained by solving nonlinear balance equations for different initial conditions. The lower one (μin = −8, 5) corresponds to a grain loaded to almost complete saturation. As the initial concentration is reduced, while keeping constant the final value, the desorption curves obtained from balance equations approach those predicted by the linear diffusion equation. A dimensionless desorption time τ = 2πδñ(0) can be obtained from the areas under the curves. Its dependence with grain size, presented in Figure 7, is characteristic of normal

Figure 5. Loading dependence of the transport diffusion coefficient. Calculated at 298 K for cavities with 6 sites, N1 = 6, N2 = 5, ε1 = −6, ε2 = −2 and different values of the interaction energies U, V.

general increase by 1 or 2 orders of magnitude, because the number of molecules able to diffuse at low concentrations is very small and there is almost no current response to a concentration gradient. As n20 increases, so does the transport diffusion coefficient, but at high coverage the behavior is again dominated by the attractive or repulsive character of the interactions between molecules at weakly bound adsorption states. This behavior is very different from that obtained in ref

Figure 7. Grain size dependence of desorption time. For cubic grains with L3 cavities τ has been calculated from the solutions of balance and diffusion equations. Diffusion is normal and desorption times obtained from diffusion equation are systematically larger. 7239


Article

The Journal of Physical Chemistry C diffusion (proportional to L2). Desorption times obtained by solving the nonlinear balance equations are systematically smaller than the value predicted by the linear integrodifferential diffusion equation. The continuous limit has no special significance, since the ratio between those desorption times is essentially independent of grain size. Thus, the difference can be attributed to nonlinear effects included in the balance equations. The linear approximation corresponds to diffusion with the low coverage value of the transport diffusion coefficient, much smaller than the average value during desorption, as shown in Figure 4. When solving balance equations, the loading dependence of the number of transitions per unit time, associated with the impossibility of a molecule to jump to an occupied state, is explicitly taken into account. As the initial concentration is reduced, desorption times obtained from balance equations approach the value predicted by the linear integro-differential diffusion equation. The results of Monte Carlo simulations of desorption, following the rules described in the previous section, are shown in Figures 8 and 9. Now, intermolecular interactions enter the

type 2, close to the grain boundary, are released very rapidly. Then, it follows desorption of the remaining molecules at states of type 2. Finally, molecules at strongly bound states are slowly desorbed. A desorption time τα = 2πδñα(0) now expressed in MC steps can be obtained from the areas under the curves. The results for L = 10 are presented in Table 1 and show the expected Table 1. Desorption Times from MCSa U

V

τ1

τ2

τ

0 0.25 −0.25 0.25 −0.25

0 0.5 0.5 −0.25 −0.25

1649 1615 1739 1615 1868

285 197 282 262 413

1029 971 1077 1000 1207

In a cubic grain with 103 cavities, desorption times τ1, τ2, τ (103 MC steps) have been obtained for different values of intermolecular interaction energies. a

dependence with U and V. We have also calculated desorption times τ for values of L between 5 and 25 in the case U = V = 0. Fitting the data to a power law L2 = Aτx leads to x ≈ 0.38 as shown in Figure 10.

Figure 8. Desorption kinetics from MCS. Relative amount of hydrogen adsorbed after t MC steps during desorption from saturation to a final equilibrium state at T = 298 K and μ = −12 in a cubic grain with 1000 cavities, each with 6 sites, and maximum occupation numbers N1 = 6, N2 = 5, calculated from MCS for ε1 = −6, ε2 = −2 and different values of the interaction energies U and V.

Figure 10. Grain size dependence of desorption time τ (in MC steps). Calculated from MCS for cubic grains with L3 cavities.

The value x ≈ 0.38 corresponds to “anomalous subdiffusive” behavior.39 It contrasts with the “normal” diffusion (x = 1) described by balance and diffusion equations. The terms “normal” or “anomalous” diffusion should be used with caution, since they correspond, strictly speaking, to the limit behavior of infinite systems, and here we are dealing (as in practical applications) with finite, quite small submicronic grains. In all three methods employed to describe desorption, the dynamics is apparently the same: diffusion is retarded due to the capture of molecules at the strongly bound states. However, when MCS are implemented, a geometrical constraint to Brownian motion, neglected in balance and diffusion equations, is automatically imposed: two successive molecular jumps in the same direction are impossible. The point is that a molecule, coming from a cavity at r ⃗ − a⃗ to a cavity at r,⃗ cannot do its next jump to the cavity at r ⃗ + a⃗ because it has no near neighbors there. It must previously jump to another empty site in the same cavity (where it has some probability of been trapped at the strongly bound state) and only then to the cavity at r ⃗ + a⃗, as sketched in Figure 1. This intermediate step increases waiting times and was not considered in balance and diffusion equations as a consequence of spatial averaging over fluctuations in the scale of a cavity (eq 6). On the other hand, fluctuations in the occupation numbers of adsorption

Figure 9. Desorption kinetics for strongly and weakly bound states. Relative amount of hydrogen adsorbed at different states after t MC steps during desorption from saturation to a final equilibrium state at T = 298 K and μ = −12 in a cubic grain with 1000 cavities, each with 6 sites, and maximum occupation numbers N1 = 6, N2 = 5, calculated from MCS for ε1 = −6, ε2 = −2 and different values of the interaction energies U, V.

energy balance to determine the acceptance of a new configuration at successive MC steps. As expected, diffusion is slower in the presence of attractive interactions (U < 0, V < 0) with respect to the repulsive case (U > 0, V > 0). The curves corresponding to mixed cases (U > 0, V < 0; U < 0: V > 0 and U = V=0) lie between those shown in the figures. From Figure 9 it can be appreciated that desorption takes place in three stages. Initially, molecules adsorbed in states of 7240



■

states lead to random values of the transition rates between them. These two differences, a broader distribution of waiting times and disorder related to fluctuations, are responsible for the “anomalous subdiffusive” behavior observed in MC simulations. This is expressed in desorption times increasing with grain size approximately as L5.3 in the considered range. A general random walk theory for diffusion in the presence of nanoscale confinement has been developed and applied in ref 40. Three time scales are identified, corresponding to normal diffusion inside the confining cage, subdiffusive behavior for intermediate times, and again normal diffusion (with smaller diffusivity) in the long time limit. These results for singleparticle diffusion cannot be easily extended to transport diffusion during desorption from a saturated grain. In the present case, molecular jumps are limited not only by intercavity potential barriers, but also by geometrical constraints, traps, interactions, and random site occupation inside cavities. We have not found the third time scale in our MCS of grains with up to 125 000 cavities. For practical purposes this means that, in the submicron range, the model predicts much slower desorption in larger grains. We do not know, for sure, whether or not subdiffusive behavior, with the same or a different exponent, should prevail for very large grains. A different situation arises when, instead of diffusion in a single grain, we consider a very large storage vessel containing many nanosized grains, together with larger intergrain pores where diffusion is faster. Now the picture described in ref 40 is partially recovered. Ordinary diffusion will take place at times much longer than the single-grain desorption time, with some effective diffusion coefficient limited by subdiffusive behavior inside grains.

Article

ASSOCIATED CONTENT

S Supporting Information *

Details of long, but straightforward calculations are given in the Supporting Information accompanying this paper. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 537 + 8762099. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the support of an “Alma Mater” project from the University of Havana and a “SECITIDFCLAF” project of the Latin-American Centre of Physics and the Government of Mexico City. Recommendations derived from discussions with Dr. Edilso Reguera Ruiz from CICATA, IPN, Mexico, are also highly appreciated.

■

REFERENCES

(1) Jena, P. Materials for Hydrogen Storage: Past, Present, Future. J. Phys. Chem. Lett. 2011, 2, 206−211. (2) Bhatia, S. K.; Myers, A. L. Optimum Conditions for Adsorptive Storage. Langmuir 2006, 22, 1688−1700. (3) Lochan, R. C.; Head −Gordon, M. Computational Studies of Molecular Hydrogen Affinities. The Role of Dispersion Forces, Electrostatic and Orbital Interactions. Phys. Chem. Chem. Phys. 2006, 8, 1357−1370. (4) Niu, J.; Rao, B. K.; Jena, P. Binding of Hydrogen Molecules by a Transition Metal Ion. Phys. Rev. Lett. 1992, 68, 2277−2280. (5) Tuan, K. H.; Antonelli, D. M. Exploiting the Kubas Interaction in the Design of Hydrogen Storage Materials. Adv. Mater. 2004, 21, 1787−1800. (6) Reguera, E. Materials for Hydrogen Storage in Nanocavities: Design Criteria. Int. J. Hydrogen Energy 2009, 34, 9163−9167. (7) Liu, Y.; Brown, C. M.; Neumann, D. A.; Peterson, V. K.; Keper, C. J. Inelastic Neutron Scattering of H2 Adsorbed in HKUST-1. J. Alloys Compd. 2007, 385, 446−447. (8) Panella, B.; Hones, K.; Muller, U.; Trukhan, N.; Schubert, M.; Putter, H.; Hirscher, M. Desorption Studies of Hydrogen in MetalOrganic Frameworks. Angew. Chem., Int. Ed. 2008, 47, 2138−2142. (9) Oskarsson, F.; Aradottir, E. S.; Jonsson, H. Calculations of Release Temperature of Hydrogen Storage Materials. Prepr. Pap. − Am. Chem. Soc., Div. Fuel Chem. 2006, 51, 1. (10) Balmaseda, J.; Reguera, E.; Gómez, A.; Roque, J.; Vázquez, C.; Autié, M. On the Microporous Nature of Transition Metal Nitroprussides. J. Phys. Chem. B 2003, 107, 11360−11369. (11) Culp, J. T.; Matranga, C.; Smith, M.; Bittner, E. J.; Bockrath, B. Hydrogen Storage Properties of Metal Nitroprussides M[Fe(CN)5NO], (M =Co, Ni). J. Phys. Chem. B 2006, 110, 8325−8328. (12) Reguera, L.; Roque, J.; Hernández, J.; Reguera, E. High Density Hydrogen Storage in Nanocavities: Role of the Electrostatic Interaction. Int. J. Hydrogen Energy 2010, 35, 12864−12869. (13) Roque-Malherbe, R. Adsorption and Diffusion in Nanoporous Materials; CRC Press, Taylor and Francis Group: Boca Raton, 2007. (14) Kärger, J.; Ruthven, D. M.; Theodorou, D. N. Diffusion in Nanoporous Materials; Wiley-VCH Verlag GmbH &Co. KGaA: Weinheim, 2008. (15) Zamora, B.; Al-Hajjaj, A. A.; Shah, A. A.; Bavykin, D. V.; Reguera, E. Kinetics and Thermodynamics of Adsorption on Titanate Nanotubes Decorated by a Prussian Blue Analogue. Int. J. Hydrogen Energy 2013, 38, 6406−6416. (16) Salles, F.; Jobic, H.; Maurin, G.; Koza, M.; Llewel-lyn, P. L.; Devic, T.; Serre, C.; Ferey, G. Experimental Evidence Supported by

IV. CONCLUSIONS The present work provides a useful tool to understand the contradictory relationship between adsorption and diffusion in materials designed for hydrogen storage. By combining a lattice model with transition state theory, a simple, unified theoretical framework to study both the thermodynamics and kinetics of hydrogen adsorption−desorption in solids with strongly and weakly bound states coexisting in nanopores is obtained. Minimum computational effort is needed to calculate quantities which can be directly compared to experimental results. It allows us to discuss the influence of a wide set of microscopic (energy and degeneracy of adsorption states, character and magnitude of intermolecular interaction energies), and macroscopic variables (pressure, temperature, grain shape, and size). Similarly, alternative models can be obtained with slight modifications to account for complementary information about the adsorption sites and energies of hydrogen or other species in nanoporous solids. The strong loading dependence of diffusion coefficients evidence the nonlinear character of diffusion, which is better described in the whole loading interval by MCS. The presence of two different microscopic time and length scales, one determined by intracavity phenomena and the other by intercavity transitions, leads to different regimes according to the relationship between them. Normal diffusion results from neglecting spatial fluctuations inside a cavity. Nonretarded diffusion is obtained as the intracavity relaxation time vanishes. The practical importance of a detailed description of the dynamics inside the pores, including the role of fluctuations, is confirmed. 7241


Article


(37) Dubbeldam, D.; Beerdsen, E.; Vlugt, T. J.; Smit, B. Molecular Simulation of Loading-Dependent Diffusion in Nanoporous Materials Using Extended Dynamically Corrected Transition State Theory. J. Chem. Phys. 2005, 122, 224712−224728. (38) Beerdsen, E.; Dubbeldam, D.; Smith, B. Molecular Understanding of Diffusion in Confinement. Phys. Rev. Lett. 2005, 95, 164505−164508. (39) Bouchaud, J. P.; Georges, A. Anomalous Diffusion in Disordered Media: Statistical Mechanisms, Models and Physical Applications. Phys. Rep. 1990, 195, 127−293. (40) Calvo-Muñoz, E. M.; Esai Selvan, M.; Xiong, R.; Ojha, M.; Keffer, D. J.; Nicholson, D. M.; Egami, T. Applications of a General Random Walk Theory for Confined Diffusion. Phys. Rev. E 2011, 83, article #01120.

Simulations of a Very High H2 Diffusion in Metal Organic Framework Materials. Phys. Rev. Lett. 2008, 100, 245901−245904. (17) Joon, G. K.; Su-Huai, W.; Yong-Hyun, K. Microscopic Theory of Hysteretic Hydrogen Adsorption in Nanoporous Materials. J. Am. Chem. Soc. 2010, 132, 1510−1511. (18) Xuebo, Z.; Bo, X.; Fletcher, A. J.; Thomas, K. M.; Bradshaw, D.; Rosseinsky, M. J. Hysteretic Adsorption and Desorption of Hydrogen by Nanoporous Metal-Organic Frameworks. Science 2004, 306, 1012− 1015. (19) Sholl, D. S. Understanding Macroscopic Diffusion of Adsorbed Molecules in Crystalline Nanoporous Materials via Atomistic Simulations. Acc. Chem. Res. 2006, 39, 403−411. (20) Dubbeldam, D.; Snurr, R. Q. Recent Developments in the Molecular Modeling of Diffusion in Nanoporous Materials. Mol. Simul. 2005, 33, 305−325. (21) Skoulidas, A. I.; Sholl, D. S. Self-Diffusion and Transport Diffusion of Light Gases in Metal-Organic Framework Materials Assessed Using Molecular Dynamics Simulations. J. Phys. Chem. B 2005, 109, 15760−15678. (22) Yang, Q.; Zhong, Ch. Molecular Simulation of Adsorption and Diffusion of Hydrogen in Metal-Organic Frameworks. J. Phys. Chem. B 2005, 109, 11862−11864. (23) Liu, J.; Lee, J. Y.; Pan, L.; Obermyer, R. T.; Simizu, S.; Zande, B.; Li, J.; Sankar, S. G.; Johnson, J. K. Adsorption and Diffusion of Hydrogen in a New Metal-Organic Framework Material: [Zn (bdc)(ted)0.5]. J. Phys. Chem. C 2008, 112, 2911−2917. (24) Liu, B.; Yang, Q.; Xue, C.; Zhong, C.; Smith, B. Molecular Simulation of Hydrogen Diffusion in Interpenetrated Metal-Organic Frameworks. Phys. Chem. Chem. Phys. 2008, 10, 3244−3249. (25) Zhiqiang, W.; Zhiping, L.; Wenchuan, W.; Yiqun, F.; Xu, N. Diffusion of H2, CO, N2, O2 and CH4 through Nanoporous Carbon Membranes. Chinese Journal of Chemical Engineering 2008, 16, 709− 714. (26) Garbeoglro, G.; Vallauri, R. Adsorption and Diffusion of Hydrogen and Methane in 2D Covalent Organic Frameworks. Microporous Mesoporous Mater. 2010, 116, 540−547. (27) Jinchen, L.; Rankin, R. B.; Johnson, J. K. The Importance of Charge-Quadrupole Interactions for H2 Adsorption and Diffusion in Cu-btc. Mol. Simul. 2009, 35 (1−2), 60−69. (28) Di-Zhang, Z.; Wei-Xiong, Z.; Feng-Lei, C.; Long, J.; Tong-Bu, L. A Three-Dimensional Microporous Metal−Organic Framework with Large Hydrogen Sorption Hysteresis. Chem. Commun. 2011, 47, 1204−1206. (29) Elliott, J. A. W.; Ward, C. A. Temperature Programmed Desorption: A Statistical Rate Theory Approach. J. Chem. Phys. 1997, 106, 5677−5684. (30) Adhangale, P.; Keffer, D. J. Obtaining Transport Diffusion Coefficients from Self-Diffusion Coefficients in Nanoporous Adsorption Systems. Mol. Phys. 2004, 102, 471−483. (31) Kamat, M.; Keffer, D. An Analytical Theory for Diffusion of Fluids in Crystalline Nanoporous Materials. Mol. Phys. 2003, 101, 1399−1412. (32) Kamat, M.; Weijing, D.; Keffer, D. Agreement between Analytical Theory and Molecular Dynamics Simulation for Adsorption and Diffusion in Crystalline Nanoporous Materials. J. Phys. Chem. B 2004, 108, 376−386. (33) Rodríguez, C.; Reguera, E.; Cabrera, R. Hydrogen Diffusion in Nanoporous Solids with Strongly Localized Adsorption States. Revista ́ 2011, 28, 37−41. Cubana de Fisica (34) Pazzona, F. G.; Demontis, P.; Suffritti, G. B. Coarse Graining of Adsorption in Microporous Materials: Relation between Occupancy Distributions and Local Partition Functions. J. Phys. Chem. C 2014, 118, 28711−28719. (35) Chandler, D. Statistical Mechanics of Isomerization Dynamics in Liquids and the Transition State Approximation. J. Chem. Phys. 1978, 68, 2959−2970. (36) Beersden, E.; Smit, B.; Dubbeldam, D. Molecular Simulation of Loading Dependent Slow Diffusion in Confined Systems. Phys. Rev. Lett. 2004, 93, 248301−248304. 7242


Two Level Model of Hydrogen Diffusion in Nanoporous Solids with

Recommend Documents