Multimodels computation for adsorption capacity of activated carbon

Abstract

The pore size distribution of activated carbon is conventionally characterized with nitrogen adsorption measurements at 77 K. The adsorption isotherms are commonly analyzed with a nonlocal density functional theory in combination with a mathematical model for the pore size and geometry. While nonlocal density functional theory is significantly more accurate than the Brunauer–Emmett–Teller theory for gas adsorption, its application to materials characterization is mostly based on a mean-field approximation for van der Waals attractions that is only qualitative in comparison with alternative versions of nonlocal density functional theory or molecular simulations. Toward development of a more reliable theoretical procedure, we compare mean-field approximation-nonlocal density functional theory with three recent versions of non-mean-field methods for gas adsorption at conditions corresponding to experiments for porous materials characterization. The potential applicability of different nonlocal density functional theory methods for pore size distribution predictions is evaluated in terms of the theoretical error bound scale analysis. We find that the weight density approximation is the most reliable for predicting the pore size distribution of amorphous porous materials. In addition to accurate isotherm, weight density approximation yields the theoretical error bound scale for pore size distribution prediction nearly 10⁴ times narrower than that corresponding to mean-field approximation. The new theoretical procedure has been used to analyze the pore size distribution of four activated carbon samples and to predict the adsorption capacities of these materials.

Introduction

Activated carbons can be classified as one kind of amorphous porous materials playing important roles in diverse engineering applications (Herdes et al., 2016; Nikolaidis et al., 2016; Ray and Roy, 2013). Because the performance of such material is critically dependent on the pore size and specific surface area, structure characterization is essential before their practical use. However, direct nanoscale analysis is difficult for industrial-level applications.

Pore size distribution (PSD) represents one of the simplest and most broadly used benchmark for characterization of amorphous porous materials (Balzer et al., 2015; Palmer et al., 2011; Wang et al., 2015). PSD is commonly acquired from the isotherm fitting method, not only because of its wide detection range, but also for its reliability and low cost. Specifically, PSD is normally determined from fitting N₂ adsorption isotherm at 77 K with the assumptions of a given pore geometry and a suitable mathematical model. While the procedure has been standard in materials characterization, controversies often arise in terms of its fidelity to represent the real material structure and properties beyond adsorption isotherms. On the one hand, there is an apparent diffusion problem at low temperature (Madani et al., 2016), which is consequently affected by the working capacity characterization under practical conditions. On the other hand, N₂ molecules are more easily subject to reacting with the surface groups than nonpolar adsorbates, leading to the overprediction of the surface area.

PSDs fitted with a particular adsorption isotherm have been successfully utilized to predict gas adsorptions at other conditions (Sitprasert et al., 2011; Sweatman and Quirke, 2005; Wang and Do, 1997). However, the effective PSD may not be able to represent the physical structure of real porous materials. As PSD is typically obtained by solving the adsorption integral equation (AIE) according to the adsorption isotherm, the outcome for such method is primarily affected by the kernel matrix, which is treated as an a priori hypothesis on the heterogeneous extent of the amorphous materials. Besides, an accurate prediction of the excess adsorption in the confined space remains a challenge since the atomistic structure at the surface is poorly understood. For gas adsorption in small pores, conventional thermodynamic methods and even some statistical-mechanics theories like the nonlocal density functional theory (NLDFT) with the mean-field approximation (MFA) are proven unreliable.

NLDFT has been integrated into the standard software for materials characterization. While NLDFT is built upon an exact foundation, i.e. unique determination of an equilibrium system in terms of particle density distributions, its practical use relies on approximations for the thermodynamic nonideality due to intermolecular interactions. For that purpose, the standard procedure employs a nonlocal weighted density approximation for short-range repulsion and the MFA for longer range attraction. To apply similar ideas for materials characterization, there have been a number of modifications. From a theoretical perspective on the one hand, one may modify the excess Helmholtz free energy functional using weight density approximation (WDA) and density expansion method (Liu et al., 2010; Yu, 2009; Zeng et al., 2009). From a mathematical perspective, on the other hand, the geometry of the confined space may be described with pores beyond the traditional slit pore model (Gor et al., 2012; Jagiello and Olivier, 2009; Landers et al., 2013). In this work, we choose the traditional slit pore for PSD applications and focus on the shortcoming of MFA for calculation of the attractive term of the excess Helmholtz free energy functional. For comparison with existing results, we consider three kinds of NLDFT methods including two WDA and one based on density expansion. For quantitative comparison, the theoretical error bound scale (TEBS) of the numerical PSD solution is proposed to rank the performance of four NLDFT methods to represent the PSDs. For further computer-aided screening of amorphous materials, we introduce a more reliable PSD for understanding the solid heterogeneity. The ultimate goal is to calculate the adsorption capacity of well-characterized activated carbon by high accuracy NLDFT.

