Variable-Pressure

Jun 6, 2019 - Introduction ... Finally, the results were analyzed based on the two introduced parameters of .... simulation and before starting of nex...
0 downloads 0 Views 7MB Size
Article Cite This: J. Phys. Chem. C 2019, 123, 15523−15533

pubs.acs.org/JPCC

Distinct Understanding of Constant-Volume/Variable-Pressure Experimental Method on CO2 Capture Using Graphtriyne Membrane through an Atomistic Approach Amin Khorsandi-Langol and Seyed Majid Hashemianzadeh* Molecular Simulation Research Laboratory, Department of Chemistry, Iran University of Science & Technology, Tehran 16844-13114, Iran Downloaded via UNIV OF SOUTHERN INDIANA on July 22, 2019 at 10:40:37 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: In this research, the nonequilibrium thermodynamics of the gas permeation process based on the constantvolume/variable-pressure experimental method was explored with the novel algorithm using atomistic simulation. The hybrid force field and the in-house FORTRAN code were used in the proposed algorithm and the pressure was considered through the exerted spring force within the NVT ensemble. The graphtriyne layers were utilized as a porous membrane for the investigation of N2 and CO2 permeation. Two parameters of free tendency and channel cross-sectional area (CCSA) were introduced to analyze the result of the simulation. The result of the simulation revealed that the effect of the CCSA on N2 permeation decreases as the number of graphtriyne layers increases, whereas the CO2 permeation through the membrane is independent of the CCSA. Also, there is a distinct permeation behavior for CO2 and N2 so that first CO2 is trapped within the graphtriyne layers and then permeation is started. In contrast, trapping of N2 increases with increasing the number of graphtriyne layers in the permeation process. This trend originated from the higher free tendency of CO2 relative to N2, which is confirmed by the time-dependent density parameter. Hence, it is possible to make the membrane permselective by adjusting the pressure and number of graphtriyne layers. that satisfy either permeability or selectivity.10 Recently, experimental and computational studies have demonstrated the potential of these novel materials. Various research studies have been carried out in the CO2 separation field using the graphene-based material, such as examination of nanoporous graphene for the separation of N2 from CO2,14 prediction of CO2 and CH4 gas molecule permeation through subnanometer graphene pores,15 study the effect of oxidation degree, interlayer spacing, and channel length for CO2/N2 separation through graphene oxide,16 investigation of diffusion coefficient of different gases in graphene oxide,10 CO2 capture by polyethylenimine-functionalized graphene oxide,17 and using Boron-functionalized graphene oxide−organic frameworks for highly efficient CO2 capture.18 Because of the unique properties of carbon, it seems, there is no end for the allotropes of this element. Apart from different types of these allotropes, two-dimensional structures similar to graphene were predicted by Baughman et al.,19 which provided a new facility for gas separation. Also, these structures have been synthesized experimentally.20,21 Graphyne is a new carbon allotrope that has been formed by joining hexagonal

1. INTRODUCTION Global warming originated from greenhouse gases is the world crisis that should be considered seriously. The rise in the sea level, climate change, and ocean acidification are some consequences of greenhouse gases in the atmosphere that can impact human societies.1,2 Nevertheless, the human activities are the main source of greenhouse gas emission in the atmosphere such as the use of fossil fuels in electricity production, transportation, and industry.3 Carbon dioxide as the main product of fossil fuel combustion plays a key role in heat-trapping and global warming. Therefore, research in the field of CO2 capture has been increased drastically since the past decades.4,5 Membrane technologies have been proved to be an efficient separation method of CO2 in the viewpoint of the energy cost.6,7 In this regard, a literature survey shows that there is a trade-off between permeability and selectivity to achieve a better separation.8 In other words, having simultaneous high permeability and selectivity properties in a membrane is the main challenge. Nonetheless, among the various membranes, the graphene-based structures such as nanoporous graphene,9 graphene oxide,10 reduced graphene oxide,11 and graphene oxide frameworks12 have shown great potential for gas separation. The single atomic layer structure13 and nanoscale porosity are two distinct characteristics of this class of materials © 2019 American Chemical Society

Received: March 1, 2019 Revised: June 3, 2019 Published: June 6, 2019 15523

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C rings through the triple carbon bond. This structure provides triangular nanopores between hexagonal rings so that the nanopore size depends on the number of the triple carbon bond.19 Therefore, properties such as the adjustable size of nanopores and single atomic layer make graphyne-based material an ideal candidate as a molecular sieve for gas separation. Accordingly, there are several attempts for theoretical studies along with practical construction of graphyne-based structures. Some computational studies include water purification,22,23 hydrogen adsorption and purification,24,25 and post-combustion CO2 capture.26,27 Apart from the few studies in the field of CO2 separation using graphtriyne structures, there is not much attention paid to the nonequilibrium thermodynamic rules that govern the permeation process. Hereupon, in this research, molecular dynamics (MD) simulation was applied for modeling of the experimental set-up of gas separation, relying on the aforementioned challenge. Several innovations were incorporated to simulate the overall membrane separation process according to the common experimental method of constantvolume/variable-pressure used in the permeation measurement. For this purpose, MD simulation was carried out by considering a box divided into three regions including feed reservoir, membrane, and dead volume of the downstream chamber similar to the practical setup commonly used in the experiments. This innovative idea was used to scrutinize the permeation of CO2 and N2 gas through the graphtriyne structure as a membrane. The different pressures (i.e., 10 bar up to 50 bar) under NVT ensemble was accomplished by adjusting the spring force. Moreover, the effect of the different number of graphtriyne layers (up to nine layers) on the permeation of gas molecules was investigated to find the bulk properties of the graphtriyne membrane. With the aid of the in-house FORTRAN code, the overall MD simulation run was split up to sequential shorttime MD (sub-MD) simulations to model the flow of gas molecules under NVT ensemble and nonperiodic conditions. Finally, the results were analyzed based on the two introduced parameters of free tendency and channel cross-sectional area (CSA) that is not well understood in the gas permeation process in the literature.

