Accelerate the Electrolyte Perturbed-Chain Statistical Associating Fluid Theory-Density Functional Theory Calculation With the Chebyshev Pseudo-Spectral Collocation Method. Part II. Spherical Geometry and Anderson Mixing.

To improve the efficiency of electrolyte perturbed-chain statistical associating fluid theory-density functional theory (ePC-SAFT-DFT) calculation of the confined system, in this work, first, the Chebyshev pseudo-spectral collocation method was extended to the spherical pores. Second, it was combined with the Anderson mixing algorithm to accelerate the iterative process. The results show that the Anderson mixing algorithm can reduce the computation time significantly. Finally, based on the accelerated ePC-SAFT-DFT program, a systematic study of the effects of the temperature, pressure, pore size, and pore shape on the CO2 solubilities in the ionic liquids (ILs) confined inside the silica nanopores was conducted. Based on the simulation results, to obtain high CO2 solubilities in the ILs confined in silica, a better option is to use the silica material with a narrow spherical pore, and the IL-anion should be selected specifically considering that it has a more significant impact on the absorption enhancement effect.

1 Introduction

Mitigating CO2 emission from fossil-fueled power plants as well as from transports has become an urgent and worldwide research topic, in which CO2 separation is often needed (MacDowell et al., 2010; Boot-Handford et al., 2014). Ionic liquids (ILs) are promising absorbents for CO2 separation due to their extremely low vapor pressure, high CO2 solubility, as well as low-energy usage for solvent regeneration (Brennecke and Gurkan, 2010; Ramdin et al., 2012; Zhang et al., 2012). However, the viscosity of pure ILs is relatively high compared with common organic solvents, causing a significant decrease in the mass and heat-transfer rates. Using supported ILs has been proposed as a promising solution for practical applications. This can take the advantage of high gas selectivity in ILs, and also, the high surface area of the supported materials can reduce the impact of high viscosity, improve the gas transfer, and hence increase the absorption rate (Zhang et al., 2006; Ren et al., 2012; Romanos et al., 2014).

Research has been conducted to address the confinement effect on the gas solubility in ILs via experiments and molecular simulations (Baltus et al., 2005; Ilconich et al., 2007; Ilconich et al., 2007; Zhang et al., 2010; Iarikov et al., 2011; Ren et al., 2012; Banu et al., 2013; Shi and Luebke, 2013; Romanos et al., 2014; Santos et al., 2014; Budhathoki et al., 2017; Shen and Hung, 2017). According to previous research, several factors will affect the CO2 solubility inside the confined ILs, for example, temperature, pressure, pore size, and shape of porous materials (Baltus et al., 2005; Zhang et al., 2006; Ilconich et al., 2007; Zhang et al., 2010; Iarikov et al., 2011; Ren et al., 2012; Banu et al., 2013; Shi and Luebke, 2013; Romanos et al., 2014; Santos et al., 2014; Budhathoki et al., 2017; Shen and Hung, 2017). However, to screen a suitable IL, optimizing the structure of supported material and operation conditions by experiment or molecular simulations is time and cost consuming, considering the fact that the huge number (1018) of possible ILs can be synthesized (Wasserscheid and Thomas, 2008), as well as the wide temperature and/or pressure range in applications. Therefore, it is desirable to develop a theoretical model to predict the properties of confined IL–CO2 systems.

