CET Volume 86 DOI: 10.3303/CET2186186 Paper Received: 23 September 2020; Revised: 24 February 2021; Accepted: 2 May 2021 Please cite this article as: Moliner C., Antonucci B., Focacci S., Maclean Heap J., Moreno Martel A., Hamzah F., Martin C., Martinez-Felipe A., 2021, Molecular Dynamics Simulation of the Interactions Between Carbon Dioxide and a Natural-based Carbonaceous Microporous Material, Chemical Engineering Transactions, 86, 1111-1116 DOI:10.3303/CET2186186 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 Molecular Dynamics Simulation of the Interactions between Carbon Dioxide and a Natural-Based Carbonaceous Microporous Material Cristina Molinera*, Beatrice Antonuccia,b, Simona Focaccia,b, Jonathan Maclean Heapb, Aldo Moreno Martelb, Fazlena Hamzahc, Claudia F. Martínb, Alfonso Martinez-Felipeb a Department of Civil, Chemical and Environmental Engineering, University of Genoa, Via Opera Pia 15, 16145 Genoa, Italy b Chemical and Materials Engineering Group, School of Engineering, University of Aberdeen. King’s College, Aberdeen AB24 3UE, UK c Faculty of Chemical Engineering, Universiti Teknologi MARA, UiTM 40450 Shah Alam, Selangor, Malaysia cristina.moliner@edu.unige.it Oil palm plantation has drastically changed the scenario of Malaysian agriculture and economy since the second half of the 20th century. World palm oil production in 1980 was around 5.0 million tonnes, doubled to 11.0 million tonnes in 1990, reaching at the beginning of the 21th century approximatively 21.8 million tonnes per year. Malaysia is responsible of about half of the world palm oil production (10.8 million tonnes), thus becoming the world’s largest producer and exporter of palm oil. Oil palm industries also generate large amounts of lignocellulosic biomass waste as a sub-product, which currently has no economic market value other than feedstock for energy valorisation in some limited cases. The utilization and reformulation of this biomass as new materials with added value, would reduce accumulation of waste, mitigating the environmental impact of such intensive industry sector, and would simultaneously generate further economic revenue. In this work, we have studied the interactions of CO2/N2 streams with one slit-like graphite layers, using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics simulations, to simulate new carbon-based adsorbents from empty fruit bunches (EFB). Simulations were conducted using the high-performance computing (HPC) service provided by the University of Aberdeen (Scotland, UK), following accuracy tests to yield minimum operational costs. The simulations provided density profiles of CO2 and N2 molecules as a function of temperature (300/400/500 K), pressure (1/5/10 atm) and micro-meso porosity (12/20/30 Å). The results showed that the number of CO2 molecules retained in the graphite framework decreases at higher temperatures and lower pressures, in accordance with the exothermicity of the process, and smaller pore sizes promote CO2 adsorption at the tested pressures, due to the superposition of the Van der Waals force given by two adjacent walls of slit pores. Our findings exhibit the potential of lignocellusic materials from palm oil tree as materials for CO2 capture. 1. Introduction Oil palm plantation has drastically changed the scenario of Malaysian agriculture and economy since the second half of the 20th century. World palm oil production in 1980 was around 5.0 million tonnes, doubled to 11.0 million tonnes in 1990, reaching at the beginning of the 21th century approximatively 21.8 million tonnes per year. Malaysia is responsible of about half of the world palm oil production (10.8 million tonnes), thus becoming the world’s largest producer and exporter of palm oil (Zaki and Rahim, 2019). Oil palm industries also generate large amounts of lignocellulosic biomass waste as a sub-product, which currently has no economic market value other than feedstock for energy valorisation in some limited cases. GHGs are also generated at various stages of production (transportation of raw material, combustion of fuel in boilers or main oil production steps) (Hosseini and Wahid, 2015). Carbon, Capture and Storage (CCS) technologies are in constant development to mitigate GHGs emission. The main existing technologies include, among other, 1111 physical and chemical adsorption, with promising adsorbents for CO2 capture as activated carbons (ACs) and zeolites (Pardakhti et al., 2019). ACs are promising adsorbents for carbon capture due to their (i) abundant availability; (ii) low-cost and sustainability; (iii) high surface area and possibility of tailoring number and size of pores during the activation process and (iv) high adsorptive and high selectivity characteristics. The maximum adsorptive capacity of ACs largely depends on the raw materials’ structure and its production, by dehydration and carbonization of the biomass followed by chemical or physical activation step (Yahya et al. 2018). Several experimental studies have dealt with the adsorption of single CO2 or CO2 mixtures on carbon-based material such as ACs and graphite (Hinkov et al. 2016). Fewer computational results are reported on the selective adsorption of pure CO2 or mixtures on ACs. Monte Carlo (MC) simulations we used to predict the thermodynamic equilibrium properties of the CO2−carbon pore system to evaluate the influence of surface functional groups of carbon pore surfaces on the adsorption of CO2 (Liu and Wilcox, 2012a). The effect of pore width and surface charges on the adsorption, diffusion and separation of CO2 in a CO2/H2 gas mixture was studied through Molecular Dynamics (MD) simulations (Trinh et al. 2015). The results implied that surface charges can be used to favor the enrichment process of CO2. Cao and Wu (2005) studied the separation of H2 and CO2 via grand canonical MC simulations. H2 showed no adsorption preference over CO2 at higher temperatures as the properties of the two gases are very similar whereas the activated carbon exhibited strong adsorption preference to CO2 at room temperature. The effect of different pore geometries on the selectivity of carbon materials for the adsorption of CO2 from mixtures of CO2/H2 was evaluated by Kumar and Rodriguez-Reinoso (2012) through molecular simulations. Simulation results showed a high sensitivity of the selectivity for CO2 to the pore structure and fluid composition on the interactions of CO2 with the pore walls. Mixtures containing more components have been also studied. The adsorption isotherms and selectivity of CO2 from CO2/CH4, CO2/N2, and CO2/H2O mixtures on oxygen-containing surfaces were calculated using MC simulations (Liu and Wilcox, 2012b). The induced polarity of the surface functionalization increased the selectivity of CO2 over CH4 and N2 whereas water vapor was preferentially adsorbed over CO2. Few studies have been done for the mixture CO2/N2. MC simulations showed the high dependence of the equilibrium selectivity of CO2 and N2 in nano-porous carbon membranes (Jia et al. 2007) on pore size and mixture composition. MD simulations were applied to study the separation of CO2 and N2 through bilayer graphene oxide (GO) membranes (Wang et al. 2017) and the results suggested that the separation performance could be optimized through regulating the microscopic structure of the GO membranes. In this work, the effect of temperature, pressure and pore width on the adsorption of an equimolar CO2/N2 mixture on a carbon-based adsorbent have been investigated through molecular dynamics (MD) simulations. The paper is organized as follows: In Section 2 we present the description of the potential model, initial configurations set-up and simulation details. Results are shown and discussed in Section 3 and general conclusions are presented in Section 4. This preliminary model will provide basic knowledge for the design of optimised carbon-based adsorbents obtained from oil palm residues and will set the basis for further functionalisation and optimisation. 2. Methodology 2.1 Theoretical background Classical MD simulations were performed using LAMMPS (Plimton et al. 2007), to generate the dynamical trajectories of a system of interacting particles by integrating Newton's equations of motion: = ; restated = ( ) (1) where is the acceleration of each particle; is the force; is the mass of each particle; is the position of each particle; is the time; and is the potential energy. The discretization of the equations in time is done through the standard Velocity Verlet Algorithm. This algorithm updates after each timestep the position of each atom with the following equation: ( + ) = ( ) + ( ) + 12 ( ) (2) where ( + ) = ( ) + 12 ( ( ) + ( + )) (3) ( ) = ( ) (4) 1112 Force field and electrostatic forces Physisorption processes of CO2 and N2 in microporous carbons are predominantly associated with Van der Waal's forces (also known as dispersion-repulsion forces) and electrostatic forces (also known as Coulombic interactions). There are several previous studies in the literature that report on the MD Simulation of the CO2 and N2 adsorption phenomena on graphite through the 12-6 Lennard-Jones (LJ) potential and the long-range Coulombic interactions formulas (Liu and Wilcox, 2012b). The inter-particle potential energy ( ) is represented through the standard 12-6 LJ potential: = 4 − < (5) where is the depth of the potential well which indicates the strength of the interaction; is the finite distance at which the inter-particle potential is zero; is the distance between two atoms; and is the cut-off distance after which the effect of these forces becomes negligible, equal to 15 Å. The CO2 molecule is modelled using the TraPPE model (Potoff and Siepmann, 2001), which is a rigid three- site model that accounts for the intrinsic quadrupole moment of CO2 using a partial charge at each site. The N2 molecule is represented as a three-site model with two sites located at two N atoms with the third located at its centre of mass (COM) (Kandagal et al. 2017). The standard 12-6 LJ potential is chosen to describe the intermolecular interactions in the CO2/graphene and CO2/N2/graphene systems. Table 1 shows the force field parameter for interactions between different types of particles involved in the simulation. Table 1: Force field parameters for particle interactions Parameter Bonds Value Parameter Bonds Value (kcal/mol) C-C 0.0537 Å C-C 3.0500 O-O 0.1569 O-O 2.8000 N-N 0.0715 N-N 3.3100 The harmonic bond style uses the potential = ( − ) (6) where is the distance of separation; is the equilibrium bond distance; is a force constant including a 0.5 factor. The harmonic angle style uses the potential : = ( − ) (7) where is the angle of separation; is the equilibrium value of the angle; is a force constant including a 0.5 factor. Table 2: Force field parameters for distance and angles Bond (kcalmol-1A2) Angle (kcalmol-1A2) C-O 5,000 1.16 O-C-O 500 180 N-ghost atom 5,000 0.55 N-ghost atom 500 180 As for the Coulombic pairwise interaction, that are the electrostatic interactions that arise from an electric field acting upon a second charge at a distance, we employed the electrical force: = (8) where is the electrostatic constant = with = 8.85419 ∗ 10 and and are the charges on the 2 atoms i and j. The partial charges on C and O atoms in CO2 are qC = 0.70e and qO = -0.35e respectively, where ‘e’ is the elementary charge and is equal to 1.6022 * 10−19 C. Each N2 molecule, instead, is assigned a negative charge on each N atom, that is, qN = -0.482e, and a positive charge at the COM site, that is, qCOM = 0.964e. 1113 2.2 Model set-up MD simulations are performed to study the effect of pressure, temperature and pore width in microporous and mesoporous carbon materials. ACs have been modelled as a series of slit-like pores with a graphite wall (Müller and Gubbins, 1998). In this work, the solid sorbent is modelled as two parallel infinite sheets of graphene treated as a rigid body. The equimolar initial molecules of CO2 and N2 were located at the centre of the slit as shown in Figure 1a. The pore is modelled as two graphene sheets of specified distance from each other and the graphite pore size is defined as the distance between the two carbon atoms on the opposite graphite sheets. Simulations were conducted using the high-performance computing (HPC) service provided by the University of Aberdeen (Scotland, UK). Molecular dynamics simulations were performed in the canonical ensemble (NVT) using the LAMMPS MD package (Plimton et al. 2007). The Nose–Hoover thermostat (Martyna and Klein, 1992) was used for maintaining the constant temperature condition. Periodic boundary condition was applied in all three directions and the electrostatic forces were calculated through Coulombic interaction. Newton’s equations of motion were integrated using the velocity Verlet algorithm with a time step of 1 fs. The gaseous mixture consisted of CO2 and N2 with a molar fraction of 0.5. The initial configuration was generated by inserting 20 gas molecules for each component, following accuracy tests to yield minimum operational costs, inside the two graphene sheets. A time step of 1 fs is used to integrate the equations of motion. The initial system was equilibrated by a NVT simulation running for 10 ns, which is followed by a production run for 10 ns, sufficiently long to obtain good statistics and consistent trajectories. During this period, coordinates of gas molecules were collected every 50 fs. The simulations provided density profiles of CO2 and N2 molecules (Figure 1b) as a function of temperature (300/400/500 K), pressure (1/5/10 atm) and micro-meso porosity (12/20/30 Å) (see Section 3). (a) (b) Figure 1: Initial configuration of the molecules inside the pore (a); final configuration of molecules inside the pore and corresponding density profile (Blue molecules: N2; Green molecules: CO2) 3. Results and discussion 3.1 Adsorption of CO2/N2 mixture The local density distribution of CO2 and N2 along the pore width (in the z direction) was calculated for all the temperatures (300/400/500 K), pressures (1/5/10 atm) and pore widths (12/20/30 Å). Figure 2 shows the concentration profile of CO2 (full line) and N2 (dashed line) as a function of the pore width (12 Å) at 300 K and 10 atm. Equivalent profiles were obtained at all tested conditions, not represented for the sake of brevity. Figure 2. Local density distribution of CO2 and N2 (300 K, 10 atm, 12 Å) Two symmetrical peaks can be distinguished denoting a two-layered adsorption, at the same positions for both components (4 and 9 Å) which indicates that adsorption is on the same layer for the components. This fact is in accordance with previous works (Jia et al., 2014) where adsorption took place in two layers when the pore size was above 10 Å. The peaks are located near the pore walls, with a higher and narrower shape for 1114 CO2, which indicates a higher attractive potential of the adsorbent for CO2, whereas molecules of N2 can be observed at the centre of the pore due to a weaker interaction. These profiles indicate that both molecules will compete for adsorption on the pore surface and suggest that charged surfaces could enhance the separation of CO2 over N2. In fact, CO2 could be preferentially attracted by negatively charged surfaces due to its strong permanent quadrupole moment of CO2 and the weaker polarity of N2 compared to that of CO2 (Liu et al. 2018). 3.2 Influence of pressure and temperature on CO2 adsorption capacity The density profiles of adsorbed CO2 along the pore (12 Å) at different operational conditions are shown in Figure 3. In general, the density of CO2 across the entire pore space increases, especially near the pore walls, suggesting an adsorption layer being formed around it. The effect of pressure (Figure 3a) is more evident with respect to that of temperature (Figure 3b) as per the more pronounced different height of peaks. The increase in pressure (from 1 to 10 atm, at 300 K) results in a relevant increase in the concentration of CO2. In contrast, the peaks are higher and narrower with the decrease of temperature (Figure 3b) as expected due to the exothermic nature of the reaction and the reduced thermal motion. Figure 3. (a) CO2 density profiles at varying pressure: 1-5-10 atm (T= 300 K, pore width = 12 Å); (b) CO2 density profiles at varying temperature: 300-400-500 K (P= 10 atm, pore width = 12 Å) Temperature and pressure have no effect on the locations of the peaks, always close to the surface of the carbon-based material. Remarkably, no molecules remained at the centre of the pore at all pressures indicating that, at 300 K, CO2 had mainly diffused to the pore walls. The increase in temperature resulted in a higher quantity of CO2 remaining in the centre of the pore, due to probably weaker attraction forces, in line with previous experimental works (Mishra and Ramaprabhu, 2011). 3.3 Influence of pore size Three different pores sizes were analyzed (12 Å, 20 Å and 30 Å) and are shown in Figure 4 at different tested conditions. Slightly lower adsorption of CO2 was observed with the increasing size of the pore, with a more evident effect at higher temperatures. In general, at low loadings as in this case, the adsorbate molecules are close to the pore walls which are the most energetically-favourable positions as the energy associated with the CO2–CO2 interaction is much smaller than the CO2–surface interaction (Liu and Wilcox, 2012c). The attractive potentials due to the overlap of the pore walls is most significant in the smallest pore, resulting in an increased superposition extent of the van der Waals force, and thus overall densities decrease with increasing pore size. CO2 was found close to the pore walls in all cases, confirming also that surface diffusion has an important contribution to the mass transport of the gas through the carbon membrane regardless the pore size, with the only exception of some remaining molecules in the centre of the pore at higher temperatures, as discussed previously. Figure 4. (a) CO2 density profiles at varying pore width (12/20/30 Å) for (a) 300 K; (b) 400 K and (c) 500 K and 10 atm 1115 4. Conclusions Molecular dynamics simulations were conducted to investigate the adsorption of an equimolar CO2/N2 mixture on a carbon-based adsorbent. The local density distribution of both components along the pore width was obtained for different temperatures (300/400/500 K), pressures (1/5/10 atm) and pore widths (12/20/30 Å). Results have shown that CO2/N2 adsorption occurs in two layers, at the same position for both components indicating a possible competitive reaction. The peaks were located near the pore walls, with higher attractive potential of the adsorbent for CO2. The CO2 adsorption capacity increased at higher pressure and lower temperature as expected, with a more evident effect of the varying pressure on adsorption. None of them had effect on the locations of the adsorption peaks, always close to the surface of the carbon-based material. Slightly lower adsorption of CO2 was observed with the increasing size of the pore, with a more evident effect at higher temperatures, due to the superposition of the Van der Waals force given by two adjacent walls of slit pores. References Cao D.P., Wu J.Z., 2005, Modeling the selectivity of activated carbons for efficient separation of hydrogen and carbon dioxide. Carbon, 43(7):1364-70. Hinkov I., Darkrim Lamari F., Langlois P., Dicko M., Chilev C., Pentchev I., 2016, Carbon dioxide capture by adsorption (review), Journal of Chemical Technology and Metallurgy, 51, 6, 609-626. Hosseini S. E., Wahid M. A., 2015, Pollutant in palm oil production process, Journal of the Air & Waste Management Association, 65,7, 773-781. Jia Y. Wang M., Wu L., 2007, Separation of CO2/N2 Gas Mixture through Carbon Membranes: Monte Carlo Simulation, Separation Science and Technology, 42: 3681–3695. Kandagal V. S., Chen F., Jónsson E., Pringle J. M., Forsyth M., 2017, Molecular simulation study of CO2 and N2 absorption in a phosphonium based organic ionic plastic crystal, The Journal of Chemical Physics, vol. 147. Kumar K. V., Rodríguez-Reinoso F., 2012, Effect of pore structure on the selectivity of carbon materials for the separation of CO2/H2 mixtures: new insights from molecular simulation, RSC Adv., 2, 9671-9678. Liu Y, Wilcox J., 2012a, Effects of Surface Heterogeneity on the Adsorption of CO2 in Microporous Carbons. Environ Sci Technol. 46(3):1940-7. Liu Y., Wilcox J. 2012b, Molecular Simulation Studies of CO2 Adsorption by Carbon Model Compounds for Carbon Capture and Sequestration Applications. Environ Sci Technol. 47(1):95-101. Liu Y., Wilcox J., 2012c, Molecular simulation of CO2 adsorption in micro- and mesoporous carbons with surface heterogeneity. International Journal of Coal Geology 104, 83-95. Liu B., Qi C., Mai T., Zhang J., Zhan K., Zhang Z., He J., 2018, Competitive adsorption and diffusion of CH4/CO2 binary mixture within share organic nanochannels. Journal of Natural Gas Science and Engineering 53, 329- 336. Martyna G.J., Klein M.L., 1992, Nosé–Hoover chains: the canonical ensemble via continuous dynamics. J Chem Phys. 97, 2635. Mishra A. K., Ramaprabhu S., 2011, Carbon dioxide adsorption in graphene sheets, AIP Advances 1, 032152 Müller, E.A. and Gubbins K.E., 1998, Molecular simulation study of hydrophilic and hydrophobic behavior of activated carbon surfaces, Carbon, 36, 10, 1433-1438. Pardakhti M., Jafari T., Tobin Z., Dutta B., Moharreri E., Shemshaki N. S., Suib S., Srivastava R., 2019, Trends in Solid Adsorbent Materials Development for CO2 Capture, Applied materials and interfaces 11, 38, 34533– 34559. Plimpton S, Crozier P, Thompson A., 2007, LAMMPS-large-scale atomic/molecular massively parallel simulator. Sandia National Laboratories. Potoff J.J., Siepmann J.I., 2001, Vapor-liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen. AIChE J; 47(7):1676-82. Trinh T. T., Vlugt T. J.H., Hägg M.-B., Bedeaux D., Kjelstrup S., 2015, Simulation of Pore Width and Pore Charge Effects on Selectivities of CO2 vs. H2 from a Syngas-like Mixture in Carbon Mesopores, Energy Procedia 64, 150–159. Wang P., Li W., Du C., Zheng X., Sun X., Yan Y., Zhang J., 2017, CO2/N2 separation via multilayer nanoslit graphene oxide membranes: Molecular dynamics simulation study, Computational Materials Science 140, 284–289. Yahya M. A., Mansor M. H., Wan Zolkarnaini W. A. A., Rusli N. S., Aminuddin A., Mohamad K., Sabhan F. A. M., Atik A. A. A., Ozair L. N., 2018, A Brief Review on Activated Carbon Derived From Agriculture By-Product. AIP Conference Proceedings 1972, 030023. Zaki H.H.M, Rahim A.R.A., 2019, Palm Oil: Malaysia – EU Trade Issue. Kuala Lumpur: Khazanah Research Institute. 1116