From Crystal to Adsorption Column: Challenges in Multiscale

Oct 5, 2018 - Multiscale material screening strategies combine molecular simulations and process modelling to identify the best performing adsorbents ...
0 downloads 0 Views 3MB Size
Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Separations

From Crystal to Adsorption Column: Challenges in Multiscale Computational Screening of Materials for Adsorption Separation Processes Amir Hajiahmadi Farmahini, Shreenath Krishnamurthy, Daniel Friedrich, Stefano Brandani, and Lev Sarkisov Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b03065 • Publication Date (Web): 05 Oct 2018 Downloaded from http://pubs.acs.org on October 9, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

From Crystal to Adsorption Column: Challenges in Multiscale Computational Screening of Materials for Adsorption Separation Processes Amir H. Farmahinia, Shreenath Krishnamurthya, Daniel Friedrichb, Stefano Brandania, and Lev Sarkisova* a

School of Engineering, Institute of Materials and Processes, The University of Edinburgh, Sanderson Building, EH9 3FB

b

School of Engineering, Institute for Energy Systems, The University of Edinburgh, Sanderson Building, EH9 3FB

Abstract Multiscale material screening strategies combine molecular simulations and process modelling to identify the best performing adsorbents for a particular application, such as carbon capture. The idea to go from the properties of a single crystal to the prediction of the material performance in a real process is both powerful and appealing, however it is yet to be established how to implement it consistently. In this article, we focus on the challenges associated with the interface between molecular and process levels of description. We explore how predictions of the material performance in the actual process depend on the accuracy of molecular simulations, on the procedures to feed the equilibrium adsorption data into the process simulator and on the structural characteristics of the pellets, which are not available from molecular simulations and should be treated as optimization parameters. The presented analysis paves the way for more consistent and robust multiscale material screening strategies.

1. INTRODUCTION A significant number of recent studies focused on the systematic search for the most energy efficient technology to remove carbon dioxide from typical industrial streams. Adsorption technology, based on Vacuum Swing or Pressure Swing Adsorption (VSA/PSA), is considered as one of the promising options1. To a significant extent the surge of interest in adsorption (and also membrane) separations is fuelled by the unprecedented developments in material science, where several new classes of porous materials have been discovered in the last twenty years. *

Corresponding author: [email protected]

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Examples of these classes of materials include metal-organic frameworks (MOFs), zeolitic imidazole frameworks (ZIFs), porous organic cages (POCs) and polymers with intrinsic microposoristy (PIMs), to name a few2, 3. The discovery of these materials has also prompted a significant shift in how we approach the design and optimization of the processes based on these materials. First of all, the rate of the discovery of new materials within these classes outpaces our ability to test them in real applications. On the other hand, the very nature of these materials (the fact that they share several common classes of topologies, but can be built from an infinite variety of building blocks) implies that (a) there are likely to be thousands of yet-to-besynthesized materials within these classes and (b) we can anticipate their crystalline structure by choosing a specific topology and building blocks. These ideas led to the emergence of the computational screening approaches, where adsorption experiments are replaced by their equivalent in the molecular simulation domain, such as grand canonical Monte Carlo simulations (GCMC). With the advances in computer power, molecular simulations can be applied to libraries of thousands and millions of materials. These materials can be both real and virtual, constructed in a lego-like approach from the building blocks following some topology. The idea of these computational screening strategies is to identify the most promising candidates for a particular application before any costly experimental effort is committed. Examples of these large scale screening studies include the pioneering works of Snurr and co-workers4-6, and Smit and co-workers7, 8. Hitherto, these studies have been based on molecular simulations, which can produce a number of equilibrium properties and metrics, such as adsorption isotherms, Henry’s constants, isosteric heats of adsorption, and binary selectivity. These properties are very important for understanding adsorption mechanisms at the molecular level and for the design of adsorbents with bespoke characteristics. However, it is now becoming apparent that on their own they are not enough to assess performance of the material in a real PSA or VSA process. Recently, this has been eloquently demonstrated by Rajagopalan and co-workers9. What further follows from the work of Rajagopalan et al.9, Khurana and Farooq10 and others is that the ranking of porous materials for carbon capture or any other application should be based on their performance in the actual separation process. Indeed, consider application of VSA process to carbon capture from a typical flue gas stream from a coal power station (15% carbon

2

ACS Paragon Plus Environment

Page 2 of 52

Page 3 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

dioxide, 85% nitrogen). The process must recover 90% of carbon dioxide with 95% purity: this is a typical constraint imposed in the design of carbon capture processes. This requirement is made more difficult by the fact that CO2 is the heavy component and is more strongly adsorbed on all practical adsorbents. This means that the CO2 purity and recovery will require an unconventional cycle and the performance will depend not only on the CO2 adsorption but also on the adsorption of the lighter components. In contrast to this, in light component purification, e.g. H2 purification, the most weakly adsorbed component is obtained pure in the adsorption step11, 12, which is the conventional configuration of VSA and PSA cycles. For a given adsorbent material, the configuration of the cycle can be optimized to minimize energy consumption of the process per unit of CO2 captured and to maximize productivity (that is how much carbon dioxide is produced in the evacuation step per unit of volume of the adsorbent per unit of cycle time). These two characteristics, energy penalty and productivity, are the true metrics of the material performance defining the economical cost of the carbon capture process. These two characteristics are also in competition with each other and the same material can meet the design constraints of 90% recovery and 95% purity in a process configured either for high productivity and high energy penalty or for low productivity and low energy penalty. Different combinations of the energy penalty and productivity values corresponding to different configurations of the cycle form a Pareto front, which delineates the best process performance that can be obtained with a given adsorbent. Optimized configurations of the cycle, location of the Pareto front and its width will depend on the material, and cannot be obtained from molecular simulation or equilibrium properties of the material alone. Assessment of the materials or processes based on the Pareto fronts has been already adopted in the process modelling community. Here we mention only some of the recent studies in this field. Haghpanah and co-workers13 developed a PSA/VSA adsorption process simulator and applied it to a 4-step adsorption process in the context of post-combustion carbon capture. Maximizing productivity and minimizing energy penalty under constraints of specified purity and recovery requires multi-objective optimization. The authors employed the nondominated sorting genetic algorithm (NSGA) to identify appropriate operating conditions. Focusing on zeolite 13X as a case study, they explored performance of the process, in terms of the Pareto fronts, as a function of process parameters, such as for example evacuation pressure. In a later

3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

study, Haghpanah and co-workers14 extended this approach to compare performance of various process configurations within 4-, 5- and 6-step cycles. In particular, they observed that under 90% purity and 90% recovery constraint, a 4-step VSA cycle with light product pressurization (LPP) gives the minimum energy penalty of 131 kWh/tonne (14.95 kJ/mol) of CO2 captured at a productivity of 0.57 mol CO2/m3 adsorbent/s. Under 95% purity and 90% recovery constraint, the minimum energy penalty rises to 154 kWh/tonne (17.57 kJ/mol). This provides a useful benchmark for comparison with other studies and materials. In Rajagopalan et al.9, the authors used the same 4-step VSA-LPP process to compare performance of zeolite 13X, two metalorganic frameworks (i.e. Mg-MOF-74, UTSA-16), and activated carbon material CS-AC. Process performance from the analysis of Pareto fronts was used to demonstrate the need for a complete process optimization and analysis, as opposed to selective metrics from molecular simulation or adsorption experiments. In the process modelling studies reviewed above equilibrium adsorption data was provided from experiments. What if, instead, one could use molecular simulations to generate the required data? If possible, this would allow us to extend computational screening methods. In fact several recent studies emerged that build on this idea. For example, Hasan et al.15, used a combination of computational structure characterization techniques to preselect the candidates from a database of zeolite materials, followed by grand canonical Monte Carlo simulations (to generate equilibrium data) and 4-step PSA/VSA simulations (to evaluate process performance) on selected materials. In a similar fashion, Nalaparaju and co-workers16, used a combination of molecular simulations, and process optimization, to explore performance of several MOFs in carbon capture from flue gas, and compared it with the reference results for zeolite 13X. In their work again, 4-step VSA-LPP process was considered, similar to the previous studies mentioned earlier. More recently, Khurana and Farooq10 used the same 4-step VSA-LPP process as a model to perform a large scale screening on 74 real and hypothetical materials, with the equilibrium adsorption data provided by either experiments (where available) or simulations.

4

ACS Paragon Plus Environment

Page 4 of 52

Page 5 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(a)

(b)

Figure 1. Schematic illustration of the concepts associated with multiscale modelling of VSA/PSA systems (see text for explanations and narrative). Figure 1 (a) schematically depicts a typical process involved in the construction of the multiscale model in these studies. In order to understand the remit of the current work, let us take the reader through the required stages. On the left, the process starts with a crystal structure of a porous material (either real or virtual). In principle, in the next stage structural characterization methods can be used to provide information on its porosity, surface area, pore limiting diameter and other characteristics17,

18

. At the next stage, molecular simulations are employed to generate

equilibrium data, such as adsorption isotherms. The equilibrium data is then passed to the process modelling and optimization stage. This is a typical template of the studies reviewed above. However, at the interface between molecular and process modelling several important issues have to be addressed: 1) Isotherm fitting: Modelling of a cyclic VSA or PSA process is based on a set of mass and energy balance equations. In this formulation, the required information on the adsorption

5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

equilibria comes in the form of an appropriate analytical adsorption isotherm, with the parameters of the isotherm obtained from the fitting to the available data. Typically, the Dual Site Langmuir (DSL) equation is assumed (both in the studies where the equilibrium data comes from experiments and from simulations)13, 16, 19-21. Several studies demonstrated that the model is adequate for prediction of adsorption from binary mixtures22-24. The DSL model based on thermodynamically consistent parameters should agree reasonably well with ideal adsorbed solution theory (IAST)22-24 provided that (1) the porous structures considered have micropore volumes which are equally accessible to all adsorbates (e.g. carbon dioxide and nitrogen), and (2) adsorbate molecules are not very dissimilar in terms of interaction energy. Applicability of the DSL model for adsorption in zeolite 13X has been demonstrated by Hu et al25. What is less evident is how to extract parameters of the DSL model from the available experimental or simulation data in a systematic and consistent way. At first glance, it may seem as a trivial issue and not surprisingly, the details of the fitting procedure are rarely reported in the literature. However, as we demonstrate in this work, accuracy of the DSL model can be sensitive to the details of the fitting procedure employed. Although from the error of the fit and from the visual inspection (both on linear and log scales) the DSL isotherms may appear as matching the reference data well (in the case of experimental data this is almost always single component adsorption data), the inaccuracies of the model will manifest themselves when trying to predict adsorption at other conditions, such as temperature, or adsorption from the binary CO2/N2 mixture. Nitrogen isotherm is particularly sensitive to how the DSL model is constructed. It has been demonstrated that the performance of the material in the PSA/VSA process is sensitive to nitrogen adsorption20, 26, 27. In heavy component purification, the purity of the heavy component depends on the amount of the light components (here N2) left in the column before the production step. Therefore, errors in predicting nitrogen adsorption from the binary mixtures will have an impact on the accuracy of the performance prediction in the PSA/VSA process. 2) Structural properties of the pellets: The process modelling requires a number of parameters that cannot be provided by the molecular simulations. To illustrate this point, consider the diagram in panel (b) of Figure 1. The adsorption column shown on the left is packed with adsorbent pellets. Since we are interested in comparison of the performance of different materials to each other, we can fix parameters of the column such as column material, diameter and length in all process simulations. The size of the pellets is another important parameter, influencing the

6

ACS Paragon Plus Environment

Page 6 of 52

Page 7 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

pressure drop across the column and mass transfer process into the pellet. As schematically shown in this figure, each pellet (treated in the figure as a spherical particle) can be seen as a composite system of solid material and macropores. The volume fraction of macropores in the pellet defines pellet porosity. The solid in the pellet, in its turn is a system of microcrystallites of the adsorbent material glued together by a binder. It is important to recognize that molecular simulations can produce properties of the ideal crystal of the adsorbent material (e.g. crystal density, equilibrium adsorption, diffusion coefficients), and therefore in the schematic provided above these properties pertain the microcrystallite regions only. Properties of the pellet (pellet size, pellet density, fraction of binder) have to be either measured experimentally or guessed. In the multiscale simulation studies reviewed above, for the materials for which this data was not available, the authors assumed pellet properties of a reference material (i.e. zeolite 13X). Of course, this is a sensible approach in the absence of any other additional information. In the computational material screening, however, it is important to demonstrate to what extent this assumption is correct. This will be possible only by calculating sensitivity of the process predictions to the properties of the pellet. 3) Equilibrium data from molecular simulations: Prediction of performance of the material in the actual process crucially depends on the accuracy of the available equilibrium data. The molecular simulation community has recently invested a significant effort in developing and improving available force fields for new porous materials28-30. However, we would like to highlight three principal issues with these efforts. Firstly, so far it has not been precisely defined what a good agreement between the experimental and simulation data is. From the perspective of the process modelling, a good agreement in the high-pressure region may not necessarily lead to accurate predictions of the material performance in the actual process, operating at lower partial pressures, where it is much more sensitive to the quality of data and agreement in the Henry’s law region. A good example to explain this point comes from carbon capture processes. The high recovery of CO2 requirement means that the loss of CO2 must be minimized during the adsorption step. This requires an accurate representation of the CO2 concentration near the outlet of the column, where the compositions will be rich in N2. Therefore, an accurate estimate of the Henry’s constant of adsorption for CO2 is crucial. In contrast, in traditional PSA/VSA processes used for light product separation where high purities of the light product are needed, the Henry’s