The classical density functional theory (DFT) is considered as an efficient theoretical method for studying the confined properties (Tripathi and Chapman, 2005; Qiao et al., 2018; Shen et al., 2013; Xu et al., 2008; Sauer and Gross, 2017; Camacho Vergara et al., 2019). In addition, in our previous work (Ji et al., 2012; Ji et al., 2014; Shen et al., 2015; Ji and Held, 2016; Sun et al., 2019), electrolyte-perturbed-chain statistical associating fluid theory (ePC-SAFT) (Cameretti et al., 2005) has been developed to represent the thermodynamic properties of IL systems. Moreover, the developed ePC-SAFT has been combined with DFT (ePC-SAFT-DFT) to describe the properties of IL and CO2/IL confined in nanopores with acceptable results (Shen et al., 2018). Recently, in order to calculate the properties of the confined ILs with ePC-SAFT-DFT efficiently, the Chebyshev pseudo-spectral collocation method (Yatsyshin et al., 2012; Nold et al., 2017) was implemented to accelerate the ePC-SAFT-DFT calculation (Sun et al., 2021). However, only the slit-shaped and cylindrical pores have been considered previously, while for the spherical cavity, the corresponding method has not been available. In addition, an advanced iteration method is required for replacing the simple Picard iteration to accelerate ePC-SAFT-DFT calculation further. Anderson mixing is an elaborate iteration method that has been used in the work by Mairhofer and Gross (2017) and Shen et al. (2021), showing desirable performance in accelerating DFT computing. Therefore, replacing the Picard iteration with Anderson mixing can be an effective strategy.

In this work, the Chebyshev pseudo-spectral collocation method was extended to the spherical geometry, where an expression of 9-3 Lennard–Jones potential for the spherical cavity was derived. In addition, Anderson mixing was used to accelerate the ePC-SAFT-DFT calculation further. Based on the modified program, the CO2 solubility of ILs confined in the silica nanopores was chosen as the representative to conduct the investigation, considering silica is a promising supporting material for ILs (Shi and Luebke, 2013), and the effects of temperature, pressure, pore structures, and IL-ions were investigated systematically.

2 Theory

2.1 Electrolyte Perturbed-Chain Statistical Associating Fluid Theory-Density Functional Theory

According to DFT, in the presence of a solid surface, the grand potential Ω at grand canonical ensemble is given by the equation: where A is the Helmholtz free energy, is the molecular density of component i at position , is the chemical potential, m is the number of segments in a chain for component i, and is the nonelectrostatic external field acting on the segment of component i. The Helmholtz free energy A in ePC-SAFT-DFT for a system can be expressed as (Qiao et al., 2018; Shen et al., 2018): where is the ideal free energy, , , , and are the excess free energies due to hard-sphere repulsions, chain connectivity, dispersive, and electrostatic interactions, respectively. The model performance has been verified in a previous work (Shen et al., 2018), where the model predictions in the density profiles of ionic fluids in the charged pores were compared with the molecular simulation results with high consistency.

The details of , , and have been described elsewhere ( : Barker and Henderson, 1967; Rosenfeld, 1989; Roth et al., 2002; Yu and Wu, 2002a; Yu and Wu, 2002b; Tripathi and Chapman, 2005; Shen et al., 2013). The is composed of two terms, the Coulomb term ( ) and the Debye–Hückel term . accounts for the short-range electrostatic interaction caused by the inhomogeneous distribution at the nearby region of ions, which is based on the Debye–Hückel (DH) theory. accounts for the Helmholtz free energy caused by the long-range electrostatic interaction due to the unsymmetrical distribution of cation and anion. The expressions of and can be found in Li and Wu (2006) and Shen et al. (2018), respectively. All these terms are involved in the ePC-SAFT-DFT calculations for the IL–CO2 systems.

Minimization of the grand potential with respect to the density profile of component i yields the following Euler–Lagrange equation: where q is the charge of component i, and is the mean electric potential at position , which can be obtained by solving the Poisson’s equation (Li and Wu, 2006). The expressions of functional derivatives are described in Rosenfeld et al. (1997), Gross (2009), Roth (2010), and Sun et al. (2021).

2.2 Evaluation of the Convolution-Like Integrals for Spherical Cavity

The expressions of the functional derivatives are described in Rosenfeld et al. (1997), Gross (2009), Roth (2010), and Sun et al. (2021). All the convolution-like integrals involved in ePC-SAFT-DFT calculation can be classified as: where R represents the weighting distances, which are:

In Eq. 4b, is the Heaviside step function, and is the Dirac delta function. The and in Eq. 4b represent the temperature-independent and -dependent hard sphere segment diameter (Barker and Henderson, 1967), respectively. The is the weighting distance for the Debye–Hückel term, which can be referred to in Shen et al. (2018). The was set as 1.5 in order to be consistent with the bulk PC-SAFT model (Gross and Sadowski, 2001). The convolution integrals in Eq. 7a map to with these integral kernels.