Table 1. Applied Force Field Parameters in MD Simulations31−35 parameters atomic type

charge (e)

σ (Å)

ε (kcal mol−1)

mass (g mol−1)

CCO2

0.6512

2.7570

0.0559

12.0112

OCO2

−0.3256

3.0330

0.1599

15.9994

NN2

−0.4820

3.3200

0.0723

14.0067

Nmid cp ct1 ct2

0.9640 0.0852 −0.0852 0.0000

0.0000 3.6171 3.6171 3.6171

0.0000 0.1479 0.1479 0.1479

1.0 × 10−20 12.0112 12.0112 12.0112

atom located in the middle of two nitrogen atoms (see Figure 1). The interactions between gas molecules and graphtriyne

Figure 1. Denoted atomic type of force field for graphtriyne, CO2, and N2.

layers were given by the consistent valence force field (CVFF)35 model. The σ, ε, and charge parameters in CVFF were assigned to graphtriyne structures considering interactions of graphtriyne structures and gas molecules (refer to Table 1). Finally, Lorentz−Berthelot combination rules were used for cross term pairwise Lennard-Jones interactions, due to the incorporation of the hybrid force field (AREBO and CVFF). 2.2. Molecular System Preparation. The MD simulation was performed for modeling of the experimental set-up for gas permeation based on the constant-volume/variable-pressure method as shown in Figure 2a. Gas permeability through the mentioned setup is calculated according to eq 136

2. COMPUTATIONAL METHOD 2.1. MD Simulation Details. All calculations have been performed using MD simulation as implemented in the LAMMPS28 suite of package. The MD simulations were carried out under NVT ensemble (constant particle number, volume, and temperature). The Nosé−Hoover thermostat was applied for controlling the temperature at 298 K. The cut-off distance for electrostatic and van der Waals forces were considered as 15 and 16 Å, respectively. The two potential walls were placed at the bottom and top edge of the simulation box (z-direction) to confine gas molecules between the two edges. The multilevel summation method29 was used to calculate long-range electrostatic interactions due to nonperiodic conditions in the z-direction. The adaptive intermolecular reactive empirical bond order (AREBO) force field30 was only employed for simulation of intramolecular interaction of the graphtriyne structure to achieve more accurate results. Additionally, force field parameters for CO2 and N2 were obtained from the literature (more details are in Table 1).31−34 A three-site model was selected for the N2 based on recent reports33 that consist of a massless dummy

P=