Method and experimental validation

Samples and isotherm measurement

The experimental isotherms of sample A, B, and C were from Sweatman and Quirke (2001). Sample D from activated carbon of Huaibei Activated Carbon Co., Ltd for the purpose of mixture gas separation was degassed prior to all gas adsorption measurements (300℃; 3 h; 10⁻⁵kPa). CH₄ and CO₂ adsorption isotherms were obtained at 298 K using Hiden IGA-2 to test the isotherm within the pressure range from 0.03 to 2 MPa.

Potential model

We treat CH₄ and CO₂ as Lennard–Jones (LJ) particles. The energy parameters are ɛ_ff/k = 148.10 and 236.14 K for CH₄ and CO₂, respectively, and the size parameters are σ_ff = 0.381 and 0.364 nm (Maddox and Gubbins, 1997). Here k stands for the Boltzmann constant. As for multiple center character of CO₂ molecule, the potential parameters of CO₂ are determined by fitting the experimental phase diagram when the temperature at 298 K with the modified Benedicit-Webb-Rubin (MBWR) equation of state (Ihara et al., 2012). In mapping the CO₂ properties into those for a LJ system, the cutoff distance is set at 3.5σ_ff.

We follow the slit pore model as in standard materials characterization to represent the pore geometry. The Steele 10-4-3 wall potential is used to describe interaction between a gas molecule and each surface (Steele, 1973). The parameters for the solid wall are ɛ_ss/k = 28.0 K, σ_ss = 0.34 nm, number density ρ_s = 114 nm⁻³, and graphite spacing interlayer Δ = 0.335 nm (Cao et al., 2002). The size and energy parameters between solid and fluid are calculated through the standard Lorentz–Berthelot rules.

NLDFT methods

We have tested four kinds of NLDFT methods, symbolized as Mean Field Approximation (MFA), First-order Mean Spherical Approximation (FMSA), WDA-Y, and WDA-L. The main differences of above methodologies lie in the functional form of attractive term in excess Helmholtz free energy ( ${\text{F}}_{\text{att}}^{\text{ex}}$ ), including,

MFA