For spherical cavities, and only vary with the radial direction, and the analytic expressions of reads (Groh and Mulder, 2000):

2.3 Evaluation of Mean Electric Potential Distribution Inside Spherical Cavity

In the spherical symmetric distribution, the Poisson’s equation reduces to: where is the total charge at position r.

For a spherical cavity with diameter R, by solving Eq. 6 with the boundary conditions and , the expression for the mean electric potential reads:

According to Eq. 7, the boundary electric potential is required, which corresponds to a specific charge density on the internal surface of a nanopore. This value can be obtained via “trial and error” until the charge of ions in the nanopore is equal in magnitude to that with the opposite sign on the internal surface of the nanopore.ü

2.4 Chebyshev Pseudo-Spectral Collocation Method

The Chebyshev pseudo-spectral collocation method for DFT modeling was developed by Yatsyshin et al. (2012). In the Chebyshev pseudo-spectral collocation method, for one-dimensional DFT calculation with an N-point discretization scheme, the density or the weighted density profile over the whole computation domain is determined from the density or the weighted density at a prescribed set of collocation points {z }, k = 1,2 … , N using the barycentric form (Baltensperger, 2002; Berrut and Trefethen, 2004): where the primes indicate that the first and last terms in the sums are divided by 2. The collocation points (z ) can be obtained by a conformal map from the Chebyshev collocation points (x ) (Baltensperger, 2002; Berrut and Trefethen, 2004):

The integrals associated with the DFT calculation can be evaluated with the Clenshaw–Curtis quadrature (Clenshaw and Curtis, 1960).

2.5 Implementation of Chebyshev Pseudo-Spectral Collocation Method in Electrolyte Perturbed-Chain Statistical Associating Fluid Theory-Density Functional Theory for Spherical Cavity

In the ePC-SAFT-DFT calculation, three domains need to be discretized for interpolation. The first domain is for the density profile, the second one is for the weighted density function profiles in the hard-sphere term, and the third one is for the weighted density function profile in the Debye–Hückel term.

For the spherical cavity with diameter R, we defined a diameter R′ as:

The density profile is considered in the domain (0, R') based on the coordinate system illustrated in Figure 1. In other words, this domain is discretized for interpolating the density profile function. In general, the density profile vibrates dramatically near the wall, which implies that more collocation points are required in these regions compared with the middle of the nanopore. In this work, the conformal map proposed by Bayliss and Turkel was used to map the Chebyshev collocation points to domain (0, R') (Bayliss and Turkel, 1992):

FIGURE 1

Schematic diagram of the domains considered in a spherical cavity.

The and in Eq. 11a read: where:

Equation 11a maps the Chebyshev points in the transformed coordinate into the points that cluster in the physical coordinate near with a density that is large when is large. In this work, the was set as: where is the mean diameter of all components.

Therefore, the collocation points are clustered in the physical coordinate near ( ), where high peaks have occurred in the nearby region. The was set as 2 in this work. Similar approaches were used for the discretization of another two domains, i.e., the domains for weighted density function profiles of hard-sphere term (0, R ) and Debye–Hückel term (0, R ).

2.5.1 Evaluation of Convolution-Like Integral

As pointed out by Yatsyshin et al. (2012), the maps (Eq. 4) involved in DFT calculation can be carried out by matrix–vector products with discrete data based on interpolation (Eq. 8) and Clenshaw–Curtis quadrature. In the spherical cavities, for a map with original space discretized into N points and the image space discretized into N points, the map can be represented by N  × N matrix (Yatsyshin et al., 2012): with where a and b are the lower and upper limits of integral for the ith discrete point, W is the corresponding integral kernel in Eq. 5, and is the qth Chebyshev grid of the ith discrete point in the image space: is the corresponding Clenshaw–Curtis weight, and the values at M Chebyshev grids are used to evaluate each integral. In this work, M was set as 45.

For the singularities in the expression of Eq. 13b (i.e., ), the corresponding elements can be evaluated based on the following expression (Sun et al., 2021):