i dp y 273.15 × 1010Vl × jjj zzz 760 × AT × ((P0 × 76)/14.7) k dt {

(1)

where P is the gas permeability across the membrane in barrer (1.0 Barrer = 1 × 10−10 cm3(STP)cm/(cm2 s cmHg)), V is the dead volume of the downstream chamber (cm3), l is the membrane thickness (cm), T is the experimental temperature (K), A is the effective membrane area (cm2), P0 is the upstream feed gas pressure (psia), and finally

dp dt

( ) is the

steady-state rate of pressure change gradient in the downstream side (mmHg s−1) (derivation of eq 1 is indicated in the Supporting Information (S1)). Also, the selectivity of two gases in the same conditions is calculated using eq 2. P a= A PB (2) 15524

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C

Figure 2. Schematic of gas separation. (a) Experimental set-up and (b) design of MD simulation based on the related set-up.

Figure 3. Height of the gas reservoir under different pressures (a−e). The presentation of hypothetical spring to exert respective pressure schematically (f).

where PA and PB are the permeability of pure gas molecules of A and B, respectively. As an initial configuration, graphtriyne layers as membranes were placed in the xy-plane that is perpendicular to the zdirection. Therefore, the two regions can be considered in the simulation box. The bottom region of the membrane was considered as the dead volume of the downstream chamber. This region was empty at the onset of simulation. Also, the upper region of the membrane was considered as the gas reservoir similar to the experimental set-up and the gas molecules were loaded inside it (Figure 2b). Since the number of loaded molecules (CO2 and N2) are the same in the reservoir under different pressures, the height of the upper region was adjusted to simulate the various applied pressures as shown in Figure 3a−e. 2.2.1. Pressure Exerting on NVT Conditions. Pressure exerting was accomplished by a spring force f(z) = k(zt − z0), where k is the spring constant (force/distance units), z0 is the equilibrium distance from the tether point (distance units), and zt is the membrane distance (z-component) at any time step (time-dependent variable) from the tether point, while the MD simulations were carried out under NVT ensemble. The lowest layer of the membrane was tethered to the upper side (positive z-direction) of the box by the hypothetical spring. According to Figure 3F, the relaxation position of spring was settled at z0 = 0, so that the initial position of the membrane in the simulation box indicates the stretched conditions of spring.

This means that the membrane (graphtriyne layers) inclines to turn back to the upper side of the box to reach the equilibrium conditions at z0 = 0. The k values as spring constants were calculated in such a way that the membrane can exert specific pressure on molecules in the gas reservoir during the permeation process. Two approaches were considered for spring constant estimation. First, based on the Hooke’s law, the force is proportional to the extension. Hereupon, the distance of the lowest layer of the membrane from the upper side of the simulation box (i.e., Δz) was considered as the extension value at the onset of simulation (Figure S1). The second approach is related to spring laws. The membrane movement toward the upper side of the simulation box (MD simulations carried out under NVT ensemble) can be imagined as the membrane is fixed and the upper side of the box moves toward the membrane. This case is similar to a piston that pushes gas molecules into the membrane (Figure S3). It should be noted that the simulation box is nonperiodic in the z-direction because of applying the potential wall. The potential wall acts as a continuous surface in our simulations. According to this argument, the pressure can be calculated through the hook’s law (see eq 3). F = −k Δz or PA = −k Δz

(3)

where F, k, Δz, P, and A are the force, spring constant, membrane distance from the upper side of the simulation box, pressure and cross-sectional area of the simulation box. The 15525

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C corresponding k values to the given pressure has been listed in Table 2. Also, the pressure drop in the reservoir resulted from Table 2. Calculated Spring Constants Corresponding to Different Pressures (a) 10 bar (b) 20 bar (c) 30 bar (d) 40 bar, and (e) 50 bar under NVT Ensemble spring constant ×10−4, N/Å number of layers

(a)

(b)

(c)

(d)

(e)

1 2 3 4 5 6 7 8 9

1.405 1.394 1.384 1.373 1.363 1.353 1.343 1.333 13.237

5.620 5.535 5.453 5.373 5.295 5.219 5.146 5.074 5.005

12.646 12.361 12.088 11.828 11.578 11.338 11.109 10.888 10.676

22.481 21.811 21.179 20.583 20.019 19.485 18.980 18.499 18.043

35.127 33.827 32.620 31.495 30.446 29.464 28.544 27.680 26.866

Figure 4. Illustration of the mechanism of gas permeation.

at the distance of 100 Å from the membrane at the end of each sub-MD simulation and before starting of next sub-MD simulation (Figure 4). Finally, a virtual reservoir was considered to simulate the constant-volume reservoir (Figure 2a) similar to the experimental set-up. For calculation of pressure changes in the constant-volume reservoir over time

gas permeation through the membrane that can be compensated by the spring force so that the pressure remains constant in the upper region during the MD simulations. 2.2.2. Permeation Algorithm. The increasing number of gas molecules in the bottom region (dead volume) resulted by gas permeation leads to an increase in the dead volume pressure, which can push back the membrane toward the upper side. This phenomenon will corrupt the gas permeation procedure under desired pressure because dead volume pressure will exert excess pressure on the membrane. In addition, the membrane movement toward the upper side of the box causes two other challenges. Besides the increase in the height of the dead volume region, the membrane movement by the spring force toward the upper side of the box decreases the Δz and the force decreases due to decreasing of Δz when the spring constant remains constant during the simulation. Therefore, three strategies were considered to maintain constant pressure (first and second strategies), and fix the height of dead volume (third strategy) similar to the experimental set-up during the simulation. At the first, MD simulations per various pressure and/or layers were divided into sequential sub-MD simulations over 50 000 fs. These gas molecules in the bottom region that follow the (Zlower layer of graphtriyne − Zgas) ≥ 50 Å conditions were deleted after each sub-MD simulation and before starting next sub-MD simulation. The selection of distance of 50 Å and beyond 50 Å away from the lowest layer of graphtriyne ensures no return of any permeated molecules. Additionally, the two potential walls were incorporated at the bottom of the simulation box. The distance of the potential walls was about 2 Å with different σ, ε, and cut-off radius. This trick causes the gas molecules to be trapped between the two potential walls (refer to Figure 4). Therefore, the first strategy helps to prevent the accumulation of molecules in the bottom region during the permeation process. As the second strategy, the inhouse FORTRAN code was incorporated to calculate the Δz and updates spring constant at the end of each sub-MD simulation and before starting of next sub-MD simulation (50 000 fs). Results show that the average of membrane movement during any sub-MD simulation is small so that the force changing is negligible at any sub-MD simulation. For example, the average of membrane movement is equal to 1.3 Å under 50 bar as the maximum studied pressure and consequently maximum possible movement (Figure S2). In the third strategy, the lower side of the box (z-direction) was located

dp dt

( ), the deleted gas molecules were collected in the virtual reservoir. The volume of a virtual reservoir was set equivalent to the volume of 1106 particle at 1.0 bar so that the ideal gas law was met. The sub-MD simulations were continued until 90% of gas molecules permeated through the membrane. Since the gas molecules in the reservoir are limited, the pressure cannot be considered constant at the end of the permeation process. Thus, to achieve reliable accuracy for results under constant pressure,

dp dt

( ) was reported based on the 50% of the

permeation process. According to eq 1, variables of V, l, A, and T are constant and the same for both gases, so permeation (P) is proportional to

dp dt

dp dt

( ). Consequently, ( ) was investigated as

the representative of the permeation parameter.

3. RESULTS AND DISCUSSION With the aid of the simulation box, the effect of pressure and the number of graphtriyne layer on the permeation of CO2 and N2 was investigated. The first question that comes to mind is how many layers of graphtriyne in atomistic scale represents the bulk properties. The cut-off was utilized for answering this question. The average distance of the graphtriyne interlayer based on the AREBO force field under periodic boundary conditions was calculated as 3.46 Å. Nonetheless, the average interlayer distance of graphtriyne in the designed simulation box (Figure 2b) under nonperiodic boundary conditions was estimated as 3.50 Å. Also, the results indicate that the graphtriyne interlayer distance was independent of different pressures. Hence, the number of layers must be nine at least by considering the selected cut-off radius for this work (15 Å) and the graphtriyne interlayer distance (3.5 Å). This means that when one gas molecule is located in the middle of fourth and sixth layers of the membrane, its distance from gas molecules in the reservoir is out of interaction range (out of the cut-off radius). In other words, only the interlayer interaction has been considered for the permeation of the respective molecule. Therefore, nine layers of graphtriyne have been investigated. 15526

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C 3.1. Free Tendency of Gases toward Graphtriyne. The free tendency of gas molecules toward graphtriyne can be one of the important parameters in the permeation process. In other words, free tendency describes how gas behaves in the presence of the membrane. It can be said that the gas behavior toward the membrane is affected by the type and magnitude of the interaction energy between them. In this regard, the chemical potential is appropriate to show the interaction between gas molecules and the membrane. In thermodynamics, the chemical potential of a species is energy that can be absorbed or released due to a change of the number of particles in the given species. Also, the self-diffusivity arises from adsorption effects and intermolecular interactions of the membrane and gas molecules as well as the size and shape of gas molecules. Accordingly, self-diffusivity and isothermal adsorption behavior were selected to clarify the free tendency of related gases. The self-diffusion coefficient is proportional to the slope of the MSD curve vs time according to Einstein’s equation37 that to reduce statistical errors in simulations, the self-diffusivity is defined from ensemble averaging Einstein’s equation over a large enough number of molecules38 (eq 4). ⎯ ⎯→ ⎯ 1 ⎯→ i 1 y Ds = jjj zzz ∑ lim ⟨| rk (t ) − rk (0)|2 ⟩ k 6N { k = 1 t →∞ t

interaction energy between gas molecules and one-layer graphtriyne are indicated in Table 4. Results of interaction energy show that electrostatic contribution to the binding energy plays a crucial role in CO2 free tendency toward the graphtriyne sheet. Also, gas adsorption was conducted under isothermal conditions (298 K) through Grand Canonical Monte Carlo (GCMC) for two gases, individually. As shown in Figure 6b, there is a significant difference between the free tendency of CO2 and N2 at the ambient temperature and pressure. However, the average loading of the N2 approach to that of CO2 at high pressures. 3.2. Effect of Channel Cross-Sectional Area on Permeation. Several open channels can be produced by overlapping the triangular surface area of multiple graphtriyne layers. These channels are considered as a pass way for gas molecules to go through the membrane. In other words, the cross-sectional area (CSA) of channels can influence the permeation of gases through the graphtriyne membrane. Nevertheless, the CSA of open channels can be confined by the sliding of layers over each other. In ideal conditions, the cross-section of an open channel can be a maximum value when the triangular area of each layer is placed alongside each other as shown in Figure 7a,b. At any other arrangement of layers, the different CSA can be formed. Some of these configurations are depicted in Figure 7c−f. To show the variation of the CSA of channels resulted from the sliding of layers over each other, MD simulations were carried out at 298 K and 1.0 bar under periodic boundary conditions. Figure 8 illustrates the CSA distribution of channels for each configuration (i.e., 1 layer, 2 layers, etc.) obtained by an inhouse FORTRAN code numerically. It can be seen that distribution curves shift toward the lower CSA with increasing the number of graphtriyne layers, which confirm the inverse relation of CSA of channels and number of graphtriyne layers. It is clear that the effective cross-sectional area (ECSA) of open channels for passing gas molecules are less than CSA due to van der Waals radius of atoms. According to Lorentz− Berthelot combination rules, the ECSA of open channels in the graphtriyne membrane depends on the type of gases which permeates through the membrane. Therefore, assigning the absolute ECSA to graphtriyne is incorrect and unreasonable. Accordingly, the ECSA was estimated for the vertical and horizontal orientation of gas molecules toward the graphtriyne layer (see Figure 9i). For this purpose, the interaction energy of the gas molecules and graphtriyne as a function of distance was assessed. According to Figure 9j, gas molecules were moved from y1 to y2 point along the a⃗ vector. The distance of the central atom of gas molecules (i.e., the C atom in the CO2 molecule and Nmid atom in the N2 molecule) from the y1 point was considered as the distance parameter in interaction energy calculation. As clearly illustrated in Figure 10a, the Δσ as a measure of ECSA of CO2 that is equal to |σ4 − σ1| is greater than that of N2 molecule (i.e., |σ3 − σ2|) in the vertical orientation. It should be noted that σ is the distance at which the binding energy between the gas molecules and graphtriyne sheet is zero (refer to Figure 9k). Furthermore, the center of mass of the triangular surface is the minimum point of interaction energy for two gas molecules. However, the minimum point of interaction energy for two gas molecules of CO2 relative to that of N2 is indicative of more tendency of CO2 in vertical orientation toward the graphtriyne sheet. In contrast, the size of two gas molecules is effective in the binding energy in the horizontal orientation of molecules

N

⎯→ ⎯ rk (t ),

(4)

⎯→ ⎯ rk (0)

where Ds, N, and are self-diffusivity, the number of molecules, the position of the kth molecule of species i at time t, and the reference position of the kth molecule of species i at t = 0, respectively. Also, the center of the mass of each particle was considered as the position of that particle. Simulations were carried out to show the effect of a number of graphtriyne layers on the diffusion coefficient under 10 bar and 298 K. The data acquisition for MSD calculation was performed in 0.2 ns after equilibration of the system. The results have been summarized in Table 3 and Figure 5. According to Table 3, the Table 3. Effect of the Number of Graphtriyne Layers on the CO2 and N2 Diffusion Coefficient diffusion coefficient (cm2 s−1) number of layers 0 1 2 3 4 5 6 7 8 9

CO2 2.98 4.44 2.60 1.92 1.14 5.77 2.08 2.21 1.95 7.02

× × × × × × × × × ×

10−2 10−3 10−3 10−3 10−3 10−4 10−5 10−5 10−5 10−6

N2 2.03 1.36 1.12 9.51 8.20 6.50 5.53 3.93 5.48 5.07

× × × × × × × × × ×

10−2 10−2 10−2 10−3 10−3 10−3 10−3 10−3 10−3 10−3

diffusion coefficient of pure CO2 and N2 are the same approximately. After placing the first graphtriyne layer in the middle of the simulation box, the difference between the diffusivity of gases becomes considerable. The impact of onelayer graphtriyne on the CO2 diffusivity is more significant than this effect on N2 diffusivity. Numerically, the CO2 diffusion coefficient has been dropped about 85% of its initial value while N2 diffusivity was declined nearly 33% relative to its initial value. Figure 6a shows the effect of a number of layers on the diffusivity of CO2 and N2. Moreover, the results of the 15527

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C

Figure 5. Effect of the number of graphtriyne layers (no layer (0L) up to 9 layers (9L)) on the free tendency of CO2 and N2.

Figure 6. Diffusion coefficient vs the number of graphtriyne layers and adsorption isotherm for (a) the presentation of how the number of layers affects diffusivity of CO2 and N2 and (b) the plot of adsorption isotherm of CO2 and N2 under room temperature, respectively.

used to calculate the MECSA in the N2 permeation process. To calculate the MECSA with appropriate accuracy, the shape of ECSA is preferred to be in circular form. Since the ECSA is elliptical due to inequality of σ2 and σ3, the average of two σ +σ parameters of σ2 and σ3 i. e. , 2 2 3 should be calculated. In this way, the elliptical form of ECSA was converted to the circular form. When the ECSA is circular, the MECSA can be calculated based on the radius of the circle inscribed in an equilateral triangle. The radius of the inscribed circle (r) in an equilateral triangle is presented in eq 5 according to geometry rules. The surface area of the equilateral triangle (Atriangle) can be calculated according to eq 6.

Table 4. Interaction Energy Components between Two Gases and Membrane binding energy (kcal mol−1) total/component of binding energy

CO2

N2

total van der Waals electrostatic

−58.777 −4.567 −54.210

−15.202 −7.169 −8.033

(

toward graphtriyne. Referring to Figure 10b in the horizontal orientation, the positive value of binding energy for CO2 shows that permeation of CO2 is prohibited while N2 can permeate through graphtriyne due to its negative binding energy. Analysis of interaction energy derived from Figure 10 is given in Table 5. Hereupon, it can be concluded that the N2 permeation is feasible under the vertical and horizontal orientation of molecules, unlike CO2 permeation that can only occur in the vertical orientation of molecules toward the graphtriyne sheet. Consequently, the CSA of channels is the most important parameter in N2 permeation based on the lower free tendency of N2 relative to CO2. Another issue is finding the number of layers of graphtriyne in which the permeation of N2 through the channel crosssectional area parameter is blocked due to the sliding of layers over each other and consequently decreasing the CSA. In this way, the minimum value of ECSA (MECSA) was determined based on the vertical orientation of linear molecules rather than horizontal orientation. This basis was selected because the required area for passage of molecules in the vertical orientation is less than horizontal orientation. N2 molecules cannot permeate through the graphtriyne membrane bearing the ECSA equal to or less than the MECSA value. The data related to molecules in the vertical orientation in Table 5 was

r=

3 a 6

A triangle =

)

(5)

3 2 a 4

(6)

where the “a” is the side-length of the equilateral triangle. As shown in Figure 10, the distances of σ2 and σ3 from the center of the mass of the equilateral triangle (rε) are equal to 0.65 and 1.05 Å, respectively. So, the sum of these two numbers represents Δσ (i.e., ECSA). The averaged Δσ was used to convert the ECSA of an elliptical shape to circular form (see eq 7). |rε − σ2| = 0.65| o geometry rule o → Δσ = ECSA = 1.70 ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ rECSA } o |rε − σ3| = 1.05 o ~ = 0.85

(7)

The β value is the vertical distance between the equilateral sides of triangle and center of mass of the equilateral triangle 15528

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C

Figure 7. Various arrangements of CSA. Top (a) and side (b) view of ideal conditions with maximum CSA and other configuration of CSA related to 1 layer (c), 3 layers (d), 6 layers (e), and 9 layers (f) from left to right, respectively.

The GCMC calculations in Section 3.1 showed that there is a significant difference between the free tendency of CO2 and N2 to graphtriyne layers at ambient temperature and pressure. Moreover, MSD calculations revealed that the graphtriyne layers has a considerable effect on the diffusivity of CO2 in comparison with N2 at 10 bar. However, the impact of the free tendency in the permeation process will diminish gradually with increasing the pressure. It means that the free tendency of gas molecules for permeation through the membrane is inversely related to the pressure. This relation can be perceived well in gas permeation through one-layer graphtriyne because the CSA parameter is omitted in a one-layer graphtriyne case. Therefore, as shown in Table 6, one-layer of graphtriyne demonstrates the better selectivity (α) under low pressure against high-pressure conditions. Under 10 bar as the lowest pressure condition in Table 7, for two graphtriyne layers, the permeation value is identical for two gases approximately, due to large enough CSA for permeation of two gases. With increasing the number of graphtriyne layers, the difference between permeation values for two gases become significant, gradually. The effect of the CSA on permeation decreases gradually as the number of graphtriyne layers increases. Although the N2 and CO2 permeation can be limited by decreasing the CSA, the permeation of CO2 is higher than that of N2 due to its higher free tendency toward the membrane. The assessment of timedependent density has been carried out to reveal the useful information about the permeation process that was defined by dividing the number of gas molecules into free space volume between graphtriyne layers. There is a maximum value for the time-dependent density of CO2 at the beginning of the permeation process for all layers as shown in Figure 12a. This means that the releasing of CO2 from the membrane does not occur until the free space between the layers is saturated with CO2. In spite of large CSA of the channels for the two-layer graphtriyne, the trapping of CO2 occurs before permeation

Figure 8. Relationship between the number of graphtriyne layers and CSA of channels. the vertical pink dashed line shows the theoretical prediction of minimum CSA required for passing gas molecules.

(see Figure 11a). By subtraction of rECSA (i.e., 0.85 Å) from β value, the surface area in which the permeation of N2 is prohibited can be obtained. A circle with a radius of 2.63 Å (i.e., green-colored circle in Figure 11b) can be drawn by subtracting the rECSA from the β value. Finally, the MECSA can be calculated by an equilateral triangle (i.e., blue-colored triangle) which circumscribe the green-colored circle as shown in Figure 11b. Consequently, MECSA (Atriangle) was calculated as 35.98 Å2 according to eq 5. By referring to Figure 8, the pink dash line that is settled at 35.98 Å2 is considered as a border for permeation of gasses through the membrane. Therefore, for a CSA of channel less than 35.98 Å2, the permeation of gasses is not possible. In other words, gas permeation through channel cross-sectional area parameters for those configurations having the number of layers equal and/or more than eight layers, is not feasible. 3.3. Estimation of Permeation. The number of graphtriyne layers and operating pressure as two effective parameters on the permeation process have been studied. 15529

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C

Figure 9. (i) Different orientation of gas molecules toward the planar graphtriyne in vertical (a, c) and horizontal (b, d) state. Blue and pink contour surfaces around the N2 and CO2 molecules are drawn based on the van der Waals surface using σN (ra) and σO (rc), receptively, (j) the direction of moving gas molecules from y1 to y2 point and (k) definition of σ as the distance at which the binding energy between the gas molecules and graphtriyne sheet is zero.

Figure 10. Binding energy in two states (a) vertical and (b) horizontal orientations.

Table 5. Calculation of the ECSA (Δσ) and rε in Angstroms for CO2 and N2 Based on Lennard-Jones Potential (12−6) CO2

N2

orientation

σ1

σ4

Δσ = σ4 − σ1



σ2

σ3

Δσ = σ3 − σ2



vertical horizontal

2.59

5.06

2.47

3.48 3.88

2.83 3.43

4.53 4.13

1.70 0.70

3.48 3.68

Under 20 bar, the N2 permeation is higher than that of CO2 for up to four layers of graphtriyne. However, this superiority is temporary and vanishes with increasing the number of layers. The decrease of the CSA of the channel becomes a serious obstacle, especially for N2 with increasing the number of graphtriyne layers. Also, the trapping trend of gases under 20 and 10 bar is similar to each other as can be inferred from Figure 12c,d. However, there is an appreciable difference between the releasing time of CO2 and N2 under 20 bar after the trapping procedure. The higher releasing time of CO2 relative to that of N2 return to the higher free tendency of CO2. Although the increase of pressure will lead to enhancing the permeation of both gases but the N2 permeation is higher due to higher trapping of the CO2 relative to N2 molecules when a smaller number of layers are used. Higher permeation of N2

through the graphtriyne layers. Therefore, it can be inferred that trapping is just attributed to the free tendency of CO2 for permeation. However, unlike the behavior of CO2 in permeation, time-dependent density for N2 permeation is proportional to the number of graphtriyne layers. As shown in Figure 12b, N2 trapping increases through the membrane with the increasing of the number of layers. Therefore, N2 trapping gets enhanced with decreasing of the CSA of channels in the membrane with increasing the number of graphtriyne layers. Furthermore, the results derived from Figure 12 confirmed the estimated MECSA of channels (i.e., 35.98 Å2) obtained using the geometry rule, and also, approved the no permeation of gases for those configurations having the number of layers equal and/or more than eight layers. 15530

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C

Figure 11. Estimation of the minimum of ECSA for gas molecules passing schematically.

Table 6. CO2 and N2 Permeation Values for One-Layer Graphtriyne under Different Pressures pressure (bar) gas parameters

10

30

7.48 × 10

1.07 × 10

dp dt N 2

2.62 × 10−7 2.85

( ) ( )

one-layer of graphtriyne

20

dp dt CO 2

a=

−7

PCO2 PN2

−6

40 −6

50

1.46 × 10

−6

1.57 × 10

1.76 × 10−6

4.85 × 10−7

6.34 × 10−7

9.15 × 10−7

1.07 × 10−6

2.21

2.30

1.72

1.65

Table 7. CO2 and N2 Permeation Values under Different Pressures and Number of Graphtriyne Layers (a) 10 bar, (b) 20 bar, (c) 30 bar, (d) 40 bar, and (e) 50 bar dp dt

( ) × 10 (a)

−7

(bar fs−1)

(b)

(c)

(d)

(e)

number of layers

CO2

N2

CO2

N2

CO2

N2

CO2

N2

CO2

N2

2 3 4 5 6 7 8 9

4.90 4.31 3.26 2.74 2.49 2.01 1.70 1.22

4.75 2.69 2.00 1.73 1.16 0.86 0.12 0.08

6.61 6.12 5.22 4.15 3.64 2.73 1.82 1.60

7.79 6.23 5.59 2.76 2.71 1.79 0.74 0.52

8.48 7.32 5.82 4.61 4.75 3.00 2.22 1.77

10.12 8.12 7.13 3.28 3.06 2.36 0.644 0.26

8.85 7.63 6.70 4.57 4.52 3.03 2.40 1.84

11.34 9.45 8.35 5.23 5.00 3.10 0.56 0.11

9.03 8.45 6.33 5.36 6.06 3.04 2.63 1.85

12.80 11.28 8.76 8.86 7.15 6.25 1.13 0.11

tendency and channel cross-sectional area were introduced to analyze the result of the simulation. The Free tendency was assessed based on the diffusion coefficient of gas molecules under 10 bar at ambient temperature and isothermal adsorption calculation under different pressures at 298 K through Grand Canonical Monte Carlo (GCMC). The result confirmed the higher free tendency of CO2 toward the membrane than that of N2. Also, on estimation of the effective cross-sectional area, it was found that the cross-sectional area (CSA) of channels is the most important parameter in N2 permeation regarding its lower free tendency relative to CO2. Nevertheless, the CSA of open channels can be confined by the sliding of layers over each other so that gas permeation through the graphtriyne membrane bearing 8 layers and beyond 8 layers is not feasible. Consequently, CO2 permeation is higher than that of N2 for those configurations having 8 layers and/or more than 8 layers under different pressure conditions. This CO2 permeation superiority is related to the ineffectiveness of the CSA parameter in the N2 permeation

relative to CO2 can be seen for the graphtriyne with up to four layers under 30 bar, up to six layers under 40 bar, and up to seven layers under 50 bar. The N2 releasing time decreases significantly under 50 bar against 10 bar, whereas the CO2 releasing time remains almost unchanged at each pressure, except for 50 bar as shown in Figure 12e,f.

4. CONCLUSIONS Molecular dynamics simulation was implemented to model the constant-volume/variable-pressure method by applying the nonequilibrium thermodynamic rule for the gas permeation process. Accordingly, CO2 and N2 gas molecule permeation through graphtriyne structures as the membrane was assessed. Permeation was investigated under various conditions such as different pressures (i.e., 10 bar up to 50 bar) and the number of graphtriyne layers (i.e., 1 up to 9). Regarding the cut-off radius, it was understood that the bulk properties could be modeled by nine layers of the graphtriyne structure in the atomistic scale. Two novel and innovative parameters of free 15531

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C

Figure 12. Time-dependent density related to the different number of layers (a, c, e) of CO2 gas molecules and (b, d, f) N2 gas molecules under 10, 20, and 50 bar.

process. The permeation in a lower number of graphtriyne layers is pressure dependent. Since the effect of CSA on the permeation of N2 and CO2 is the same for a lower number of graphtriyne layer, pressure and free tendency play the main role in permeation. Moreover, the distinct permeation behavior of CO2 and N2 gas molecules is such that the CO2 permeation from the membrane does not occur until the free space between the layers is saturated with CO2 gas molecules. The CO2 is trapped without any dependency on the number of graphtriyne layers so that even for the two-layer graphtriyne, the trapping of CO2 occurs before permeation through the graphtriyne layers. In contrast, N2 trapping is proportional to the number of graphtriyne layers and get enhanced with the increasing of the number of layers.



MD simulation (50 000 fs) by the spring force under 50 bar; membrane movement toward the upper side of the simulation box can be imagined as the membrane is fixed and the upper side of the box moves toward the membrane (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected], hashemianzadeh@gmail. com. Phone: (+98) 2177240287. Fax: (+98) 2177491204. ORCID

Seyed Majid Hashemianzadeh: 0000-0003-1480-9014 Notes

ASSOCIATED CONTENT

The authors declare no competing financial interest.

S Supporting Information *



The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.9b01988.

ACKNOWLEDGMENTS The authors are grateful to Dr. Foad Mehri, Dr. Seyed Vahid Mousavi, and Dr. Seyed Mostafa Rahimian for discussions and their friendly support.

Hypothetical spring to exert respective pressure schematically; membrane movement during any sub15532

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533

Article

The Journal of Physical Chemistry C



(21) Johnson, C. A., II; Lu, Y.; Haley, M. M. Carbon Networks Based on Benzocyclynes. 6. Synthesis of Graphyne Substructures via Directed Alkyne Metathesis. Org. Lett. 2007, 9, 3725−3728. (22) Lin, S.; Buehler, M. J. Mechanics and Molecular Filtration Performance of Graphyne Nanoweb Membranes for Selective Water Purification. Nanoscale 2013, 5, 11801−11807. (23) Bartolomei, M.; Carmona-Novillo, E.; Hernández, M. I.; Campos-Martínez, J.; Pirani, F.; Giorgi, G.; Yamashita, K. Penetration Barrier of Water through Graphynes’ Pores: First-Principles Predictions and Force Field Optimization. J. Phys. Chem. Lett. 2014, 5, 751−755. (24) Bartolomei, M.; Carmona-Novillo, E.; Giorgi, G. First Principles Investigation of Hydrogen Physical Adsorption on Graphynes’ Layers. Carbon 2015, 95, 1076−1081. (25) Sang, P.; Zhao, L.; Xu, J.; Shi, Z.; Guo, S.; Yu, Y.; Zhu, H.; Yan, Z.; Guo, W. Excellent Membranes for Hydrogen Purification: Dumbbell-Shaped Porous γ-Graphynes. Int. J. Hydrogen Energy 2017, 42, 5168−5176. (26) Bartolomei, M.; Giorgi, G. A Novel Nanoporous Graphite Based on Graphynes: First-Principles Structure and Carbon Dioxide Preferential Physisorption. ACS Appl. Mater. Interfaces 2016, 8, 27996−28003. (27) Apriliyanto, Y. B.; Faginas Lago, M. N.; Lombardi, A.; Evangelisti, S.; Bartolomei, M.; Leininger, T.; Pirani, F. Nanostructure Selectivity for Molecular Adsorption and Separation: The Case of Graphyne Layers. J. Phys. Chem. C 2018, 122, 16195−16208. (28) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (29) Hardy, D. J.; Wu, Z.; Phillips, J. C.; Stone, J. E.; Skeel, R. D.; Schulten, K. Multilevel Summation Method for Electrostatic Force Evaluation. J. Chem. Theory Comput. 2015, 11, 766−779. (30) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A Reactive Potential for Hydrocarbons with Intermolecular Interactions. J. Chem. Phys. 2000, 112, 6472−6486. (31) Harris, J. G.; Yung, K. H. Carbon Dioxide’s Liquid-Vapor Coexistence Curve And Critical Properties as Predicted by a Simple Molecular Model. J. Phys. Chem. 1995, 99, 12021−12024. (32) Makrodimitris, K.; Papadopoulos, G. K.; Theodorou, D. N. Prediction of Permeation Properties of CO2 and N2 through Silicalite via Molecular Simulations. J. Phys. Chem. B 2001, 105, 777−788. (33) Murthy, C. S.; Singer, K.; Klein, M. L.; McDonald, I. R. Pairwise Additive Effective Potentials for Nitrogen. Mol. Phys. 1980, 41, 1387−1399. (34) Jiang, J.; Sandler, S. I. Separation of CO2 and N2 by Adsorption in C168 Schwarzite: A Combination of Quantum Mechanics and Molecular Simulation Study. J. Am. Chem. Soc 2005, 127, 11989− 11997. (35) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Structure and Energetics of Ligand Binding to Proteins:Escherichia Coli Dihydrofolate ReductaseTrimethoprim, a Drug-Receptor System. Proteins: Struct., Funct., Genet. 1988, 4, 31−47. (36) Czichos, H.; Saito, T.; Smith, L. E. Springer Handbook of Metrology and Testing; Springer Science & Business Media, 2011; pp 371−387. (37) Einstein, A. On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat. Ann. Phys. 1905, 322, 549−560. (38) Malek, K.; Coppens, M.-O. Knudsen Self- and Fickian Diffusion in Rough Nanoporous Media. J. Chem. Phys. 2003, 119, 2801−2811.

REFERENCES

(1) Fine, M.; Tchernov, D. Scleractinian Coral Species Survive and Recover from Decalcification. Science 2007, 315, 1811. (2) Moore, S. K.; Trainer, V. L.; Mantua, N. J.; Parker, M. S.; Laws, E. A.; Backer, L. C.; Fleming, L. E. Impacts of Climate Variability and Future Climate Change on Harmful Algal Blooms and Human Health. Environ. Health 2008, 7, S4. (3) Davis, S. J.; Caldeira, K.; Matthews, H. D. Future CO2 Emissions and Climate Change from Existing Energy Infrastructure. Science 2010, 329, 1330−1333. (4) Vitillo, J. G.; Smit, B.; Gagliardi, L. Introduction: Carbon Capture and Separation. Chem. Rev. 2017, 117, 9521−9523. (5) Huaman, R. N. E.; Lourenco, S. A Review on: CO2 Capture Technology on Fossil Fuel Power Plant. J. Fundam. Renewable Energy Appl. 2015, 05, 1−11. (6) Lin, H.; He, Z.; Sun, Z.; Vu, J.; Ng, A.; Mohammed, M.; Kniep, J.; Merkel, T. C.; Wu, T.; Lambrecht, R. C. CO2-Selective Membranes for Hydrogen Production and CO2 Capture − Part I: Membrane Development. J. Membr. Sci. 2014, 457, 149−161. (7) Ramasubramanian, K.; Zhao, Y.; Winston Ho, W. S. CO2 Capture and H2 Purification: Prospects for CO2 -Selective Membrane Processes. AIChE J. 2013, 59, 1033−1045. (8) Robeson, L. M. The Upper Bound Revisited. J. Membr. Sci. 2008, 320, 390−400. (9) Sun, C.; Wen, B.; Bai, B. Application of Nanoporous Graphene Membranes in Natural Gas Processing: Molecular Simulations of CH4/CO2, CH4/H2S and CH4/N2 Separation. Chem. Eng. Sci. 2015, 138, 616−621. (10) Jiao, S.; Xu, Z. Selective Gas Diffusion in Graphene Oxides Membranes: A Molecular Dynamics Simulations Study. ACS Appl. Mater. Interfaces 2015, 7, 9052−9059. (11) Lin, L. C.; Grossman, J. C. Atomistic Understandings of Reduced Graphene Oxide as an Ultrathin-Film Nanoporous Membrane for Separations. Nat. Commun. 2015, 6, No. 8335. (12) Razmkhah, M.; Moosavi, F.; Mosavian, M. T.; Ahmadpour, A. Tunable Gas Adsorption in Graphene Oxide Framework. Appl. Surf. Sci. 2018, 443, 198−208. (13) Meyer, J. C.; Geim, A. K.; Katsnelson, M. I.; Novoselov, K. S.; Booth, T. J.; Roth, S. The Structure of Suspended Graphene Sheets. Nature 2007, 446, 60−63. (14) Liu, H.; Dai, S.; Jiang, D. E. Insights into CO2/N2 Separation through Nanoporous Graphene from Molecular Dynamics. Nanoscale 2013, 5, 9984. (15) Yuan, Z.; Govind Rajan, A.; Misra, R. P.; Drahushuk, L. W.; Agrawal, K. V.; Strano, M. S.; Blankschtein, D. Mechanism and Prediction of Gas Permeation through Sub-Nanometer Graphene Pores: Comparison of Theory and Simulation. ACS Nano 2017, 11, 7974−7987. (16) Li, W.; Zheng, X.; Dong, Z.; Li, C.; Wang, W.; Yan, Y.; Zhang, J. Molecular Dynamics Simulations of CO2/N2 Separation through Two-Dimensional Graphene Oxide Membranes. J. Phys. Chem. C 2016, 120, 26061−26066. (17) Li, X.; Cheng, Y.; Zhang, H.; Wang, S.; Jiang, Z.; Guo, R.; Wu, H. Efficient CO2 Capture by Functionalized Graphene Oxide Nanosheets as Fillers To Fabricate Multi-Permselective Mixed Matrix Membranes. ACS Appl. Mater. Interfaces 2015, 7, 5528−5537. (18) Haque, E.; Islam, M. M.; Pourazadi, E.; Sarkar, S.; Harris, A. T.; Minett, A. I.; Yanmaz, E.; Alshehri, S. M.; Ide, Y.; Wu, K. C. W.; Kaneti, Y. V. Boron-Functionalized Graphene Oxide-Organic Frameworks for Highly Efficient CO2 Capture. Chem. - Asian J. 2017, 12, 283−288. (19) Baughman, R. H.; Eckhardt, H.; Kertesz, M. Structure-property Predictions for New Planar Forms of Carbon: Layered Phases Containing sp 2 and sp Atoms. J. Chem. Phys. 1987, 87, 6687−6699. (20) Marsden, J. A.; Haley, M. M. Carbon Networks Based on Dehydrobenzoannulenes. 5. Extension of Two-Dimensional Conjugation in Graphdiyne Nanoarchitectures. J. Org. Chem. 2005, 70, 10213−10226. 15533

DOI: 10.1021/acs.jpcc.9b01988 J. Phys. Chem. C 2019, 123, 15523−15533