7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

constant of the second weakest adsorbing component is important, since this component will be the main impurity in the product. Furthermore, in PSA/VSA process, the lowest pressure attained during the regeneration step (production step in carbon capture) is governed by the practical limitations of the existing vacuum pump equipment. As a result, very high Henry’s adsorption constants for the heavy adsorbate may lead to lower working capacity of the material and reduce overall process performance11. These practical considerations are often neglected in the computational screening studies, but they are important for the accuracy of the predictions of the process performance. Finally, so far the molecular simulation community has focused predominantly on the force fields for carbon dioxide adsorption. At the same time, as has been already discussed above, it is known in the process modelling community that it is the nitrogen adsorption and the resulting binary selectivity (among other factors) that governs overall process performance. In this article, we use a well-known and important case of CO2/N2 adsorption in zeolite 13X to explore these issues. This system has been extensively studied in both experiments and molecular simulations thus providing a wealth of reference data. Furthermore, zeolite 13X is the current industrial benchmark and provides an important reference for comparison of the performance of other materials. To further relate our studies to the previous publications, a 4-step vacuum swing adsorption cycle with light product pressurization (VSA-LPP) is considered here. The results section of the article is organized around the three issues raised in this introduction. Briefly, to address the first issue we use a set of adsorption isotherms obtained from molecular simulations. The isotherms are generated for both single component adsorption of carbon dioxide and nitrogen and for their binary mixture (15% CO2 + 85% N2) at different temperatures. This allows us to probe the accuracy of DSL parameters obtained from alternative fitting procedures in predicting adsorption of a binary mixture; and then investigate the impact of the different fitting procedures on the predictions of the process simulation and the resulting Pareto fronts. To address the second issue we use the parameters of the DSL model reported previously based on the experimental data for a specific sample of zeolite 13X. In the first group of simulations we obtain a series of Pareto fronts, corresponding to different values of pellet porosity (while all other characteristics of the material, pellets and the column remain the same). As these Pareto

8

ACS Paragon Plus Environment

Page 8 of 52

Page 9 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

fronts show, substantial variability in the process performance is observed. We set to explore the origin of this effect, which seems to stem from a complex interplay between kinetic and equilibrium factors. In the second group of simulations we obtain a series of Pareto fronts as functions of pellet size (while other parameters remain the same). Again, the optimal size seems to be a result of an interplay between two factors: as the particle size increases, the pressure drop in the column decreases (lowering energy penalty), but the mass transfer into the pellet becomes the controlling factor. To probe the third issue, we return to the molecular simulation studies for the source of equilibrium data and explore how the deviations in the prediction of the nitrogen adsorption isotherms (originated either from a poor analytical fit to the adsorption data, or from inaccuracy of the available atomic force fields in molecular simulation) affect predictions of process performance. In the next Methodology section and the supporting information (SI) we provide technical details of the models and simulations performed in this study.

2. METHODOLOGY AND COMPUTATIONAL DETAILS 2.1.

Process simulation and optimisation

In this section, we describe the process cycle and its characteristics, the process simulator used to model it and the methodology for optimisation of the process cycle configuration. 2.1.1. 4-step VSA cycle with light product pressurization (LPP) The traditional Skarstrom cycle was designed for light product purification and cannot provide optimal performance for heavy component purification processes, such as carbon capture from CO2/N2 mixtures. In contrast, the 4-step VSA cycle with light product pressurization (LPP) shown in Figure 2 is tailored for this task. This cycle consists of the following four steps: adsorption at atmospheric pressure with feed, co-current blowdown to an intermediate pressure to remove nitrogen, counter-current evacuation to low pressure to obtain CO2 product and counter-current pressurisation with light product. This cycle has been shown to recover 90% of CO2 with purity greater than 95% with both process simulations and pilot scale experiments31. The high performance of this cycle can be attributed to the counter current pressurisation with light product, which concentrates the CO2 front closer to the column inlet, thereby minimizing the CO2 losses in the adsorption and blowdown steps. We note, however, that although the 4-step VSA LPP cycle can achieve the desired targets for CO2 purity and recovery it does it only by employing very low evacuation pressures (the lowest limit in this study is 0.01 bar). This does 9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 52

not correspond to the industrial practices where typically pumps would not be able to provide pressures below 0.13 bar32. Consequently, the performance of this cycle strongly depends on the vacuum pressure of the evacuation step, which controls the purity and is responsible for the majority of the energy consumption (70-90%). The blower for the adsorption step and the vacuum pump for the blowdown step are responsible for around 10-30% of the total energy penalty. Within this allocation, the blowdown step is responsible for approximately two thirds of the energy consumption (~67%). Although this 4-step VSA LPP cycle is unlikely to be used by industry due to the very low required evacuation pressure, it became a standard benchmark process for material ranking and for comparison of more advanced process configurations10, 14, 16.

N2

Ads (1 bar)

Feed

Evac (0.010.05 bar)

Bd (0.050.5 bar)

LPP (1 bar)

CO2

Figure 2. Schematic depiction of the 4-step VSA process with LPP. From left to right, adsorption, blowdown, evacuation and light product pressurization steps are shown. 2.1.2. Process Simulation of the 4-step VSA Cycle The 4-step VSA cycle with LPP is simulated using the in-house adsorption cycle simulator, CySim33. Within CySim, process flowsheets can be built using units and modules, such as stream, tank splitter, valve, column, each described by a specific set of mass and energy balance equations. According to Figure 3, CySim model of the 4-step VSA cycle with LPP consists of four streams, two splitters, four valves and an adsorption column. During the adsorption step 1, valves V1 and V2 are open, V3 and V4 are closed. Feed stream F is sent to the adsorption column, AC. The product of this step is stored in a stream unit AP (adsorption product). In the blowdown step, V1 and V2 close, V4 open and the product of the step is stored in the stream unit BP (blowdown product). In the evacuation step, V4 closes and V3 opens. The product of this step

10

ACS Paragon Plus Environment

Page 11 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

is stored in the unit EP (evacuation product). Finally, in LPP step, only V2 is open and adsorption product is used to pressurize the column. AP V2 S2

BP V4

AC

EP

S1 V3 V1

F

Figure 3. Schematic diagram for implementation of the 4-step VSA process with LPP within CySim. This detailed adsorption process simulation is different to the simulation approach employed by Farooq and co-workers (see for example, Ref16). Their simulation model considers the adsorption column only and defines different steps through changes in the boundary conditions of the column. For example, at the switch from the adsorption to the blowdown step, the model by Farooq and co-workers sets the flowrate at the bottom of the column to zero and changes the pressure at the top from atmospheric pressure to the blowdown pressure. Both approaches have benefits and drawbacks: simulating only the column simplifies and speeds-up the simulation but requires the explicit specification of boundary condition profiles representative of the real system; simulating the whole flowsheet is computationally more expensive, but is closer to the real system and the column boundary conditions are given implicitly. One consequence of the more realistic approach of CySim is the need for a number of valves in the flowsheet. These valves control the flowrate between different units of the adsorption process and consequently the rate of pressure change based on the size of the valve (e.g. valve Cv value), its current

11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 52

position (e.g. open or closed) and the pressure difference across the valve. From this, it follows that the Cv values of the valves need to be tailored to the feed flowrate and step times and, thus, that the Cv values need to be included in the optimisation process. This increases the number of optimisation parameters in comparison to the approach by Farooq and co-workers. Despite these differences in the column boundary conditions, CySim can be configured to reproduce reasonably closely reference results from the study by Nalaparaju et al.

16

(see Figure S1 of the

SI document). The behaviour of the adsorption column and the auxiliary units, e.g. valves and feed units, are described by mass, momentum and energy balances. These form a system of differential algebraic equations (DAEs) which in the case of CySim are integrated with the SUNDIALS library33, 34. The complete set of mass, momentum and energy equations, boundary conditions and other technical aspects of the model are provided elsewhere35. Here we highlight the most important governing equations of the model and key assumptions employed. The following are the model assumptions: 1. Ideal gas law is valid. 2. Axial dispersed plug flow in the column takes place. 3. The pressure drop across the column is described by the Ergun equation. 4. Non-isothermal model accounting for energy balances inside the column and the heat transfer between the column, column wall and ambient environment is employed. 5. The adsorption equilibrium is described by the dual-site Langmuir isotherm. 6. The mass transfer mechanism is governed by diffusion through macropores with contributions from Knudsen and molecular diffusion mechanisms25, 36. The mass transfer rate equation is described by linear driving force (LDF) approximation. These assumptions are generally consistent with the previously published process models for the 4-step VSA with LPP, including adsorption in zeolite 13X13, 16, 37. As discussed in detail by Hu et al.25, this system is macropore diffusion controlled. This will be the most likely case for all microporous materials that have pore openings greater than 0.4 nm11, 38. The mass balance in the adsorption column is given by

12

ACS Paragon Plus Environment

Page 13 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

 1 −    ∙  + ∙ + + =0    

(1)

 =   + 1 −  

(2)

   + 1 −   =   − 

   =  ∗ − 

 

(3) (4)

where  ,  and  are the gas phase, macropore gas phase, and total pellet concentrations of component i, respectively.  is the concentration of component i in the crystals of porous

materials within the pellet. The equilibrium adsorbed concentration ∗ is given by an adsorption

isotherm, which is described in Section 3.2. The column void fraction is given by ;  and  are

the pellet porosity and tortuosity; and the macropore and micropore LDF coefficients are  and  , respectively. The diffusive flux is given by  and the interstitial flow velocity by . The

independent parameter  represents the spatial coordinate along the length of the column. The model for the LDF coefficient is based on Glueckauf’s approximation39, which is shown to be

adequate for simulation of equilibrium driven separations:  

15 # =   

where



!,

and

respectively; 

15  =     %

!,

are

!,

15  1 1 =  $  + %&    

molecular

# and 

!,

and

'(

Knudsen

(5) diffusivities

of

component

i,

are macropore and effective macropore diffusivities of

component i, respectively. The micropore LDF coefficient is set to an arbitrary high value so that the sorbate concentration is in equilibrium with the macropore gas phase concentration, which is consistent with the fact that adsorption in zeolite 13X is macropore controlled25.

In our simulations, the system is modelled as non-isothermal with heat transfer allowed between the column and its wall, but pellets and gas phase are kept at the same temperature. The energy balance is formulated as following:

+, + +, ∙  / 1  * * .  . 6 0= + 1 −

++ + 0 + ℎ5 8, − 85       7 2

0 = 95 ̂;,5

85 

+ ℎ5