2.5.2 Evaluation of the Mean Electric Potential

Equation 7 can be represented with matrix–vector products with the discrete data: where and are composed of and q at the position of collocation points, respectively.

The two N × N matrixes and read: where and are composed of: where and are the lower and upper triangular matrices with all elements being unity, respectively. The elements in the (N-1) × N matrix read: with

The elements in (N-1) × N read: with where can be evaluated with Eq. 13c with the a and b defined in Eq. 14h. The singularities in equations Eq. 14e and Eq. 14g can also be evaluated based on Eq. 13d.

and are the linear operators implementing the piecewise integrations between two adjacent collocation points, while and are used to sum up the results of these piecewise integrations for obtaining the final integration results of Eq. 7.

2.6 Solving Euler–Lagrange Equation With Anderson Mixing

In the Chebyshev pseudo-spectral collocation method, the Euler–Lagrange equation (Eq. 3) is transformed into a vector equation that can be solved numerically. In our previous work, the Picard iteration was used. In this work, the Anderson mixing is implemented to solve Eq. 3 iteratively (Anderson, 1965). In the Anderson mixing, the Euler–Lagrange equation can be rewritten as: where k represents the Boltzmann constant, represents the bulk phase density of component i, and the superscript res represents the residual quantities.

Eq. 15a can be solved iteratively by (Anderson, 1965): where S is the relaxing factor, , m was set as 50 in this work, the subscript i represents the ith component, and is determined by:

The constrained optimization can be transformed to unconstrained optimization (Fang and Saad, 2008; Walker and Ni, 2011): where

The successive optimization (Eq. 15d) can be solved efficiently by the updating QR factorization, and the necessary Matlab code is given in Walker (2011). In order to avoid divergence, Ns steps of Picard iterations can be performed at the beginning and then switched to the Anderson mixing procedure.

2.7 General Scheme Combined With Anderson Mixing

A general scheme for calculating the density profile of confined systems with ionic contribution was proposed in Sun et al. (2021), as illustrated in Scheme 1. In a loop of Scheme 1, the Euler–Lagrange equation needs to be solved twice. For the first time, solving the Euler–Lagrange equation in a loop, Ns (steps of Picard iteration performed before the Anderson mixing procedure) can be set between 1,200 and 1,500, while for the second time, due to the initial guess of density profile being not far away from the equilibrium density profile, we do not perform Picard iterations prior to the Anderson mixing procedure. (i.e., Ns = 0).

SCHEME 1

Calculating the density profile of ionic liquid (IL)–CO2 system inside nanopores with specific surface charge using the general scheme combined with the Anderson mixing.

In Scheme 1, T, p, , and x refer to those in the bulk phase, which were obtained with ePC-SAFT. H (or R) and Q are the features of nanopores. represents the initial guess of the density profile, and represents the calculated equilibrium density profile with the boundary electric potential . The charge absorbed per area of nanopore Q can be computed with: where the amount of component i adsorbed per surface area for the slit-like pores can be evaluated by: for the cylindrical pores: and for the spherical cavities:

2.8 Model Ionic Liquid–CO2 Confined in Nanopores

Following ePC-SAFT for the ILs in the bulk systems (Ji et al., 2012), an ionic liquid molecule is composed of one IL cation and one IL anion. Each individual IL ion was modeled as a nonspherical species with repulsion, dispersive attraction, and Coulomb interactions. The 9-3 Lennard–Jones potential was used to represent the nonelectrostatic interaction between silica and fluid. The nanopore was modeled as an infinitely large slit, infinitely long cylinder, or spherical cavity. The 9-3 Lennard–Jones potential for large slit and infinitely long cylinder have already been presented in the literature (Fitzgerald et al., 2003; Siderius and Gelb, 2011; Lee, 2016), while for a spherical cavity with diameter R, the 9-3 Lennard–Jones potential can be obtained from: where r is the distance of the fluid molecule from the center of the spherical cavity, is the solid atom density, is the effective solid-fluid diameter, and is the potential representing the interaction between surface and fluid segment. and can be determined with the Berthelot–Lorentz combining rules: where and are the size and potential parameters of a solid surface, respectively, and is the potential parameters of fluid segment i. Here, the used potential parameters for silica are , , and (Pinilla et al., 2005). The ePC-SAFT parameters for ILs were taken from the previous work (Ji et al., 2014), while those for CO2 were taken from Gross and Sadowski (2001).