$${\text{F}}_{\text{att}}^{\text{ex}}[\rho (\mathbf{z})]=\frac{1}{2}\int \int \rho (\mathbf{z})\rho (\mathbf{z}\mathbf{\text{'}}){\mu }_{\text{att}}(|\mathbf{z}-\mathbf{z}\mathbf{\text{'}}|)\text{d}\mathbf{z}\text{d}\mathbf{z}\mathbf{\text{'}}$$

(1)

FMSA

$$\beta {\text{F}}_{\text{att}}^{\text{ex}}[\rho (\mathbf{z})]=\beta {\text{F}}_{\text{att}}^{\text{ex}}[{\rho }_{\text{b}}]-\int C_{\text{att}}^{(1)}(\rho (\mathbf{z})-{\rho }_{\text{b}})\text{d}\mathbf{z}-\frac{1}{2}\int \int C_{\text{att}}^{(2)}(|\mathbf{z}-\mathbf{z}\mathbf{\text{'}}|)(\rho (\mathbf{z})-{\rho }_{\text{b}})(\rho (\mathbf{z})-{\rho }_b)\text{d}\mathbf{z}\text{d}\mathbf{z}\mathbf{\text{'}}$$

(2)

$$C_{\text{att}}^{(1)}=-{\mu }_{\text{att}}^{\text{ex}}=-\delta \beta {\text{F}}_{\text{att}}^{\text{ex}}/\delta \rho (\mathbf{z})$$

(3)

$$C_{\text{att}}^{(2)}(|\mathbf{z}-\mathbf{z}\mathbf{\text{'}}|)=-{\delta }^2\beta {\text{F}}_{\text{att}}^{\text{ex}}/\delta \rho (\mathbf{z})\delta \rho (\mathbf{z}\mathbf{\text{'}})$$

(4)

WDA

$${\text{F}}_{\text{att}}^{\text{ex}}[r(\mathbf{z})]{=F}^{\text{MFA}}[r(\mathbf{z})]{+F}^{\text{cor}}[r(\mathbf{z})]$$

(5) where $\beta =1/\mathrm{kT}$ , F^MFA in equation (5) standing for mean field contribution is the result of equation (1); for WDA-Y, free energy induced by correlation interaction (F^cor) is calculated with the unit of volume; for WDA-L, it is per molecule number. The mathematical details of these methods can be found in the previous publications (Fu et al., 2015a, 2015b). The NLDFT methods allow us to calculate the excess adsorption density at 298 K for CH₄ and CO₂ adsorption in the slit pores with the pore width changing from 0.65 to 5.00 nm (ΔH = 0.05 nm incremental) and pressure from 0.2 to 2.0 MPa (ΔP = 0.2 MPa).

Grand canonical Monte Carlo simulation (GCMC)

We adopt Towhee-7.0.4 for GCMC calculations (Martin, 2013) and compare the simulation results with NLDFT predictions. The total number of molecules in the simulation box is 4000. It takes 10⁶ configurational steps for equilibrium and 10⁶ steps for production in each simulation run. GCMC involves three types of trail moves: insertion, deletion, and displacement of gas molecule in the slit pore.

Numerical solution of the AIE

The AIE provides an explicit connection between the PSD function and the adsorption isotherm of a realistic porous material

$$N(P_j)={\int }_{Hmin}^{Hmax}\mathbf{\theta }(P_j,H_i)f(H_i)\text{d}H(j=1,2,3,\dots \text{t};i=1,2,3\dots \text{m})$$

(6) where N is the experimental adsorption isotherm, which is smoothed by a linear interpolation according to the theoretical pressure; θ is the kernel matrix representing the excess density of a specific gas in a series of slit pores at certain pressure condition; $f(H)$ is the PSD function. In this work, we adopt the B-spline (one-dimensional and four degrees) function as the PSD function. It is to combine the experiment isotherm and the simulated isotherm by a trail-and-error procedure following ${\rho }_{\text{z}}{\int }_{H_{min}}^{H_{max}}f(H)\text{d}H=1$ , where ${\rho }_z$ is the true density of the characterized samples (Vega, 2007). Parameters are optimized by the objective function (equation (7)) with a simulated annealing method from MATLAB R2014a. For qualitative comparison of the characterized samples (Li et al., 2014), the total pore volume must satisfy $V_{\text{psd}}={\int }_{H_{min}}^{H_{max}}f(H)\text{d}H$ , and the surface area is $S_{\text{psd}}=2{\int }_{{\text{H}}_{\text{min}}}^{{\text{H}}_{\text{max}}}f(H)/H\text{d}H$ . The PSD is obtained by minimization of the objective function

$$F_{\text{obj}}=\sum _{j=1}^{\text{t}}[N^{\exp }(P_j)-N^{\text{cal}}(P_j){]}^2(j=1,2,3\dots \text{t})$$

(7) where $N^{\text{cal}}(P)$ stands for the fitting isotherm, and $N^{\exp }(P)$ is the adsorption isotherm obtained from experiments.

TEBS

AIE can also be simplified into linear least square problem, $\mathrm{Af}=N$ , where f is the numerical solution of PSD and A is calculated by $A=\theta \times \Delta H$ , where ΔH is the step length in changing the pore width.

The residue of equation (7) is the main part of forward error control by optimization algorithm. For such ill-posed problem, the error bound of numerical solution from backward error analysis can be calculated from

$$\frac{1}{\mathrm{cond}(\mathbf{A})}·\frac{\parallel \mathbf{r}\parallel }{\parallel \mathbf{N}\parallel }\le \frac{\parallel {\mathbf{f}}^{\sim }-\mathbf{f}\parallel }{\parallel \mathbf{f}\parallel }\le \mathrm{cond}(\mathbf{A})·\frac{\parallel \mathbf{r}\parallel }{\parallel \mathbf{N}\parallel }$$

(8) where $\mathrm{cond}(\mathbf{A})$ is the condition number of matrix A , which can be calculated numerically by the ratio of the biggest and the smallest singular value by singular value decomposition of matrix A ; r is the residue vector, $\mathbf{r}=\mathbf{Af}-\mathbf{N}$ ; and $\stackrel{\sim }{\mathbf{f}}$ is the analytical solution of PSD. It is obvious that the error bound is strongly dependent on the condition number of the matrix A , therefore, we define the TEBS by the right part of equation (8) dividing the left part, viz. the square of the condition number of the matrix A , through which the effect of optimization algorithm is consequently excluded for elimination of r .

Multimodels calculation of adsorption capacity

The AIE suggests that, ideally, the PSD function acquired from one gas isotherm is applicable for prediction of the adsorption isotherms of other gases. Figure 1 shows schematically the numerical procedure. Based on the PSD function, we can also predict the adsorption isotherm for mixtures.

Figure 1.

Schematic of parameter fitting and prediction route. Red solid line: isotherm fitting; red dash line: pure gas isotherm prediction. Violet line: binary gas isotherm prediction with data from experimental isotherm (solid line) and predicted isotherm (dash line).

The conventional method to predict the adsorption isotherm for a binary gas mixture is based on those for corresponding pure gases within the framework of the ideal adsorbed solution theory (IAST) (Do, 1998). Because the numerical performance of IAST has been well established, it provides a convenient benchmark for calibration of our theoretical predictions. The average deviation between the two methods can be calculated from

$$\text{D}=\frac{1}{\text{t}}\sum _{j=1}^{\text{t}}\frac{|N^{\exp }(P_j)-N^{\text{pre}}(P_j)|}{N^{\exp }(P_j)}\times 100\%(j=1,2,3\dots \text{t})$$

(9)

For binary systems considered in this work, we calculate $N^{\text{exp}}(P)$ based on isotherms from experiment, and $N^{\text{pre}}(P)$ from one experimental adsorption isotherms and one predicted isotherm.

With the experimental and predicted isotherms, computation for adsorption capacity of adsorbent in application of pressure swing adsorption (PSA) can also be made through the published model (Maring and Webley, 2013) in large scale. The working capacity (WC) of the model is defined as the mole concentration of CO₂ in the product mixture and the specific work (SW) is calculated from equation (10), where $W_{\text{vacuum}}$ and $W_{\text{blower}}$ are the work consumed by vacuum and the blower

$$\mathrm{SW}=\frac{W_{\text{vacuum}}+W_{\text{blower}}}{\mathrm{WC}}$$

(10)

TEBS comparison

In previous study (Fu et al., 2015a), results in single slit pore of four NLDFT methodologies were all close to the calculation of molecular simulation, for example, the average deviation between MFA and GCMC within reduced bulk density range from 0.002 to 0.5 is merely 0.82%. For further comparison, we calculate the excess density of CH₄ and CO₂ in series of slit pore at different pressure conditions when the temperature is 298 K to let the inconspicuous differences among methodologies accumulate as the calculation increases. The kernel matrix for each gas can be built consequently.

The numerical solution of PSD is affected by both optimization algorithm in solving AIE for isotherm fitting and theoretical methodologies like NLDFT and GCMC. With the effect of numerical algorithm excluded, the TEBS is used for quantitative comparison to rank the performance of four NLDFT methods with the same potential parameters in calculation for the same size of kernel. The narrower the TEBS is, the more accurate the solution of AIE will be. For a well-posed problem with TEBS equal to 1, such quality can be treated as a baseline for measurement of the distance between numerical and analytical solutions of PSD.

Figure 2 shows the TEBS of PSD with the kernel of both CH₄ and CO₂. The results from WDA-Y predictions are the narrowest among four theoretical methodologies, nearly 10⁴ times narrower than MFA and very close to that from molecular simulation. Since WDA-Y is the modification of MFA in calculating the attractive term of the excess Helmholtz free energy functional with the correlational factor into consideration (Yu, 2009), which means such method is likely to provide more reliable PSD for amorphous material characterization.

Figure 2.

TEBS comparison of methodologies. Kernel of gases: CH₄, blue column; CO₂, yellow column. TEBS: theoretical error bound scale.

Reasons to account for the distinction between different adsorbates are complex. However, for those methods with TEBS less than 10¹², FMSA, WDA-Y, WDA-L, and GCMC, it is obvious that the TEBS of CH₄ kernel is a little higher than that for CO₂. With the influence of temperature taking into consideration, 298 K is higher than the critical temperature of CH₄, while lower than that for CO₂. According to a previously published work by Liu et al. (2011), it is random distribution for fluid molecule in supercritical condition and layer-by-layer structure of those subcritical fluids in the slit pore. We conclude it is helpful to adopt subcritical fluid for the calculation accuracy improvement of excess adsorption in the single slit pore and its TEBS value can be reduced consequently. Therefore, in the following calculations the PSD functions are all based on the kernel of CO₂.

PSD from isotherm fitting

The CO₂ isotherm of four activated carbons is fitted with the kernel calculated by WDA-Y. As shown in Figure 3, the fitting results are in good accordance with the experiment data of four samples, through which the required PSD is also validated for the minimum fitting deviation.

Figure 3.

Fitting of CO₂ isotherm. Experiment isotherm of samples: ○, A; □, B; , C; ⋄, D. Calculation result: −, A; -·-, B; –, C; ċċċ, D.

The form of PSD listed in Figure 4 is varied from the reference since different distribution functions. The advantages of such method include no numerical boundary, smooth and nonnegative constraints. Moreover, since no prior function form, the distributed value of the single slit pore is independent of each other.

Figure 4.

PSD of samples. Samples: •, A; ▪, B; , C; ♦, D. PSD: pore size distribution.

Adsorption capacity prediction

Based on the PSD function, we are able to make a qualitative comparison of the structural parameters. Table 1 lists those structural properties of the samples calculated from the PSD functions and those from previous calculations. It is remarkable that the pore volume of activated carbon is positive related with its surface area and adsorption capacity. From such comparison, the adsorption performance of each sample can be predicted qualitatively.

Table 1.

Structural properties analysis of samples.

Samples	Pore volume Vpsd/cm3g−1	Pore volume (Sweatman and Quirke, 2001) Vref/cm3g−1	Surface area Spsd/m2 g−1
A	2.17	1.58	2018.9
B	0.97	0.72	1056.9
C	1.03	0.65	1022.0
D	0.59	–	665.7

For further application, Table 2 shows the CH₄ isotherms for each sample predicted based on corresponding PSD function obtained from fitting CO₂ isotherms. It is noted that the prediction is very close to experimental isotherm with average prediction deviation at most 13%. Moreover, it shows the corresponding binary adsorption using the IAST. The excellent agreement means that it is alternatively as a predictive tool with the required PSD.

Table 2.

Average deviation in prediction.

Samples	Dp / %	Dbi *  / %
A	5.28	2.51
B	9.46	6.48
C	13.33	9.07
D	8.42	6.15

a

Mole fraction ratio of binary gases, $x_{{\text{CO}}_2}:x_{{\text{CH}}_4}=0.88:0.12$ .

In Figure 5, the adsorption performances of sample D in PSA process are predicted with the adsorbent quality hypothesizes of weight 1 kg, packing density 617.6 kgm⁻³, porosity 0.69, and specific heat 1050 J(kg K)⁻¹ at the condition of feed gas pressure 0.1 MPa, adsorption pressure 0.12 MPa, specific heat ratio 1.4, isentropic efficiency 75%, and temperature of gas feeding 298 K. Apparently, the results predicted from the PSD function are consistent with the experiment data. It means half of experimental research cost in previous can be reduced by the calculation based on the PSD of activated carbon with WDA-Y, the accuracy and stability of which is close to GCMC for TEBS comparison while more efficient. Therefore, such procedure is promising for adsorbent screening in large scales.

Figure 5.

Multiscale adsorption performance prediction of sample D. (Mole fraction ratio of binary gases, x(CO₂):x(CH₄) $=0.88:0.12$ ). Symbol: calculated by two experimental isotherm. Solid line: predicted by one experimental and one predicted isotherm.

Conclusion

First, a novel routine is proposed to make computation for adsorption capacity with the PSD of activated carbon in application of PSA. It is to make use of the characterized PSD for the macroscopic quality prediction of activated carbon.

Second, a term, TEBS, is adopted for theoretical methodologies comparison. The TEBS is better for its quantitative ranking functional. Compared with MFA, it is remarkable that WDA-Y version of NLDFT with the narrowest TEBS of PSD is better for reliable PSD of amorphous materials.

With such methodology, it is not only for qualitative characterization, but also promising in computer-aided screening of activated carbon materials.