3(

6 6! 8 − 8<

85 − 8,  + * 7 75 5

13

ACS Paragon Plus Environment

(6)

(7)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 52

+ and . + are the volumetric internal energy and enthalpy, respectively, and . 1 is the partial here, *

molar enthalpy of component i in the fluid phase. The subscripts f, p, w and ∞ indicate the fluid

phase, pellet, wall and surroundings respectively. The temperature, thermal diffusive flux, heat

transfer coefficients and external heat transfer are given by 8, / , ℎ and *. The wall properties are given by the density 95 and the specific heat ̂;,5 .

The momentum balance is described by the Ergun pressure drop equation:

1.751 −  9, >? 150@1 −  − =

+

| | >  4 2

where @ is the fluid viscosity and 9, the fluid density.

(8)

2.1.3. Optimisation of the 4-step VSA Cycle In CySim, we ensure all simulations are checked against specific criteria for achieving a cyclic steady state (CSS) condition. We assume the system has reached CSS once the number of cycles exceeds 300 (this number comes from our experience with similar systems, processes and from a number of preliminary tests) or until the following mathematical criteria are met: FGHI = N

 (/

M M'( ∑2  L3(KL − KL

P

Q

< TUU VWXYZ[X

FJ = \Z]^KLM − KLM'( ^, _ = 1, … , P  < 10 × TUU VWXYZ[X

(9) (10)

where FGHI is root mean-square error, FJ is maximum deviation error and KLM are the fluid and

pellet temperatures, and fluid and solid concentrations along all columns in the system at k-th cycle and for j-th component. The CCS tolerance is set to 10-6 in our simulations.

The optimisation of the 4-step PSA process is carried out by coupling CySim with the Platypus framework in Python40, which uses the third generation of non-dominated sorting genetic algorithm (NSGA-III)41. Briefly, in Platypus we define the optimisation parameters and their range, the number of CySim evaluations as well as the objective functions and constraints. In each iteration, the optimiser picks a set of particular values of the optimisation parameters and these are used to generate an input file which is used to simulate a particular configuration in CySim. The output of CySim is used to check the constraints and to calculate the objective values. These values are passed back to Platypus which then starts the next iteration.

14

ACS Paragon Plus Environment

Page 15 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Optimization using genetic algorithms is a stochastic process, and the exact shape and the points of the Pareto front will somewhat depend on the initial seeding of the decision variables and other factors. The current methodology has been tested in the early stages of the research by exploring the sensitivity of the results to multiple independent runs with different seeding parameters and other conditions. The Pareto front is monitored periodically after 10000 simulations for changes in the minimum energy and maximum productivity and when no significant change in the spread of the Pareto front is observed, it is said to be converged. In line with the requirements for CO2 capture, the optimisation is run with purity and recovery constraints of at least 95% and 90%, respectively. The two objectives are the minimisation of the amount of energy required to recover a mole of CO2 and the maximisation of the productivity of the system. The purity and recovery are calculated from the output of the stream units, which store the M cumulative amount of every component passing through the unit: [L, is the cumulative amount

of component i which has passed through unit j until the end of the k-th cycle. From these values we can calculate purity as the ratio between the amount of CO2 and the total amount of all components in the evacuation product, EP: M ?bYcKde = f

M M'( [g;,de − [g;,de f f

M M'( ∑2 3([g;, − [g;, 

(11)

where N is the number of components. Similarly, the recovery is calculated as the ratio between the CO2 recovered during the evacuation step and the CO2 in the feed stream, F: M XV XYKde = f

M M'( [g;,de − [g;,de f f M M'( [h,de − [h,de f f

(12)

The molar work for the recovery of CO2 is the ratio between the required work per cycle which is given by the pumping work at the different stream units, and the amount of recovered CO2 M iVWZY jVYde = f

M M'( M M'( khM − khM'( + kl; − kl; + kg; − kg; M M'( [g;,de − [g;,de f f

(13)

here, F, BP and EP indicate stream units associated with feed, blowdown product and evacuation product as before, and the pumping work is calculated in CySim through the differential equation

15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

>k = m >

Page 16 of 52

(14)

k0 = 0

(15)

where p is the pumping power given by

m=

q o o

z'( z

r ?vw 8s tu y r−1 ?x

− 1{

z~ '( ?x z~

p rvw 8 s tu y orvw − 1 vw ?vw o n 0.0

− 1{

, c| ?vw > ?x Z[> Vb|WVj , c| ?vw < ?x Z[> c[|WVj

(16)

, VℎXYjc€X

where  is the gas constant and r is the ratio of the molar heat capacities at constant pressure and

constant volume. The subscripts [b and Z\ indicate the neighbouring unit and atmospheric

conditions, respectively. As can be noted from equation (16), we assume 100% efficiency in the energy consumption in the isentropic compression processes. Previous studies by Farooq and coworkers would typically assume 72% efficiency16, whereas experimental studies of Krishnamurthy et al.31 indicate about only 30% efficiency in pilot plant experiments. To compare our results to those by Farooq and co-workers later in the results section we rescaled the results

by Farooq and co-workers to 100% efficiency. Given that the representation of the Pareto front according to a specific realistic efficiency simply requires a single factor, which will move the front up along the y axis, we believe it is more useful for the reader and further use of data to present the unscaled data at 100% efficiency which can be converted into the other cases depending on the assumptions. The productivity is calculated from the recovered CO2, the amount of adsorbent \ and the total cycle time  M ?YV>bc cKde f

M M'( [g;,de − [g;,de f f = \

(17)

For a given recovery and purity and constant adsorbent amount, the productivity depends only on the feed amount and the total cycle time. For example, for the recovery and purity constraints of 90% and 95% at least 90% of the CO2 in the feed is recovered and 95% CO2 concentration in the evacuation step is achieved; thus M ?YV>bc cKde ∝ f

s × ‚ƒ 

(18)

16

ACS Paragon Plus Environment

Page 17 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

The decision variables for the optimization are the durations of the adsorption (tads.), blowdown (tbd.) and evacuation (tevac.) steps, pressures in the evacuation (Pevac.) and the blowdown (Pbd.) steps, feed flow rate (F) and co-efficient of the valves (Cv) used in blowdown and evacuation steps. The typical boundary range for the decision variables are provided in Table 1. For the valves of the blowdown and the evacuation steps, the bounds for the valve coefficients are given by the required minimum and maximum flow rates. For example, the upper limits need to be large enough so that the maximum required pressure change can be achieved in the minimal step time. As described above, the valves control the mass flow rate across unit boundaries and this wide range is necessary so that the operation is not artificially constrained by the valves. For consistency, we have used a column length of 1.0 m, a column radius of 0.1445 m and a pellet diameter of 2.0 mm, as in the studies of Farooq and co-workers16, 37. The feed pressure is fixed at slightly above atmospheric pressure and the adsorption product stream AP is set to 1.0 atm. Valve V2 has a large valve coefficient which ensures that the outlet of the column is at atmospheric pressure. We have also fixed the duration of the LPP step to 20 s, which is the lower bound for the adsorption step time. The optimizations are performed with maximum of 20000 GA evaluations (i.e. CySim simulations). Other parameters relevant to our process simulations are provided in the SI (Table S1). Table 1. Lower and upper bounds for the decision variables used in optimization of the 4-step VSA-LPP process Decision variable Lower bound Upper bound

tAds. (s)

tBd. (s)

tEvac. (s)

PBd. (bar)

PEva. (bar)

Feed flow rate (mol/s)

Cv (Bd. Valve)

Cv (Evac. Valve)

20

1

10

0.05

0.01

0.01

0.05

0.1

200

200

200

0.50

0.03

2.50

0.50

2.0

As mentioned before, we have reproduced the case presented for separation of CO2 and N2 by Farooq and co-workers16, 37 using the same VSA-LPP process cycle. Details of this consistency test is provided in Section 2 of the SI.

2.2.

Molecular Simulations

The grand canonical Monte Carlo (GCMC) molecular simulation technique is employed to predict adsorption of CO2 and N2 in zeolite 13X at four different temperatures. This requires an accurate atomistic model of the zeolite structure and a set of appropriate force field parameters 17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 52

for modelling interatomic interactions. In the next two sections, we explain the procedure adopted to construct the atomistic model of zeolite 13X and provide computational details of the GCMC simulations. 2.2.1. Atomistic Modelling of Zeolite 13X The structure of dehydrated NaX (Na-exchanged faujasite) is modelled based on the crystallographic data reported by Olson42 as summarized in Table 2. Table 2. Crystallographic data for dehydrated NaX Composition

Na88Al88Si104O384

Space group

Fd3

Unit cell parameters

Coordinates of symmetric atoms

A

b

C

α

β

γ

25.099

25.099

25.099

90.0

90.0

90.0

Si1 Al1 O1 O2 O3 O4

Si4+ Al3+ O2O2O2O2-

-0.05381 -0.05524 -0.10990 -0.00110 -0.03460 -0.06930

0.12565 0.03639 0.00030 -0.00280 0.07580 0.07260

0.03508 0.12418 0.10560 0.14160 0.07110 0.18000

In NaX zeolites, which have faujasite (FAU) topology, the Si/Al ratio varies, depending on the sample. It can be determined experimentally, however the distribution of aluminium atoms remains unknown. Using molecular simulations, it has been shown that CO2 adsorption properties in Na-exchanged FAU zeolites are not very sensitive to the distribution of aluminium atoms43. Construction of an atomistic model for Na-FAU is started from a NaX structure with 96 aluminium atoms and 96 silicon atoms (Si/Al ratio of unity), obeying Lowenstein's rule. To construct faujasite-like crystal structures with higher ratios of Si/Al, positions of aluminium atoms should be fractionally occupied by silicon. For Na-FAU, this is achieved by using the procedure proposed by Dubbeldam and co-workers44. In this procedure, firstly the 0th, 5th, 10th, 15th and 20th aluminium atoms in the initial NaX structure file (in which we have 96 aluminium atoms) are substituted by silicon atoms, and then three further aluminium atoms are selected at random within the structure and replaced by silicon atoms. In this way, there are only 88 aluminium atoms left in the unit cell along with 104 silicon atoms, as desired. Following this

18

ACS Paragon Plus Environment

Page 19 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

step, 88 sodium cations are added to the unit cell in random positions. It is known that cations are located close to the aluminium positions45, nevertheless the precise crystallographic location of some sodium cations is still uncertain for NaX and NaY crystals46. Several molecular simulation studies explored mobility and distribution of these cations as well as their impact on adsorption properties45, 47, 48. Here we employ the approach of Smit and co-workers47, where configuration of mobile cations is then equilibrated using the conventional canonical Monte Carlo (MC) simulation. The MC simulations are performed at 323.15 K with 5×105 MC cycles for initialization and 106 MC cycles for production stage where MC moves comprise of traditional translation trials only. The resulting NaX structure with the desired composition of Na88Al88Si104O384 is illustrated in Figure S2. We validated this procedure for the construction of the NaX zeolite models by reproducing adsorption isotherms from reference simulations (Figure S3). 2.2.2. GCMC Simulation of CO2 and N2 Adsorption in Zeolite 13X GCMC molecular simulations are performed at 263.15 K, 298.15 K, 308.15 K and 323.15 K to predict single-component sub-atmospheric adsorption isotherms of CO2 and N2 in zeolite 13X. All GCMC simulations are performed using RASPA 2.0 simulation package49. A maximum of (3×106) MC cycles for CO2, and (5×106) MC cycles for N2 are performed at each pressure point. In these simulations each cycle consists of NMC MC trial moves where NMC is equal to the average number of adsorbate particles in the system. The number of cycles for both equilibration and production stages of the simulations are equal. Three different trial MC moves including insertion, deletion and translation/rotation are employed to obtain equilibrium concentration of CO2 and N2, while Na+ ions remain mobile and their positions are sampled via canonical translation moves only. It has been already demonstrated that the Sodalite (β-cage) cages of NaFAU are not accessible to adsorption from outside the framework47, 48, 50. Therefore, while small Na+ cations are allowed to access these cages, adsorption of larger gas molecules including CO2 and N2 into Sodalite cages is prevented by blocking them in the simulations. In all GCMC simulations reported in this study the zeolite structure is treated as a rigid framework consisting of a 2×2×2 Na-FAU lattice with periodic boundaries in all three dimensions. The 12-6 Lennard-Jones (LJ) potential model is employed to calculate dispersion interactions, while long-range electrostatic interactions are treated using the Ewald summation.

19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Potential cut-off distances for both dispersion and electrostatic interactions (in real space) are set at 24 Å. Both of the above types of interactions are calculated using grid potentials with 0.1 Å resolution. The use of pre-calculated potential grids substantially decreases the computational cost of the simulations. LJ parameters of framework atoms as well as their atomic partial charges are taken from a force field proposed by Vujić and Lyubartsev48. Consistent with the work of Vujić and Lyubartsev, for CO2 molecule, the EPM251, 52 model is used, and N2 is modelled as a 3-site molecule following Murthy et al.53, and Jiang and Sandler54. Details of all force fields used in this study are provided in the SI document. In all GCMC simulations, pressure values are consistently converted to fugacity using the Peng-Robinson equation of state.

3. RESULTS AND DISCUSION 3.1.

Molecular Simulation of Adsorption Process

In Figure 4, we compare the results from molecular simulations to the experimental data reported by Gibson et al.55. We note that the experimental data at 263.15 K is not available in this study. The experimental results are obtained on zeolite 13X pellets containing binder. Although, adsorption data is also available for the crystal samples of 13X, the reference process model of Gibson et al.55 is built on the data collected on the specific pellets, and hence for consistency we use it as a reference. As a result, the simulated adsorption isotherms must be scaled down for direct comparison with the experimental data. The exact amount of binder is not known in this case, however it is typically estimated to be around 15-25%56. Here, we use this property as a scaling factor. Specifically, the degree of scaling (accounting for the presence of the binder) is chosen so that the simulated isotherm agrees with the experimental data at the highest pressure point at 298.15K. Figure 4 shows that simulated adsorption isotherm at 298.15 K matches its experimental counterpart at 1.0 bar once scaled down by 20%, which is well between the amounts expected for mass fraction of binder mentioned earlier. We acknowledge that other protocols can also be used to scale simulated isotherms. Comparison of scaled adsorption isotherms with original GCMC data are shown in Figure S4. As shown in Figure 4 (a) and (b), the resulting simulated isotherms for CO2 are in reasonable agreement with the experimental data especially at higher pressures. In contrast, it is evident from Figure 4 (c) that experimental adsorption isotherms are significantly underpredicted for nitrogen particularly at lower temperatures. In a similar study, Vujić and Lyubartsev48 have also

20

ACS Paragon Plus Environment

Page 20 of 52

Page 21 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

reported considerable underprediction of experimental N2 adsorption in FAU- and MFI-types zeolites using the same set of force fields adopted here. In an attempt to more accurately reproduce experimental nitrogen isotherms, they adjusted parameters of the LJ potential for nitrogen (effectively making the interaction with the zeolite stronger)48. Although this nominally achieves better agreement with the experimental data, their new model lacks the transferability required to predict adsorption of nitrogen across similar classes of porous materials at temperatures below 300 K and/or higher than 1.0 bar48. This observation also emphasizes the need for development of more specialised force fields derived from ab-initio calculations for adsorption of N2 in zeolites and MOFs, which in contrary to CO2, is largely overlooked so far.

(a)

(b)

(c)

21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 52

Figure 4. Adsorption isotherms for CO2 (a, b) and N2 (c) from GCMC simulations in this study (closed symbols). Simulation results are scaled down by 20% to account for mass fraction of the binder. These results are compared to the experimental data from Gibson et al.55 (open symbols).

Here, we also highlight another issue concerning accuracy of molecular simulation predictions. As illustrated in Figure 4 (a) and (b), molecular simulations reproduce experimental CO2 adsorption isotherms more accurately at higher pressures (0.6 – 1.0 bar), whereas predictions between 0.0 to 0.2 bar are less accurate. As already mentioned in the introduction, in the specific process under consideration accurate estimation of the Henry’s constants for CO2 is crucial for the reliability of the process model.

3.2.

Effect of Fitting of the Dual-Site Langmuir Model on Accuracy of Predictions of the Binary Adsorption and Process Performance

We adopt the dual-site Langmuir (DSL) adsorption model to obtain an analytical description of our simulated adsorption isotherms reported in the previous section. For a single component system, the DSL isotherm for species i is defined by:

∗



ƒ = 0 „L, L3(

…L, × ? † 1 + …L, × ?

(19)

ƒ here, L, is saturation capacity of site j with respect to species i, bj,i is affinity of each site

described by the van't Hoff equation: …L, = …!L, exp Š

'‹Œ,Ž G/

. In the van't Hoff equation, .L, is

heat of adsorption at adsorption site j and …!L, is pre-exponential factor. As seen here, there are

ƒ ƒ six parameters ((, , , , …!(, , …!, , ∆.(, and ∆., ) for each gas component i which can be

obtained through fitting of equation (19) to the reference adsorption data using non-linear leastsquare regression. Adsorption of species A from a binary gas mixture of A and B at fixed temperature is described by the extended version of the Dual-Site Langmuir model (extended DSL)57 which is given by: ƒ ‘∗ = ’(,‘

“”,• ×;ז•

ƒ ™+’,‘

(—“”,• ×;ז• —“”,˜ ×;ז˜

“f,• ×;ז•

™

(—“f,• ×;ז• —“f,˜ ×;ז˜

(20)

where yA and yB are mole fractions of components A and B in the gas phase. To obtain physically meaningful parameters for the DLS model, normally the fitting algorithm is guided through a set of mathematical constraints, which also help the algorithm to converge.

22

ACS Paragon Plus Environment

Page 23 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

We note, that the details provided in the literature regarding the set of constraints used in the fitting algorithm are not consistent, and various studies have adopted different approaches for fitting the DSL model to the adsorption data10, 14, 58-61. This is particularly evident in a series of studies conducted on adsorption of CO2 and N2 in zeolite 13X10, 14, 60. In this section, we use molecular simulations to explore three different procedures for fitting DSL parameters. The accuracy of the models is assessed by direct comparison of the DSL predictions for the adsorption from a binary mixture of components at different temperatures to the results obtained from molecular simulations. We demonstrate that unless specific additional constraints are imposed, the DSL predictions may not be accurate. Furthermore, we show that the difference and inaccuracy in DSL models will manifest itself in the process modelling as well. This highlights the importance of adoption of well-documented and validated fitting procedures in order to establish consistency in computational screening and ranking of porous materials. 3.2.1. Three Different Fitting Procedures for the DSL Model The fitting procedure is essentially a least-squares method which minimizes the relative differences between simulated adsorption data and that predicted by the DSL model. Throughout this study, we use the Levenberg-Marquardt algorithm† in Origin Lab‡ for fitting non-linear adsorption isotherm data which includes single-component adsorption isotherms of CO2 and N2 obtained from molecular simulations at 263.15 K, 298.15 K, 308.15 K and 323.15 K (Figure 4). A fitting tolerance of 10-9 is consistently used throughout this study. Fitting procedure 1: The basic procedure (designated in the graphs of the results as Fit 1) consists of two steps: Step 1: In this procedure, the six parameters-form of the DSL model is fitted to the entire set of CO2 adsorption data simultaneously with only two constraints imposed on the fitting algorithm:

…!(,def ≥ …!,def and ∆.( Z[> ∆. ≤ 53 kJ/mol. The former constraint here facilitates easier

convergence of the fitting algorithm (ensuring that there is only one solution), the latter

† ‡

Numerical Recipes in C, Ch. 15.5, Nonlinear Models. Origin (OriginLab, Northampton, MA).

23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 52

constraint sets the upper limit for heat of adsorption in the system as obtained from GCMC simulations. Step 2: Following the previous step, a similar fitting procedure is carried out for nitrogen where the six parameters-form of the DSL model is fitted to the entire set of single-component N2 adsorption data, however this time the following fitting constraints are imposed: ƒ ƒ (,2 = (,de f f

(21)

ƒ ƒ ,2 = ,de f f

(22)

…(,2f ≥ …,2f

(23)

∆.( Z[> ∆. < 25 kJ/mol

(24)

The first two conditions are required for the thermodynamic consistency of the DSL model. In an early study, Myers showed that these conditions are essential for the accuracy of the multicomponent DSL model57. The application (or rather, interpretation) of these consistency criteria has not been very consistent in the recent publications on VSA/PSA process modelling. Here, we use this opportunity to emphasize that correct consistency criteria for the binary system are the ones shown in equations (21) and (22). In a recent study, these criteria have been applied by Gibson et al.55. From the physics of the process point of view, these criteria also make sense: both CO2 and N2 adsorb on both available sites within the structure of zeolite 13X, so that saturation capacities of both sites will be effectively the same. For fitting of the DSL model to nitrogen adsorption data, equation (23) ensures that only one solution is obtained, with the stronger interacting site corresponding to the cations accessible to the adsorbates. This should be the same for both CO2 and N2, given the small but non-negligible quadrupole moment of N2. Equation (24) imposes the upper bound for heat of adsorption, as obtained from GCMC simulations. Fitting procedure 2 (designated in the graphs of the results as Fit 2): The second fitting procedure is similar to procedure 1. The first step is essentially the same. However, in the second step, in addition to equations (21) and (22), which are imposed to fulfil the thermodynamic consistency requirement, the following fitting constraint is also applied:

…(,2f = …,2f

(25)

24

ACS Paragon Plus Environment

Page 25 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

In this procedure, equation (25) emerges as a result of the specific behaviour of nitrogen isotherms in the temperature and pressure range under consideration. As will be shown below, under conditions of interest nitrogen isotherms are essentially straight lines, corresponding to the Henry’s law regime of adsorption. To understand the mathematical need for an additional constraint in this situation, consider DSL model for nitrogen adsorption isotherm in this regime at a specific temperature

ƒ ƒ 2∗ f = (,2 × …(,2f + ,2 × …,2f  × ? f f

(26)

ƒ ƒ Values of (,2 and ,2 are fixed as a result of the thermodynamic consistency constraints f f

(equations (21 and (22)). This still leaves two free parameters, …(,2f and …,2f , to describe a line,

which leads to a mathematically underdetermined problem so that any combination of the two

parameters that gives the same Henry’s law constant can be a viable solution. Several studies10, 58 alleviated this problem by considering equation (25) as an additional constraint. It is tempting to interpret this constraint in terms of homogeneity of nitrogen adsorption compared to carbon dioxide. However, this would not be appropriate. In the range of pressures of interest in postcombustion carbon capture, adsorption isotherms for nitrogen simply do not contain enough information to make any conclusions about site heterogeneity (this is also evident from the smooth variation of the isosteric enthalpy of adsorption for nitrogen below 1.0 bar, as shown in Figure S5 (a) and (b) of the SI. This behaviour is analogous to what is reported for similar porous adsorbents such as Na-, Mg- and Al-rho-ZMOFs by Nalaparaju et al16. Other studies suggested that adsorption of CO2 in many zeolites (including zeolite 13X) is more heterogeneous than nitrogen adsorption under sub-atmospheric pressures62, 63. The complete picture on adsorption heterogeneity of nitrogen can only be exhibited, if one explores N2 adsorption up to very high pressures or lower temperatures where curvature of the adsorption isotherm can be detected. An example of such an isotherm is shown in Figure S6 of the SI. In this procedure, equation (25) is implemented by imposing:

…!(,2f = …!,2f

(27)

∆.( = ∆.

(28)

Results of the DSL model fitted to the single-component CO2 and N2 adsorption isotherms using both fitting procedures 1 and 2 are shown in Figure S7.

25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 52

Fitting procedure 3 (designated in the graphs of the results as Fit 3): Finally, we report on a more systematic fitting procedure, where we recognize the importance of the Henry’s constant of CO2 adsorption in the accuracy of the 4-step VSA LPP model; and that the lower temperature adsorption isotherms contain more information about the parameters of the model than higher temperature data. Hence, the idea of this procedure is to start the fitting from the lowest temperature data available (and obtain a good estimate for as many parameters as possible) and use higher temperature data to extract temperature dependence of the parameters (or in other words, .( and . to ultimately refine the model. This invariably depends on what data and at what temperature is available. Here, the fitting procedure works as follows:

Step 1. We start the fitting procedure from an isotherm obtained at an adequately low temperature where isotherm plateau can be observed at almost 1.0 bar. In our example, the lowest temperature is 263.15 K for which the corresponding adsorption isotherm becomes nearly flat at ~1.0 bar. This indicates that the isotherm is close to the total saturation capacity of the adsorbent at this pressure. In addition to total saturation capacity, the fitted model must also provide an accurate prediction of the Henry’s region of the isotherm, consequently following fitting constraint is imposed:

Œ,def = (,def × …(,def  + ,def × …,def 

(29)

where KH is the Henry’s constant equivalent to the slope of adsorption isotherm when pressure approaches zero. KH can be estimated from the virial adsorption model given by64: ∗ ∗ ∗ ∗ ? = de × X]mžT( + T de + TŸ de

 + T  de

Ÿ … ¡ f f f f

(30)

∗ here, de stands for adsorption at a particular pressure as before and Ci represent virial f

coefficients. In the linear form, this equation can be written as

? ∗ ∗ ∗ ln $ ∗ & = T( + T de + TŸ de

 + T  de

Ÿ … f f f def

(31)

∗ ∗ and at de → 0, de = Œ,def × ?, where Œ,def = exp−T( is the Henry’s constant of f f

adsorption at a particular temperature. Thus, by extrapolating on the virial plot

∗ ¥¦§ f

;

to zero

loading, Œ,def is obtained as the intercept with y-axis. This is shown for adsorption of CO2 at

263.15 K in Figure 5.

26

ACS Paragon Plus Environment

Page 27 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 5. Estimation of the Henry’s constant of adsorption from y-intercept of the virial plot for CO2 adsorption in 263.15 K.

To facilitate easier convergence of the fitting algorithm equation (32) is imposed:

…(,def > …,def

(32)

ƒ ƒ Now, we can obtain total saturation capacity ((,de + ,de ) of the adsorbent by fitting the 4f f

parameter formulation of the DLS model (equation (19)) to adequately low temperature

ƒ ƒ adsorption data. In this case, the fitting parameters include …(,def , …,def , (,de and ,de . We f f

obtain a reasonable initial estimate of the total saturation capacity by looking at the simulated adsorption isotherm close to its highest pressure points (~1.0 atm.), as shown in Figure 6 (a). The more accurate estimate of the total saturation capacity is obtained from the final fitting iteration while both equations (29) and (32) are satisfied.

(a)

(b)

Figure 6. Estimation of the total saturation capacity of CO2 adsorption isotherm at 263.15 K

27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 52

(a); temperature dependence of adsorption affinity of site 1 for CO2 (…(,def ) in the DSL model (b)

Step 2. Following step (1), exactly the same procedure is repeated for adsorption isotherms at 298.15 K, 308.15 K and 323.15 K with the same set of constraints (equations (29) and (32)). The only extra constraint imposed for fitting of the higher temperature isotherms is that saturation capacity of each site in the DSL model must remain equal to the corresponding values of ƒ ƒ (,de and ,de at the lowest temperature isotherm (i.e. 263.15 K), thus f f

ƒ ƒ (,de and ,de = const. f f

(33)

At the end of this step, 2 fitting parameters …(,def and …,def are obtained for each temperature resulting in 10 fitting parameters in total for all four temperatures.

ƒ ƒ Step 3. Once (,de , ,de , …(,def and …,def are calculated for all four temperatures, we can then f f

use the Van't Hoff equation …L, = …!L, exp Š

'‹Œ,Ž G/

 to obtain values of …!(,def , …!,def ,

.(,def and .,def . We also noticed that for calculation of the …!L, , use of the lowest

temperature in the Van't Hoff equation can ultimately lead to a better fitting result. Figure 6 (b) shows how .(,def is calculated from slope of the temperature-dependence plot according to the

ƒ ƒ Van't Hoff equation. At the end of this step, six parameters of (,de , ,de , …!(,def , …!,def , f f

.(,def and .,def are calculated for CO2.

ƒ ƒ Step 4. At this step, we use previously calculated values of (,de , ,de , …!(,def , …!,def , f f

.(,def and .,def from step (3) as our initial values for a new fitting iteration where the DSL

model is fitted to the entire set of adsorption isotherms simultaneously. In this new iterations equations (33), (34) and (35) are imposed as fitting constraints.

…!(,def ≥ …!,def

(34)

.(,def ≥ .,def

(35)

ƒ ƒ The outcome of this iteration is a new set of six parameters ((,de , ,de , …!(,def , …!,def , f f

.(,def , .,def ) constituting the final DSL parameters required for reproducing CO2 adsorption

isotherms. The final fit for CO2 adsorption is shown in Figure 7 (a).

28

ACS Paragon Plus Environment

Page 29 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(a)

(b)

Figure 7. DSL-fitted adsorption isotherms (dashed lines) from the fitting procedure 3, compared to the reference GCMC simulation data (symbols) for CO2 (a) and N2 (b). After calculation of the DSL fitting parameters for CO2, a similar stepwise procedure is carried out for N2 as outlined below: Step 1. For nitrogen, we repeat the same four steps as explained for CO2 before (steps 1 to 4) bearing in mind that in all steps equations (36) and (37) must be imposed on the fitting algorithm to fulfil the thermodynamic consistency requirement. Further, equations (32), (34) and (35) are replaced by equations (38), (39) and (40) respectively. ƒ ƒ (,2 =(,de f f

ƒ ƒ ,2 =,de f f

(36)

.(,2f = .,2f

(39)

(37)

…(,2f = …,2f

(38)

…!(,2f = …!,2f

(40)

The final fit for N2 adsorption data is depicted in Figure 7 (b). 3.2.2. Prediction of the Binary Adsorption Data Here, as in the rest of this article, we consider a case study of adsorption from a binary gas mixture of CO2 (15%) and N2 (85%). Table 3 summarizes three sets of DSL parameters obtained for both CO2 and N2 based on the three different fitting procedures discussed in previous section. Each set of DSL parameters reported in this table are separately fed into equation (20) to predict

29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 52

adsorption isotherms of CO2 and N2 at 283.15 K, 303.15 K and 333.15 K from the binary mixture. The resulting analytical isotherms are compared with the reference binary data obtained from GCMC simulations as shown in Figure 8. Table 3. Fitting parameters of the DSL model for three different fitting procedures CO2

N2

Fitting

Fitting

Fitting

Fitting

Fitting

Fitting

Procedure 1

Procedure 2

Procedure 3

Procedure 1

Procedure 2

Procedure 3

q1 (mmol/g)

3.109

3.109

3.156

q1 (mmol/g)

3.109

3.109

3.156

q2 (mmol/g)

2.271

q2 (mmol/g)

2.271

-1

2.271 -7

9.104×10

2.244 -7

9.233×10

-7

bo1 (bar )

5.577×10

bo2 (bar-1)

-1

2.271 -4

2.244 -5

5.850×10-5

2.557×10-5

5.867×10-5

5.850×10-5

bo1 (bar )

9.104×10

5.867×10

bo2 (bar-1)

7.170×10-7

7.170×10-7

4.884×10-7

∆H1 (J/mol)

47716.56

47716.56

47612.86

∆H1 (J/mol)

8917.61

16062.03

16059.73

∆H2 (J/mol)

37751.67

37751.67

38545.63

∆H2 (J/mol)

19355.37

16062.03

16059.73

As illustrated in figure Figure 8 (a), all three different fitting procedures of the DSL model can equally well predict binary adsorption isotherms of CO2 across all temperatures. In contrast, Figure 8 (b), (c) and (d) show that binary N2 isotherms are sensitive to the way we fit the DSL model to single-component adsorption data. It is clear from these figures that fitting procedure 1 largely overestimates adsorption of nitrogen, however fitting procedures 2 and 3 can both very well predict N2 GCMC isotherms. This is also evident from the fitting parameters provided in Table 3. As can be seen in this table, fitting parameters for CO2 are very similar across three different cases. Given small impact of N2 parameters on CO2 binary isotherm, the predicted binary adsorption isotherms for CO2 are almost identical for all fitting procedures. For nitrogen, DSL parameters of procedure 2 and 3 are very similar, resulting in a close agreement between these two procedures in prediction of nitrogen adsorption. In contrast, for the fitting procedure 1,

…!,2f and .,2f are quite different with their counterparts in procedures 2 and 3 leading to an

inaccurate prediction of N2 binary data.

Procedures 2 and 3 still somewhat overpredict adsorption of nitrogen at pressures above 0.6 bar at the lowest temperature. In general, one would expect that this should lead to worse performance of the process based on the DSL model, compared to the performance of the process according to the molecular simulation predictions. We note, however, that this effect

30

ACS Paragon Plus Environment

Page 31 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

should be rather small. This is because the partial pressure range of 0.6 - 1.0 bar for N2 is only reached during the adsorption step. The evacuation step is responsible for the majority of the energy penalty and this is only indirectly affected by the adsorption in 0.6 - 1.0 bar range. Purity and recovery of the product are governed by the accuracy of the binary desorption isotherms at much lower pressure ranges where the fit of the procedures 2 and 3 is very good and superior to procedure 1.

(a)

(b)

(c)

(d)

Figure 8. Adsorption isotherms of CO2 (a) and N2 (b-d) from binary mixture of 15% CO2, 85% N2 as predicted from DSL model based on three different fitting procedures (lines and star symbols) compared to the reference GCMC simulation data (open symbols).

31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 52

To further analyse accuracy of the three different fitting procedures, we have calculated CO2/N2 selectivities for the isotherms predicted by DSL and compared the results with those obtained from molecular simulation. Here, adsorption selectivity is calculated by

Udef /2f

∗ de $ ∗ f& 2f ¨K = de u K fy 2f

(41)

Figure 9 demonstrates the results of this comparison.

(a)

(b)

(c)

Figure 9. Selectivity Udef /2f predicted from the DSL model as compared with the reference GCMC simulation data. In line with Figure 8, this figure also proves that fitting procedure 1 lacks the accuracy required to correctly predict reference data. Both fitting procedures 2 and 3 adequately reproduce selectivity values obtained from molecular simulation in the whole pressure range, despite the inaccuracies in prediction of adsorption of nitrogen from binary mixtures at pressures above 0.6 bar. 3.2.3. Impact of the Fitting Procedure on Performance of the CO2 Separation Process Here, we examine the effect of different fitting procedures of the DSL model on the ultimate performance of the CO2/N2 separation process using the 4-step VSA-LPP cycle explained earlier in this paper. We feed the DSL parameters reported in Table 3 to our process simulator while keeping all other process parameters the same. The results are shown in the form of three different Pareto fronts in Figure 10.

32

ACS Paragon Plus Environment

Page 33 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 10. Comparisons of the Pareto fronts for CO2/N2 separation in the 4-step VSA-LPP process obtained based on three different fitting approaches of the DSL model. As illustrated in this figure, different procedures adopted for fitting of the DSL model can ultimately lead to different process outcomes. The two ends of the Pareto front obtained from fitting procedure 1 are 30% to 43% different in productivity compared to the two other Pareto fronts, which is not negligible from a process performance point of view. Nevertheless, the Pareto fronts obtained from fitting procedures 2 and 3 show very similar performances as expected from their predictions of binary adsorption data shown in Figure 8 and Figure 9. Figure 10 clearly demonstrates that variation of binary nitrogen isotherm can significantly impact separation performance at the process level. This in turn reiterates importance of accurate prediction of weaker adsorbing component (i.e. nitrogen) both by using rigorous analytical fitting procedures and application of accurate atomic force fields in molecular simulation in order to correctly predict binary adsorption data. As a result of this discussion, we set to adopt the fitting procedure 3 as our standard method for obtaining the DSL parameter in this study, and use the corresponding DSL parameters from Table 3 in the remaining parts of this paper wherever process simulations are performed based on GCMC simulated adsorption isotherms.

3.3.

Effects of the Pellet Morphology on the Process Performance

In this section, we investigate how morphology of the pellets affects process performance. In particular, effects of pellet porosity and pellet size are evaluated. Before presenting the results

33

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

associated with the analysis of the pellet morphology, it is important to establish a reference case study. For this, we focus on the original DSL parameters already reported by Gibson et al.55 based on their experimental adsorption isotherms which is also provided in the SI (Table S5). As discussed previously, this model assumes …(,2f = …,2f . For the same material (Zeolite 13X from

UOP) Hu et al.25

used a combination of mercury porosimetry, volumetric adsorption

measurements, zero length column experiments and process modelling, using CySim to establish structural, equilibrium and kinetic properties of the pellets. Their data was used to develop a model of the pellet which along with other parameters are summarized in Table S1 of the SI. Figure 11 depicts the Pareto front obtained for this material. For comparison, this Figure also provides the result reported by Nalaparaju et al.16 for their Zeochem sample of zeolite 13X. We note, that the study of Nalaparaju et al.16 was based on assuming 72% efficiency of the isentropic compression and pump work. Here we rescale the reference results to 100% efficiency, for consistency, which moves the reference Pareto front to the lower values of energy penalty. Given that we have already established the consistency of the process simulators (see the Methodology section and SI), the differences in the Pareto fronts can be attributed to different kinetic and equilibrium properties of the pellets arising from differences of the pellet porosity, tortuosity and binder mass fraction. Hence, the same material (zeolite 13X) formed into different pellets may lead to different Pareto fronts. Using our data for UOP zeolite 13X as a reference case, we now turn to a more detailed analysis of the effect of structural properties of the pellets on process performance.

Figure 11. Comparison of the 4-step VSA LPP process, using UOP 13X and Zeochem

34

ACS Paragon Plus Environment

Page 34 of 52

Page 35 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

13X pellets. The Zeochem 13X data is from Nalaparaju et al16, rescaled to correspond to 100% efficiency of compression and pump work. 3.3.1. Effect of pellet porosity on the process performance In this section we consider the following scenario. As discussed in the introduction, the pellet is made from crystals, which are glued together by an inert binder, and the structure of a pellet contains the macropores and the microporous crystals. Consider we have a method to prepare pellets with different macroporosities (or in other words, with different pellet void fraction, εp). For simplicity, we consider that the pellet size does not change and proportion of binder to the amount of crystals present in the pellet remains constant. What will be the effect of changing pellet porosity on the process performance? Before we consider the actual results for model systems, let us first try to anticipate possible effects that may come into play: 1) Under all other conditions remaining the same, increasing pellet porosities will lead to lower adsorbent material present in the column. This reduction in the adsorption capacity of the column together with the recovery requirement reduces the upper limit of the amount of gas feed into the system before too much CO2 breaks through. It is expected that the upper limit is approached for the high productivity cases so that as much as possible of the adsorbent material is used. Since the feed amount is proportional to the feed flowrate times the adsorption step time (s × ‚ƒ ),

this product reduces in line with the reduction in adsorption capacity. Due to the dependence of the productivity on the ratio of this product and the total cycle time (equation (18)), the total cycle time has to reduce by the same percentage to keep the productivity constant. For the high productivity case, this means that the optimiser has to reduce the feed amount and the cycle time accordingly. This reduction is expected to increase the specific energy because the length and energy requirement of the blowdown and re-pressurisation steps should stay reasonably constant regardless of the total cycle length and thus the specific energy will increase. In addition, the feed and evacuation steps have a larger decrease and thus the ratio of producing steps to the total cycle length is reducing. Additional effects for fast cycles could also arise from greater pressure drops across the bed which also lead to higher energy penalties. For the low productivity case, the optimiser is not constrained by the upper feed amount limit and thus has more flexibility in adjusting to the lower adsorption capacity. While for the same

35

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 52

productivity a small increase in specific energy is expected, this should not be as large as for the high productivity case. 2) As has been already discussed in the methodology section, in case of the LDF model, the

pellet mass transfer coefficient is given by  = # diffusivity for a macropore is given by 

!,

=



f Gª

­ª ¬

# 



!,

!,

=

(© «ª

f ¬ Gª



!, .

The effective

as presented earlier. Given that

tortuosity can be well approximated by the Bruggeman relation  = 1/  , one can immediately

see that in the macropore controlled adsorption process, changing porosity will have a strong impact on the kinetics of the adsorption process. In particular, increasing porosity should increase the rate of mass transfer into the material. The increased mass transfer will lead to a sharpening of the adsorption front, which, in turn, will increase the allowed feed amount before breakthrough occurs with similar consequences mentioned above. 3) The third effect of changed porosity is due to the changed macropore gas volume. Higher porosity leads to large volume of the gas phase present in the system. This will lead to higher CO2 losses during the blowdown step, which reduces the recovery. In addition, more N2 remains in the column after the blowdown step, which reduces the purity of the recovered CO2. By lowering the blowdown pressure, the purity can be increased at a cost of reduced recovery and

vice versa. Thus, the optimiser needs to adjust the operational conditions to achieve the purity and recovery constraints and it may have complex effect on productivity and specific energy. Having looked at the three main effects, it is clear that points (1) and (2) have the opposite influence. This suggests that an optimal pellet porosity exists which offers the best performance for a given material and separation. With this picture in mind, in Table 4 we present several systems corresponding to different degrees of porosities compared to the reference system. Table 4. Summary of porosity values studied and their corresponding macropore diffusivities for both CO2 and N2 Porosity ( ε P)

Tortuosity (τ)

0.2

5.0

Macropore diffusivity (®¯°±²³ ) [m2/s]

1.24×10-5

36

ACS Paragon Plus Environment

Effective macropore diffusivity (®´¯°±²³ =

µ¶

·

®¯°±²³ ) [m2/s]

0.70×10-6

Page 37 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

*

0.27*

2.6

0.35

2.0

0.5

2.0

0.75

2.0

1.24×10-5

1.30×10-6

1.24×10-5

3.10×10-6

1.24×10-5 1.24×10-5

2.20×10-6 4.80×10-6

Reference system.

In this table for εp = 0.2 the Bruggeman relation was used to set the tortuosity, whereas for other values of porosity (0.35 - 0.75), we use the lowest τ value (i.e. 2.0) typically reported for porous materials38. The Pareto fronts for the materials in Table 4 are shown in Figure 12.

Figure 12. Pareto fronts as function of pellet porosity. The optimisation parameters corresponding to minimum specific energy and maximum productivity conditions of this figure are provided in Tables S6 and S8 of the SI. Additionally, productivity and specific energy values for the cases corresponding to the minimum specific energy and maximum productivity conditions are provided in Tables S7 and S9 of the SI, respectively. Figure 12 clearly shows the competing effects of reduced adsorption capacity and increased mass transfer for a reduction in pellet porosity, with an optimal pellet porosity in the range 0.35 - 0.5. For smaller pellet porosity the system becomes more and more mass-transfer limited which results in a reduction in performance. The reduced mass transfer kinetics leads to a

37

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

broader adsorption front which can be clearly seen in the CO2 gas phase concentration shown in Figure S8. On the other hand, for the larger pellet porosity the effect of reduced adsorption capacity starts to outweigh the increase in mass transfer kinetics. For example, at εp = 0.75 the Pareto front shifts significantly to the left. Interestingly, the performance of the system at εp = 0.75 is almost the same as at εp = 0.27, although the latter column contains significantly more adsorbent. This may have implications on how the pellets of porous materials are prepared with a view to optimize the amount of adsorbent needed for a particular process. Evidently, Figure 12 represents the competition between the enhanced transport in pellets with higher porosity and general deterioration of the process as the amount of active adsorbent in the system decreases. To decouple these two effects, it is instructive to consider a scenario where only one of these factors takes place. Specifically, here we consider a case where the effective macropore diffusivity for all materials is fixed to the highest value, 4.73×10-6 m2s-1. The Pareto fronts for this case are shown in Figure 13.

-6 2 # Figure 13. Pareto fronts for different porosity values at a fixed  ! of 4.73×10 m /s

This figure shows that for systems less constrained by mass transfer kinetics the process performance increases with decreasing pellet porosity. As discussed above, the lower pellet porosity increases the adsorption capacity, which in turn reduces the fraction of the nonproducing cycle steps, i.e. blowdown and re-pressurisation, in the total cycle time. This will increase the productivity and reduce the specific energy because the ratio of producing to non-

38

ACS Paragon Plus Environment

Page 38 of 52

Page 39 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

producing steps is increased. Figure S9 illustrates CO2 gas phase concentration for the case with # fixed  ! . As depicted here, concentration profiles of carbon dioxide at different pellet

porosities are now parallel to each other as the role of one of the competing factors, the mass transfer limitations, is reduced. This is in contrast to the previous case study shown in Figure S8. The optimisation parameters corresponding to minimum specific energy and maximum productivity conditions of Figure 13 are provided in Tables S10 and S12 of the SI. Additionally, productivity and specific energy values for the cases corresponding to the minimum specific energy and maximum productivity conditions are provided in Tables S11 and S13 of the SI, respectively. 3.3.2 Effect of Pellet Size on the Process Performance In this section we consider the effect of the pellet size on the process performance. Before we consider the actual numerical results and analysis, it would be instructive to provide some context on factors influencing the design of PSA/VSA systems and discuss what behaviour we should anticipate given the structure of the equations describing the system which are provided in section 2.1.2 . The size of the adsorption column in PSA/VSA process is inversely proportional to the working capacity of the material and directly proportional to the cycle time11. As a result, a large working capacity and quick pressure swing is desirable to achieve higher productivity and smaller footprint of the process (the regeneration step should be achieved in a few minutes or even seconds)11. At very fast cycles, mass and heat transfer limitations, as well as the gas feed rate and pressure drops across the column become significant11, 65. As explained by Abanades et al11, one would require low pressure drops to achieve quick pressurization, pressure equilibration and evacuation steps. In the most general case, particle size and shape define the bed void fraction, ε. Pressure drop across the bed is governed by the Ergun equation (within the model adopted here), and smaller particle size and lower porosity will lead to higher pressure drop and as a result additional energy penalty. On the other hand, lower porosity of the bed implies higher adsorbent content, higher adsorption capacity, and lower loss of the stronger adsorbing component during the blowdown step. Another effect of pellet size on the performance of the process is associated with the

39

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

kinetics of the adsorption process. As already discussed, the mass transfer of CO2 in zeolite 13X pellets is controlled by diffusion in macropores. According to equation (5), the pellet LDF coefficient is inversely proportional to the square of pellet radius. This means that smaller particles will have higher LDF coefficient and therefore lower resistance to mass transfer. This reduced mass transfer resistance moves the process closer to an equilibrium-controlled process, which should enhance efficiency of the process. As such, there will be a trade-off between the two competing effects of mass transfer resistance and pressure drop corresponding to an optimal particle size. At the limit of large pressure drops or at very fast regeneration cycles, it is advisable to use “structured adsorbents” instead of “traditional packed beds” to achieve efficient systems11, 66. To decouple the effects described above here we consider a more constrained scenario, where the size of the particles is changed, however the mass of the adsorbent and the porosity remain constant. One may see it as a process of agglomerating the smaller particles into pellets of larger size or splitting the current pellets into smaller particles: in these scenarios the total volume of the pellets, the pellet porosity and the adsorbent mass remain the same. At this stage, we are not concerned with the geometrical arrangement of the pellets within the column. Having the mass of the adsorbent and the porosity of the bed fixed, allows us to focus on the competition of just two effects: higher pressure drop across the bed containing smaller particles (as this will still be an effect even if the two systems correspond to the same bed porosity, according to equation (8) versus enhanced mass transfer into pellets of smaller size under all other conditions remaining constant. Figure 14 shows a series of Pareto fronts as a function of pellet size.

40

ACS Paragon Plus Environment

Page 40 of 52

Page 41 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 14. Pareto fronts for pellets of different sizes. In this figure, the reference case is shown by black circles. The results presented in Figure 14 can be summarized as follows: firstly, even under a constrained scenario the pellet size has a strong impact on the process performance. Starting from 1.0 mm pellet size, as we increase the size, the Pareto front widens and moves to higher productivity and lower energy penalties, indicating the overall improvement in the process efficiency. Interestingly, further increase in the pellet size from 1.5 mm to 2.0 mm does not change the position of the Pareto front but significantly expands it, particularly in the region of lower productivity and specific energy values. Furthermore, it seems at 2.0 mm we have reached the optimal pellet size, as further increase in size leads to diminished process performance as indicated by the Pareto front shifting significantly to the left for 3.0 mm pellets. This must be associated with the transport effects across particles of larger size becoming a dominating factor. This can be illustrated in a number of ways and here we focus on just one particular aspect, namely the mole fraction profile of carbon dioxide in the gas phase at the end of the adsorption step. Figure S10 shows these profiles for the points on the Pareto fronts corresponding to the minimum specific energy (a) and maximum productivity (b) for each pellet size, while Tables S14 - S17 of the SI provide more detailed information on the parameters of the cycles and the performance metrics corresponding to these points. From this Figure, it is clear that as we increase the pellet size, carbon dioxide composition profile develops from a relatively sharp front for 1.0 mm pellets to progressively more diffuse profiles,

41

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

indicating the increasing role of the mass transfer effects. This is particularly evident for 4.0 mm pellets under all conditions.

3.4.

Impact of the Accuracy of Nitrogen Adsorption Isotherms on the Performance of the Separation Process

In this section, we return to the molecular simulations studies comparing them again to the predictions based on the experimental adsorption isotherms. This time we focus on the effect of accuracy of nitrogen adsorption isotherms on overall performance of the separation process. In section 3.2.3, we showed that prediction of performance of the material in the process level considerably depends on accurate prediction of adsorption isotherms. We demonstrated how inaccurate analytical fitting of true adsorption isotherms for N2 can adversely change the picture we obtain for performance of the process. Likewise, inaccuracy in predictions from molecular simulations, before even any analytical fitting is performed, can significantly affect observations at the process modelling stage. As shown in Figure 4, molecular simulation is unable to provide a reasonable prediction for adsorption of nitrogen in zeolite 13X. This observation however is not limited to the particular nitrogen model used in this study or to adsorption in zeolite 13X: the problem is more universal covering almost all types of generic nitrogen force fields including two of the most popular N2 models, TraPPE67 and UFF68. Several authors have reported inadequacy of these models for prediction of nitrogen adsorption in many zeolitic and metalorganic frameworks. Examples include adsorption of N2 in FAU and MFI type zeolites with different Si/Al ratio48, adsorption of N2 in open metal-site Mg-MOF-7469, and nitrogen adsorption in TIF-A1, MOF-177, CALF-15, UMCM-1, ZIF-6870 just to name a few. This problem arises from the fact that the molecular simulation community has predominantly focused on developing specialised force fields for adsorption of CO2 in new porous materials. Accurate force fields derived from ab-initio methods for nitrogen adsorption in porous materials is extremely scarce. Here, we illustrate the effect of inadequate nitrogen force fields on the overall performance of the CO2 capture process. For this we consider performance of the 4-step VSA-LPP process based on two different sets of adsorption isotherms (test cases): (1) simulated CO2 and N2 adsorption isotherms as obtained from GCMC molecular simulations (shown in Figure 4), (2) CO2 adsorption isotherms obtained from GCMC simulation, and experimental N2 isotherms reported by Gibson et al.55 (also shown in Figure 4). Here, we note that experimental

42

ACS Paragon Plus Environment

Page 42 of 52

Page 43 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

isotherms are not reported for 263.15 K by Gibson et al.55, therefore to consistently use adsorption data across the same range of temperatures from 263.15 K to 323.15 K, we have reconstructed the experimental isotherm at 263.15 K using the DSL parameters reported by Gibson et al.55. For test case 1, the DSL parameters obtained from the fitting procedure 3 are already provided in Table 3. For test case 2 however, a new set of DSL parameters had to be obtained for the experimental N2 isotherms. These parameters are provided in Table S18 of the SI. Overall performance of the separation process for both test cases are depicted in Figure 15.

Figure 15. Overall performance of the separation process based on experimental and simulated nitrogen adsorption isotherms Figure 15 shows that process performance predictions where molecular simulation isotherms for nitrogen were replaced with the model based on the experimental data agree much better with the reference Pareto front. Although, the carbon dioxide adsorption isotherms from molecular simulations are also not in a perfect agreement with the experimental data, the current observation unambiguously highlights the importance of the nitrogen adsorption data.

4. CONCLUSION The original motivation of this article stems from an elegant idea: what, if given just a crystal structure of either real or hypothetical material, one could predict the behaviour of this material in an actual adsorption process? If it was possible, it would significantly advance the current state-of-the-art in the design and optimization of adsorption separation processes for various applications, such as carbon capture. Not surprisingly, several recent articles that started to explore this idea attracted substantial interest from both molecular and process modelling

43

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

communities. Our attempt to implement this idea encountered several unaddressed challenges at the interface between the molecular and process levels of description, which must be resolved in order to make multiscale computational screening strategies possible. The main outcomes of this article can be summarized as follows: 1) Accuracy of the DSL model does depend on the details of the protocol for fitting single component isotherms. Visual inspection of the CO2 adsorption isotherms alone, even on the log plot in the Henry’s law regime is not a sufficient test for the ability of the model to provide accurate predictions of the adsorption from binary mixtures. For the computational screening, based on equilibrium data from molecular simulations, it is vital to check the accuracy of the DSL model (or any other model used, such as IAST, etc.) against direct simulation of adsorption from binary mixtures at selected conditions. Protocols for obtaining accurate parameters of the DSL models are proposed in this study: they utilize constraints associated with the thermodynamic consistency criteria, additional constraints for fitting the parameters of the model for nitrogen, and systematically build DSL parameters starting from the lower temperature adsorption data first, through several stages of refinement. Although this might seem trivial, a simple literature review quickly reveals that the current practice has overlooked significance of employing a consistent and proven fitting protocol at least for the DSL model. The process and molecular simulation community are encouraged to adopt the protocols presented in this article or at least report the details of the fitting procedure to ensure consistency across the results. Is it possible to bypass the issues associated with the DSL model and feed the simulation data directly into the process simulations? Process simulations require an analytical function for a binary adsorption isotherm so that adsorption can be predicted at an arbitrary value of pressure and temperature during the cycle. In principle, one could obtain a large number of GCMC simulation data points at various conditions of composition, temperature and pressure and fit the data to some analytical function, avoiding the DSL altogether. This, however, may prove to be computationally very costly to achieve the same accuracy of prediction as the DSL model – this is something to be investigated in further studies. 2) Performance of a material in a PSA/VSA process does depend on the porosity of the pellet. A value of the pellet porosity exist at which performance of this material is optimal, representing a compromise between mass-transfer limitations in the macropore space of the pellet and lower

44

ACS Paragon Plus Environment

Page 44 of 52

Page 45 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

efficiency of the process as a result of lower amount of adsorbent material present in the column (the sources of this lower efficiency are discussed in the Results section). Discussion of how to control pellet porosity (and to what extent it is possible) is outside of the scope of this article. Here we would like to emphasize that it is important to at least adequately understand the effect of pellet porosity on the process performance in order for computational screening of porous materials to become a reliable and robust strategy. 3) Performance of a material in a PSA/VSA process does depend on the pellet size. A value of the pellet size exists at which performance of the material is optimal, representing a compromise between mass-transfer limitations in the macropore space of the pellet and energy cost associated with the pressure drop across the bed. 4) The above observations, pertaining the morphology of the pellet, have been obtained in this article for a specific material (zeolite 13X) and for a specific pellet (UOP). This naturally poses a question whether the optimal pellet porosity and pellet size are specific for each material or represent some general figures of merit characteristic for the process. In the former case, this implies that both pellet porosity and pellet size should be included as decision variables in the optimization protocols in the material screening procedures. This will be investigated in further studies. Conclusions 2 and 3 also suggest that it might be worthwhile for material science community to further invest some effort in development of new techniques for production of adsorbents with tunable pellet size and porosity, which in turn will help to improve performance of the separation process. 5) Over the last ten years, molecular simulation and computational chemistry community has dedicated substantial amount of effort to developing accurate force fields for carbon dioxide adsorption in novel porous materials, such as MOFs and ZIFs. Yet, separation process by definition, involves at least two components. In the process modelling community, it has been understood for some time that adsorption of the weaker adsorbing component and the resulting selectivity has a strong impact on the performance as a whole and this property of the adsorbent material must be captured accurately. A single case of zeolite 13X studied here, clearly demonstrates that even for a very well-known and industrially used material, accurate prediction of the adsorption isotherms (especially for nitrogen) is still challenging. It will be therefore imperative to invest special effort in developing accurate and transferable molecular models and

45

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

force field parameters (based on comprehensive DFT calculations) for adsorption of nitrogen and other relevant light components in MOFs, ZIFs, and zeolites. 6) There are several materials databases available which collect experimental adsorption data. However, for the experimental equilibrium adsorption data to be usable in process simulation, a complete set of experimental measurements needs to be carried out consistently on the same sample of a material. Otherwise, the inconsistencies arising from measurements conducted on various samples of the same materials will lead to confusion and unrealistic prediction of process performance. In fact, lack of completeness of the reported experimental data in the literature is a major obstacle for performing reliable process simulations. As shown in this study, full characterization of the pelletized adsorbent, including pellet size, porosity, equilibrium adsorption data for a wide range of temperatures and pressures (including low temperature data where saturation capacity of the adsorbent can be measured) are required for these types of process models. We acknowledge here that some of the reported patterns in the behaviour of the processes in response to material properties and adsorption isotherms must have been already observed or recognized in the process modelling community, but not necessarily in the molecular simulation community. In contrast, the process modelling community may not be aware of either capabilities or challenges in molecular simulations of adsorption. One aspiration of this article is to make a step towards consistent communication between both communities and levels of description, without which multiscale process optimization strategies are not be possible.

CONFLICT OF INTEREST There are no conflicts to declare.

ACKNOWLEDGEMENTS This work has been supported by the UK Engineering and Physical Sciences Research Council (EPSRC), grant EP/N007859/1. This work has made use of the computational resources provided by the Edinburgh Compute and Data Facility (ECDF) (http://www.ecdf.ed.ac.uk/). We would also like to acknowledge Dr. Ishan Sharma for his insights and help in testing Genetic Algorithms for VSA process optimization.

AUTHOR CONTRIBUTIONS All authors contributed to discussions and project development. In the context of this article, the concept of combining molecular simulation with process simulation for multiscale screening of 46

ACS Paragon Plus Environment

Page 46 of 52

Page 47 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

porous materials was proposed by L. Sarkisov and S. Brandani. The CySim simulation code was previously developed by D. Friedrich and S. Brandani. The code was modified and re-optimised for the purpose of this study by D. Friedrich. Data analysis and interpretation was carried out by all authors. Molecular simulations performed in this study including (structural modelling of Zeolite 13X, pure and binary grand canonical Monte Carlo simulations, calculations of CO2/N2 selectivities) and analytical fitting of the molecular simulation data using the DSL model were performed by A. H. Farmahini. Process simulations reported in sections 3.2.3 and 3.4 are also performed and analysed by A. H. Farmahini. Process simulations reported in sections 3.3, 3.3.1 and 3.3.2 are performed and analysed by S. Krishnamurthy.

47

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

REFERENCES 1. Boot-Handford, M. E.; Abanades, J. C.; Anthony, E. J.; Blunt, M. J.; Brandani, S.; Mac Dowell, N.; Fernandez, J. R.; Ferrari, M.-C.; Gross, R.; Hallett, J. P.; Haszeldine, R. S.; Heptonstall, P.; Lyngfelt, A.; Makuch, Z.; Mangano, E.; Porter, R. T. J.; Pourkashanian, M.; Rochelle, G. T.; Shah, N.; Yao, J. G.; Fennell, P. S., Carbon capture and storage update. Energy & Environmental Science 2014, 7, (1), 130-189. 2. Furukawa, H.; Cordova, K. E.; O’Keeffe, M.; Yaghi, O. M., The Chemistry and Applications of Metal-Organic Frameworks. Science 2013, 341, (6149). 3. Das, S.; Heasman, P.; Ben, T.; Qiu, S., Porous Organic Materials: Strategic Design and Structure– Function Correlation. Chemical Reviews 2017, 117, (3), 1515-1563. 4. Wilmer, C. E.; Leaf, M.; Lee, C. Y.; Farha, O. K.; Hauser, B. G.; Hupp, J. T.; Snurr, R. Q., Large-scale screening of hypothetical metal-organic frameworks. Nat Chem 2011, 4, (2), 83-89. 5. Yazaydın, A. Ö.; Snurr, R. Q.; Park, T.-H.; Koh, K.; Liu, J.; LeVan, M. D.; Benin, A. I.; Jakubczak, P.; Lanuza, M.; Galloway, D. B.; Low, J. J.; Willis, R. R., Screening of Metal−Organic Frameworks for Carbon Dioxide Capture from Flue Gas Using a Combined Experimental and Modeling Approach. Journal of the American Chemical Society 2009, 131, (51), 18198-18199. 6. Colon, Y. J.; Snurr, R. Q., High-throughput computational screening of metal-organic frameworks. Chemical Society Reviews 2014, 43, (16), 5735-5749. 7. Lin, L.-C.; Berger, A. H.; Martin, R. L.; Kim, J.; Swisher, J. A.; Jariwala, K.; Rycroft, C. H.; Bhown, A. S.; Deem, M. W.; Haranczyk, M.; Smit, B., In silico screening of carbon-capture materials. Nat Mater 2012, 11, (7), 633-641. 8. Huck, J. M.; Lin, L.-C.; Berger, A. H.; Shahrak, M. N.; Martin, R. L.; Bhown, A. S.; Haranczyk, M.; Reuter, K.; Smit, B., Evaluating different classes of porous materials for carbon capture. Energy & Environmental Science 2014, 7, (12), 4132-4146. 9. Rajagopalan, A. K.; Avila, A. M.; Rajendran, A., Do adsorbent screening metrics predict process performance? A process optimisation based study for post-combustion capture of CO2. International Journal of Greenhouse Gas Control 2016, 46, 76-85. 10. Khurana, M.; Farooq, S., Adsorbent Screening for Postcombustion CO2 Capture: A Method Relating Equilibrium Isotherm Characteristics to an Optimum Vacuum Swing Adsorption Process Performance. Industrial & Engineering Chemistry Research 2016, 55, (8), 2447-2460. 11. Abanades, J. C.; Arias, B.; Lyngfelt, A.; Mattisson, T.; Wiley, D. E.; Li, H.; Ho, M. T.; Mangano, E.; Brandani, S., Emerging CO2 capture systems. International Journal of Greenhouse Gas Control 2015, 40, 126-166. 12. Ruthven, D. M.; Farooq, S.; Knaebel, K. S., Pressure Swing Adsorption. John Wiley & Sons: New York, 1993; p 352+XXIV. 13. Haghpanah, R.; Majumder, A.; Nilam, R.; Rajendran, A.; Farooq, S.; Karimi, I. A.; Amanullah, M., Multiobjective Optimization of a Four-Step Adsorption Process for Postcombustion CO2 Capture Via Finite Volume Simulation. Industrial & Engineering Chemistry Research 2013, 52, (11), 4249-4265. 14. Haghpanah, R.; Nilam, R.; Rajendran, A.; Farooq, S.; Karimi, I. A., Cycle synthesis and optimization of a VSA process for postcombustion CO2 capture. AIChE Journal 2013, 59, (12), 4735-4748. 15. Hasan, M. M. F.; First, E. L.; Floudas, C. A., Cost-effective CO2 capture based on in silico screening of zeolites and process optimization. Physical Chemistry Chemical Physics 2013, 15, (40), 17601-17618. 16. Nalaparaju, A.; Khurana, M.; Farooq, S.; Karimi, I. A.; Jiang, J. W., CO2 capture in cationexchanged metal–organic frameworks: Holistic modeling from molecular simulation to process optimization. Chemical Engineering Science 2015, 124, 70-78. 17. Sarkisov, L.; Harrison, A., Computational structure characterisation tools in application to ordered and disordered porous materials. Molecular Simulation 2011, 37, (15), 1248-1257.

48

ACS Paragon Plus Environment

Page 48 of 52

Page 49 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

18. Willems, T. F.; Rycroft, C. H.; Kazi, M.; Meza, J. C.; Haranczyk, M., Algorithms and tools for highthroughput geometry-based analysis of crystalline porous materials. Microporous and Mesoporous Materials 2012, 149, (1), 134-141. 19. Xiao, P.; Zhang, J.; Webley, P.; Li, G.; Singh, R.; Todd, R., Capture of CO2 from flue gas streams with zeolite 13X by vacuum-pressure swing adsorption. Adsorption 2008, 14, (4), 575-582. 20. Maring, B. J.; Webley, P. A., A new simplified pressure/vacuum swing adsorption model for rapid adsorbent screening for CO2 capture applications. International Journal of Greenhouse Gas Control 2013, 15, 16-31. 21. Leperi, K. T.; Snurr, R. Q.; You, F., Optimization of Two-Stage Pressure/Vacuum Swing Adsorption with Variable Dehydration Level for Postcombustion Carbon Capture. Industrial & Engineering Chemistry Research 2016, 55, (12), 3338-3350. 22. Chen, Y. D.; Ritter, J. A.; Yang, R. T., Nonideal adsorption from multicomponent gas mixtures at elevated pressures on a 5A molecular sieve. Chemical Engineering Science 1990, 45, (9), 2877-2894. 23. Pan, Z.; Connell, L. D., Comparison of adsorption models in reservoir simulation of enhanced coalbed methane recovery and CO2 sequestration in coal. International Journal of Greenhouse Gas Control 2009, 3, (1), 77-89. 24. Bartholdy, S.; Bjørner, M. G.; Solbraa, E.; Shapiro, A.; Kontogeorgis, G. M., Capabilities and Limitations of Predictive Engineering Theories for Multicomponent Adsorption. Industrial & Engineering Chemistry Research 2013, 52, (33), 11552-11563. 25. Hu, X.; Mangano, E.; Friedrich, D.; Ahn, H.; Brandani, S., Diffusion mechanism of CO2 in 13X zeolite beads. Adsorption 2014, 20, (1), 121-135. 26. Rajagopalan, A. K. Material selection and process design for adsorptive CO2 capture. University of Alberta, Edmonton, Canada, 2015. 27. Rajagopalan, A. K.; Avila, A. M.; Rajendran, A., Do adsorbent screening metrics predict process performance? A process optimisation based study for post-combustion capture of CO2. International Journal of Greenhouse Gas Control 2016, 46, (Supplement C), 76-85. 28. Fang, H.; Demir, H.; Kamakoti, P.; Sholl, D. S., Recent developments in first-principles force fields for molecules in nanoporous materials. Journal of Materials Chemistry A 2014, 2, (2), 274-291. 29. McDaniel, J. G.; Li, S.; Tylianakis, E.; Snurr, R. Q.; Schmidt, J. R., Evaluation of Force Field Performance for High-Throughput Screening of Gas Uptake in Metal–Organic Frameworks. The Journal of Physical Chemistry C 2015, 119, (6), 3143-3152. 30. Mercado, R.; Vlaisavljevich, B.; Lin, L.-C.; Lee, K.; Lee, Y.; Mason, J. A.; Xiao, D. J.; Gonzalez, M. I.; Kapelewski, M. T.; Neaton, J. B.; Smit, B., Force Field Development from Periodic Density Functional Theory Calculations for Gas Separation Applications Using Metal–Organic Frameworks. The Journal of Physical Chemistry C 2016, 120, (23), 12590-12604. 31. Krishnamurthy, S.; Rao, V. R.; Guntuka, S.; Sharratt, P.; Haghpanah, R.; Rajendran, A.; Amanullah, M.; Karimi, I. A.; Farooq, S., CO2 capture from dry flue gas by vacuum swing adsorption: A pilot plant study. AIChE Journal 2014, 60, (5), 1830-1842. 32. Sircar, S.; Golden, T. C., Purification of Hydrogen by Pressure Swing Adsorption. Separation Science and Technology 2000, 35, (5), 667-687. 33. Friedrich, D.; Ferrari, M.-C.; Brandani, S., Efficient Simulation and Acceleration of Convergence for a Dual Piston Pressure Swing Adsorption System. Industrial & Engineering Chemistry Research 2013, 52, (26), 8897-8905. 34. Hindmarsh, A. C.; Brown, P. N.; Grant, K. E.; Lee, S. L.; Serban, R.; Shumaker, D. E.; Woodward, C. S., SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans. Math. Softw. 2005, 31, (3), 363-396. 35. Beck, J.; Friedrich, D.; Brandani, S.; Fraga, E. S., Multi-objective optimisation using surrogate models for the design of VPSA systems. Computers & Chemical Engineering 2015, 82, 318-329.

49

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

36. Giesy, T. J.; Wang, Y.; LeVan, M. D., Measurement of Mass Transfer Rates in Adsorbents: New Combined-Technique Frequency Response Apparatus and Application to CO2 in 13X Zeolite. Industrial & Engineering Chemistry Research 2012, 51, (35), 11509-11517. 37. Khurana, M.; Farooq, S., Simulation and optimization of a 6-step dual-reflux VSA cycle for postcombustion CO2 capture. Chemical Engineering Science 2016, 152, (Supplement C), 507-515. 38. Ruthven, D. M., Principles of Adsorption and Adsorption Processes. John Wiley & Sons: 1984. 39. Glueckauf, E., Theory of chromatography. Part 10.-Formulae for diffusion into spheres and their application to chromatography. Transactions of the Faraday Society 1955, 51, (0), 1540-1551. 40. Hadka, D. Platypus - Multiobjective Optimization in Python, https://github.com/ProjectPlatypus/Platypus, 2017. 41. Deb, K.; Jain, H., An Evolutionary Many-Objective Optimization Algorithm Using Reference-PointBased Nondominated Sorting Approach, Part I: Solving Problems With Box Constraints. IEEE Transactions on Evolutionary Computation 2014, 18, (4), 577-601. 42. Olson, D. H., The crystal structure of dehydrated NaX. Zeolites 1995, 15, (5), 439-443. 43. Fang, H.; Kamakoti, P.; Ravikovitch, P. I.; Aronson, M.; Paur, C.; Sholl, D. S., First principles derived, transferable force fields for CO2 adsorption in Na-exchanged cationic zeolites. Physical Chemistry Chemical Physics 2013, 15, (31), 12882-12894. 44. Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q. User Manual - RASPA 2.0: Molecular Software Package for Adsorption and Diffusion in (Flexible) Nanoporous Materials, 2.0; 2015. 45. Plant, D. F.; Maurin, G.; Jobic, H.; Llewellyn, P. L., Molecular Dynamics Simulation of the Cation Motion upon Adsorption of CO2 in Faujasite Zeolite Systems. The Journal of Physical Chemistry B 2006, 110, (29), 14372-14378. 46. García-Sánchez, A.; Ania, C. O.; Parra, J. B.; Dubbeldam, D.; Vlugt, T. J. H.; Krishna, R.; Calero, S., Transferable Force Field for Carbon Dioxide Adsorption in Zeolites. The Journal of Physical Chemistry C 2009, 113, (20), 8814-8820. 47. Joos, L.; Swisher, J. A.; Smit, B., Molecular Simulation Study of the Competitive Adsorption of H2O and CO2 in Zeolite 13X. Langmuir 2013, 29, (51), 15936-15942. 48. Bojan, V.; Alexander, P. L., Transferable force-field for modelling of CO 2 , N 2 , O 2 and Ar in all silica and Na + exchanged zeolites. Modelling and Simulation in Materials Science and Engineering 2016, 24, (4), 045002. 49. Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q., RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials. Molecular Simulation 2016, 42, (2), 81-101. 50. Hutson, N. D.; Reisner, B. A.; Yang, R. T.; Toby, B. H., Silver Ion-Exchanged Zeolites Y, X, and LowSilica X:  Observations of Thermally Induced Cation/Cluster Migration and the Resulting Effects on the Equilibrium Adsorption of Nitrogen. Chemistry of Materials 2000, 12, (10), 3020-3031. 51. Harris, J. G.; Yung, K. H., Carbon Dioxide's Liquid-Vapor Coexistence Curve And Critical Properties as Predicted by a Simple Molecular Model. The Journal of Physical Chemistry 1995, 99, (31), 1202112024. 52. Makrodimitris, K.; Papadopoulos, G. K.; Theodorou, D. N., Prediction of Permeation Properties of CO2 and N2 through Silicalite via Molecular Simulations. The Journal of Physical Chemistry B 2001, 105, (4), 777-788. 53. Murthy, C. S.; Singer, K.; Klein, M. L.; McDonald, I. R., Pairwise additive effective potentials for nitrogen. Molecular Physics 1980, 41, (6), 1387-1399. 54. Jiang, J.; Sandler, S. I., Separation of CO2 and N2 by Adsorption in C168 Schwarzite:  A Combination of Quantum Mechanics and Molecular Simulation Study. Journal of the American Chemical Society 2005, 127, (34), 11989-11997. 55. Gibson, J. A. A.; Mangano, E.; Shiko, E.; Greenaway, A. G.; Gromov, A. V.; Lozinska, M. M.; Friedrich, D.; Campbell, E. E. B.; Wright, P. A.; Brandani, S., Adsorption Materials and Processes for

50

ACS Paragon Plus Environment

Page 50 of 52

Page 51 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Carbon Capture from Gas-Fired Power Plants: AMPGas. Industrial & Engineering Chemistry Research 2016, 55, (13), 3840-3851. 56. Brandani, S.; Mangano, E.; Sarkisov, L., Net, excess and absolute adsorption and adsorption of helium. Adsorption 2016, 22, (2), 261-276. 57. Myers, A. L., Activity coefficients of mixtures adsorbed on heterogeneous surfaces. AIChE Journal 1983, 29, (4), 691-693. 58. Mathias, P. M.; Kumar, R.; Moyer, J. D.; Schork, J. M.; Srinivasan, S. R.; Auvil, S. R.; Talu, O., Correlation of Multicomponent Gas Adsorption by the Dual-Site Langmuir Model. Application to Nitrogen/Oxygen Adsorption on 5A-Zeolite. Industrial & Engineering Chemistry Research 1996, 35, (7), 2477-2483. 59. Goel, C.; Bhunia, H.; Bajpai, P. K., Prediction of Binary Gas Adsorption of CO2/N2 and Thermodynamic Studies on Nitrogen Enriched Nanostructured Carbon Adsorbents. Journal of Chemical & Engineering Data 2017, 62, (1), 214-225. 60. Mason, J. A.; Sumida, K.; Herm, Z. R.; Krishna, R.; Long, J. R., Evaluating metal-organic frameworks for post-combustion carbon dioxide capture via temperature swing adsorption. Energy & Environmental Science 2011, 4, (8), 3030-3040. 61. Bae, T.-H.; Hudson, M. R.; Mason, J. A.; Queen, W. L.; Dutton, J. J.; Sumida, K.; Micklash, K. J.; Kaye, S. S.; Brown, C. M.; Long, J. R., Evaluation of cation-exchanged zeolite adsorbents for postcombustion carbon dioxide capture. Energy & Environmental Science 2013, 6, (1), 128-138. 62. ERTAN, A. CO2, N2 and Ar Adsorption on Zeolites İzmir Institute of Technology, İzmir, Turkey 2004. 63. Fischer, M.; Bell, R. G., Influence of Zeolite Topology on CO2/N2 Separation Behavior: ForceField Simulations Using a DFT-Derived Charge Model. The Journal of Physical Chemistry C 2012, 116, (50), 26449-26463. 64. Sing, K. S. W.; Rouquerol, F.; Rouquerol, J., 5 - Classical Interpretation of Physisorption Isotherms at the Gas–Solid Interface. In Adsorption by Powders and Porous Solids (Second Edition), Academic Press: Oxford, 2014; pp 159-189. 65. Webley, P. A., Adsorption technology for CO2 separation and capture: a perspective. Adsorption 2014, 20, (2), 225-231. 66. Ruthven, D. M.; Thaeron, C., Performance of a parallel passage adsorbent contactor. Gas Separation & Purification 1996, 10, (1), 63-73. 67. Potoff, J. J.; Siepmann, J. I., Vapor–liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen. AIChE Journal 2001, 47, (7), 1676-1682. 68. Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M., UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. Journal of the American Chemical Society 1992, 114, (25), 10024-10035. 69. Dzubak, A. L.; Lin, L.-C.; Kim, J.; Swisher, J. A.; Poloni, R.; Maximoff, S. N.; Smit, B.; Gagliardi, L., Ab initio carbon capture in open-site metal–organic frameworks. Nature Chemistry 2012, 4, 810. 70. Provost, B. An Improved N2 Model for Predicting Gas Adsorption in MOFs and Using Molecular Simulation to Aid in the Interpretation of SSNMR Spectra of MOFs. University of Ottawa, uOttawa Theses, 2015.

51

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Multiscale screening strategies aspire to go from a crystal structure to material performance in a real adsorption process 177x64mm (150 x 150 DPI)

ACS Paragon Plus Environment

Page 52 of 52