For the electroneutral silica nanopore, the amount of cation and anion adsorbed per surface area should be equal. The solubility of CO2 in the confined IL with a neutral surface is defined by: where , according to the charge neutrality condition.

3 Results and discussion

3.1 Efficiency of the General Scheme Combined With Anderson Mixing

Calculating the density profile of [C6mim] [Tf2N]-CO2 confined in electronic neutral silica pore with different structures at 333 K and 16.1 bar was selected as an example here to demonstrate the performance of the general scheme combined with the Anderson mixing, and the calculation efficiency was compared with the general scheme only using the Picard iterations. In the calculation, the width of the slit-shaped pore and the diameters for the cylindrical pore and spherical cavity are 5 nm. We used 120 collocation points to represent the density profile inside the slit-shaped pore (functional derivatives only need to be calculated in half of the collocation points due to the symmetry of the density profile), while 60 collocation points were used for the cylindrical pore and spherical cavity. The relaxing parameters used in all these three cases are 0.001.

The calculations were performed on a computer with an AMD core Ryzen 7 PRO 4750U and X64 processor. The version of the compiler is Matlab 2018b. Figure 2 compares the time required for calculations.

FIGURE 2

Time required for calculating equilibrium density profile (not including the time for calculating the matrices required in the Chebyshev pseudo-spectral collocation method).

As illustrated in Figure 2, using the Anderson mixing can improve the calculation efficiency significantly. The time needed for calculating the matrices required in the Chebyshev pseudo-spectral collocation method is listed in Table 1.

TABLE 1

Time needed for calculating the matrices required in the Chebyshev pseudo-spectral collocation method.

	Slit-shaped pore	Cylindrical pore	Spherical cavity
Time, s	0.19	39.05	0.11

As listed in Table 1, the time required for calculating the matrices for the slit-shaped pore and spherical cavity can be ignored. However, the time for the cylindrical pore is pronounced. The reason is that too much time is used for evaluating the Legendre complete elliptic integral of the first and second kinds. This may be a problem in massive ePC-SAFT-DFT calculations, for example, adjusting the model parameters from experimental data, in which the equilibrium density profile needs to be solved multiple times.

However, for systems at the same temperature and pore diameter (cylindrical pore), the elements of these matrices are the same (Sun et al., 2021). Therefore, for the sake of saving computation time, when modeling systems confined in the cylindrical pore, these matrices can be used repeatedly for the systems at the same temperature and pore diameter.

3.2 Model CO2 Solubilities of Confined Ionic Liquids

Based on the efficient algorithm discussed above, ePC-SAFT-DFT can be used in a wide range of IL systems to obtain results efficiently. In this work, several [Cnmim]-based IL-CO2 systems confined in silica nanopore were selected to investigate the effects of temperature, pressure, IL ions, as well as the size and shape of the pore. As the nanopores of different silica materials have been roughly assumed as slit-like (Li et al., 2005; Yang and Yue, 2007), cylindrical (Kresge et al., 1992; Beck et al., 1992; Maddox et al., 1997), and spherical (Nandiyanto et al., 2009) in the previous theoretical work, in this work, these three pore-shape models were adopted.

3.2.1 General View of CO2 Confined In Silica Nanopores

The calculated density profile of [C6mim] [Tf2N]-CO2 confined in 25 Å slit silica pore at 10 bar and 323.15 K is presented in Figure 3.

FIGURE 3

The calculated density profiles of [C6mim] [Tf2N]-CO2 system confined in 25 Å slit silica pore at 10 bar and 323.15 K.

As pointed out by Ho et al. (2013), two competitive mechanisms affect the solubility in a nanopore: One is the adsorption near the pore surface, and the other is the absorption inside the IL. These can be seen in the density profile of CO2. According to the results illustrated in Figure 3, the peak in the density profile near the surface of nanopores corresponds to the first mechanism, while in the middle of the nanopore, the density profile of CO2 tends to the bulk density, which is the reflection of the second mechanism. The absorption enhancement in the confined ILs is mainly contributed by the first mechanism. In this work, the additional solubility x is used to identify the absorption enhancement:

3.2.2 The Influence of Temperature

In general, with the temperature increase, the CO2 solubility in the bulk ILs will decrease. According to our calculation results, CO2 solubility in [C6mim] [Tf2N] confined in slit-shaped pore also decreases with the increase in temperature. This is consistent with the observation by Mirzaei et al. (2017). Typical examples are shown in Figures 3A and 4A. However, with the increase in temperature, the solubility of CO2 in [C6mim] [Tf2N] confined in SiO2 decreases greater than that of the bulk IL. For example, at 1 bar, when the temperature increase from 298.15 K to 373.15 K, the solubility of CO2 in the bulk [C6mim] [Tf2N] decreases by about 0.02 mol CO2/mol IL. Under the same condition, the calculated solubilities of CO2 in the confined [C6mim] [Tf2N] reduce by about 0.03 mol CO2/mol IL. It indicates that confined ILs may make it easier to desorb CO2.

FIGURE 4

Calculated CO2 solubility (A) and additional solubilities (B) of ILs confined in 25 Å slit-shaped pore at 1 bar and different temperatures.

In order to study the effect of temperature on the absorption enhancement effect, we analyzed the relation between the additional solubility and temperatures. As illustrated in Figures 4B and 5B, the additional solubility decreases with increasing temperature. In order to interpret this, the calculated density profiles of [C6mim][Tf2N]-CO2 confined in 25 Å slit silica pore at 10 bar and different temperatures are illustrated in Figure 6. To quantitatively describe the effect of the first mechanism, the ratios of the CO2 density at the first peak near the pore surface (which is also the maximum density inside the nanopore) to its bulk density in ILs were calculated as illustrated in Figure 7. According to the results shown in Figure 7, the ratio decreases with the increase in the temperature when the pressure keeps constant, and consequently, the additional solubility decreases. The same phenomenon can be observed at other pressures and in other pore structures, as listed in Supplementary Table S1.

FIGURE 5

Calculated CO2 solubility (A) and additional solubilities (B) of ILs confined in 25 Å slit-shaped pore at 10 bar and different temperatures.

FIGURE 6

The calculated density profiles of [C6mim] [Tf2N]-CO2 confined in 25 Å slit silica pore at 10 bar and (A) 298.15 K, (B) 323.15 K, (C) 353.15 K, and (D) 373.15 K.

FIGURE 7

Ratios of the maximum density of CO2 in the confined [C6mim] [Tf2N] to its bulk density in [C6mim] [Tf2N].

3.2.3 The Influence of Pressure

Most commercial ILs are physical absorbents for CO2. In general, the CO2 solubility in ILs will increase significantly with increasing pressure. Figure 8A illustrates the CO2 solubility at 298.15 K and different pressures in the 25 Å slit-shaped pore. From Figure 8B, The CO2 solubilities in confined ILs also increase with the increase in the pressure. In addition, the increase in CO2 solubility in the confined [C6mim] [Tf2N] is more significant than that of the bulk [C6mim] [Tf2N] under the same condition. For instance, when the pressure increases from 1 to 50 bar, the solubility of CO2 in the confined [C6mim] [Tf2N] increases by about 0.78 mol CO2/mol IL, while that in bulk [C6mim] [Tf2N] increases by about 0.70 mol CO2/mol IL.

FIGURE 8

Calculated CO2 solubility (A) and additional solubilities (B) of confined ILs in 25 Å slit-shaped pore at 298.15 K and different pressures.

The additional solubility also increases with increasing pressure, as illustrated in Figure 8B. To have a deep insight into this observation, the density profiles of [C6mim] [Tf2N]-CO2 confined in 25 Å slit silica pore at 298.15 K and different pressures illustrated in Figure 9 were taken as one example for the detailed analysis. According to the results shown in Figure 9, with the increase in pressure, the adsorption of CO2 near the pore surface increases, while the adsorption of IL ions near the pore surface decreases. This can be interpreted that at high pressure, a large amount of CO2 inside the nanopore leads to a stronger competitive adsorption capacity of CO2 in the pore surface than that of IL ions. Consequently, the additional solubility increases with the increase in pressure. The same phenomenon can be observed under other conditions, as listed in Supplementary Table S1.

FIGURE 9

Calculated density profiles of [C6mim] [Tf2N]-CO2 system confined in 25 Å slit silica pore at 298.15 K and (A) 1 bar, (B) 10 bar, (C) 30 bar, and (D) 50 bar.

3.2.4 The Influence of Pore Size and Shape

The CO2 solubilities of confined ILs will also be affected by the pore size and shape. According to the calculation results, the adsorption will be enhanced more significantly in the smaller nanopore, which is consistent with the results of Shi and Luebke (2013). Two typical examples are shown in Figure 10. As illustrated in Figure 10, the solubilities of CO2 increase with the decrease in the slit-shaped pore.

FIGURE 10

Calculated CO2 solubility of [C6mim] [Tf2N] at 298.15 K and 1 bar (A) and 30 bar (B) in slit-shaped pores with different widths (solid symbols, data for confined ILs; open symbols, data for bulk ILs.).

The effect of the shape of nanopores was also investigated. According to our calculation results, for nanopores with the same size (i.e., the width for slit-shaped pore and the diameter for cylindrical pore and spherical cavity), ILs confined in the spherical cavity have the highest CO2 solubilities, while those confined in the slit-shaped pore have the lowest. Typical examples are illustrated in Figure 11. This indicates that using supported material with spherical nanopores may lead to better absorption capacity. The calculated density profiles of the three shapes at 298.15 K and 30 bar are presented in Figure 12. It shows that the pore shape affects the density profile near the surface of the pore significantly. For IL ions, the density at the first peak follows the trend that ρ (slit-shaped) > ρ (cylindrical) > ρ (spherical), while for CO2, the opposite trend can be observed. This is because the inward curved surfaces impede the accumulation of large molecules (e.g., IL ions) near the surface. Therefore, under the same situation, the CO2 solubilities (x) follow the trend that x (spherical) > x (cylindrical) > x (slit-shaped).

FIGURE 11

Calculated CO2 solubility of [C6mim] [Tf2N] confined in different-shaped pores at 298.15 K and 1 and 30 bar (the width of slit-shaped pore and the diameter of cylindrical pore and spherical cavity are all 5 nm).

FIGURE 12

Calculated density profiles [C6mim] [Tf2N]-CO2 confined in 50 Å (A) cylindrical, (B) spherical, and (C) slit-shaped silica pore at 298.15 K and 30 bar.

3.2.5 The Influence of Cation

In most cases, the IL-cations are chain-like substances. Therefore, in this part, the effect of alkyl-chain length was studied. Following this, ePC-SAFT-DFT was used to model the CO2 solubilities of [C4mim][Tf2N], [C6mim] [Tf2N], and [C8mim] [Tf2N] inside silica nanopores at different temperatures, pressures, pore widths, and pore shapes. The calculated results are listed in Supplementary Table S1.

As illustrated in Figures 13A and 14A, for the ILs in the same homologous series, the CO2 solubilities in the confined ILs increase with increasing alkyl-chain length in cation, which is the same as that in the bulk ILs. According to Figure 13B and Figure 14B, at low pressures, the additional solubility also increases with the increase in the alkyl-chain length. However, this is not true at high pressures, as illustrated in Figure 14B. This indicates that the absorption enhancement of ILs with longer alkyl-chain decreases faster with the increasing pressure than the ILs with shorter alkyl-chain.

FIGURE 13

Calculated CO2 solubility (A) and additional solubilities (B) of ILs confined in 25 Å slit-shaped pore at 1 bar and different temperatures. (A) Solid symbols, data for confined ILs; open symbols, data for bulk ILs.

FIGURE 14

Calculated CO2 solubility (A) and additional solubilities (B) of ILs confined in 25 Å slit-shaped pore at 298.15 K and different pressures. (A) Solid symbols, data for confined ILs; open symbols, data for bulk ILs.

As demonstrated in Figure 15, with the increase in pore width, the additional solubility decreases, especially for the IL with longer alkyl-chain length in the IL cation. The additional solubilities of confined ILs in different pore shapes are illustrated in Figure 16. In the spherical cavity, the additional solubilities increase more significantly with increasing alkyl-chain length in IL cation than that in the slit-shaped and cylindrical pores.

FIGURE 15

Calculated CO2 solubility (A) and additional solubilities (B) of ILs confined in slit-shaped pore with different widths at 298.15 K and 1 bar. (A) Solid symbols, data for confined ILs; open symbols, data for bulk ILs.

FIGURE 16

Calculated additional solubilities of ILs confined in different shaped pores at 298.15 K and 1 bar in 50 Å. Solid symbols, data for confined ILs; open symbols, data for bulk ILs.

3.2.6 The Influence of Anion

In the bulk phase, the IL anion affects the CO2 solubilities more significantly than IL cation (Anthony et al., 2005). According to Figure 17, the IL anion also has a more significant effect than IL cation in the absorption enhancement. From 1 to 50 bar, the deviation of the calculated additional solubility between [C4mim] [Tf2N] and [C8mim] [Tf2N] changes from 0.003 to 0.009, while for [C6mim] [Tf2N] and [C6mim] [PF6], the deviation of the calculated additional solubility changes from 0.002 to 0.071 from 1 to 50 bar.

FIGURE 17

Calculated additional solubilities of ILs confined in 25 Å slit-shaped pore at 298.15 K and different pressures.

In Figure 18, the calculated additional solubilities of several confined ILs in different pore shapes are presented. The absorption enhancement of anion is more significant in the spherical cavity as shown in Figure 18. For example, in the slit-shaped pore at 30 bar, the deviation of the calculated additional solubilities between [C4mim] [Tf2N] and [C8mim] [Tf2N] is 0.001, and that between [C6mim] [Tf2N] and [C6mim] [PF6] is 0.007. While in the spherical pore, these two values are 0.007 and 0.086. In addition, for [C6mim] [BF4] and [C6mim] [PF6], the additional solubilities are considerably higher than those of [C6mim] [Tf2N] in the spherical cavity at high pressures. This indicates that changing the IL anion of ILs confined in the spherical cavity may be promising to obtain a large absorption enhancement.

FIGURE 18

Calculated additional solubilities of ILs in different-shaped pores at 298.15 K and 30 bar (the width of slit-shaped pore and the diameter of cylindrical pore and spherical cavity are all 5 nm).

Up to here of this section, the effects of temperature, pressure, IL ions, as well as the size and shape of pores on confined CO2 solubilities were investigated. In order to further improve model performance toward practical applications, the ePC-SAFT-DFT model will be combined with different pore models (pore shape and pore size distribution) with the parameters adjusted from experimental data in our future work based on the newly measured experimental data.

4 Conclusion

In this work, the Chebyshev pseudo-spectral collocation method was combined with the Anderson mixing algorithm to further accelerate the ePC-SAFT-DFT calculation. The results show that the computing time can be significantly reduced with the Anderson mixing algorithm. This makes the ePC-SAFT-DFT model more effective in screening promising ILs. However, calculating matrices used in the Chebyshev pseudo-spectral collocation method for the cylindrical pores requires a certain computing time, while in a massive computation, these matrices can be reused to save computation time.

Using the ePC-SAF-DFT model, the CO2 solubilities of ionic liquid confined in silica nanopores were studied. The results show that the CO2 solubilities of confined ILs are always higher than that in bulk ILs under the same condition, and the absorption enhancement effect will be significantly affected by pressure, pore widths, pore shapes, and IL anion. Based on the simulation results, to obtain high CO2 solubilities in silica-confined ILs, a better option is to use silica material with a narrow spherical pore. In addition, the IL anion should be selected specifically considering that the category of IL anion has a significant impact on the absorption enhancement effect.