Monte Carlo Simulation of Near- and Supercritical ... - ACS Publications

Subscriber access provided by Northern Illinois University

Article

Monte Carlo Simulation of Near- and Super-Critical Hexane Fluid and Physisorption Phase Behavior Kenneth M. Benjamin, Alireza Asiaee, Carrie Veer, Casey Losinski, Samuel Gunderson, and Trevor Larson J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b03108 • Publication Date (Web): 29 Aug 2016 Downloaded from http://pubs.acs.org on August 30, 2016

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry C 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 29

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

The Journal of Physical Chemistry

“Monte Carlo Simulation of Near- and SuperCritical Hexane Fluid and Physisorption Phase Behavior”

Kenneth M. Benjamin1,* , Alireza Asiaee1, Carrie Veer1, Casey Losinski2, Samuel Gunderson1, Trevor Larson3 1

Department of Chemical and Biological Engineering, South Dakota School of Mines and Technology, 501 E. Saint Joseph St., Rapid City, SD 57701, USA 2

Cargill, Inc., Lafayette, Indiana

3

3M, Decatur, Alabama

1 ACS Paragon Plus Environment


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 2 of 29

Abstract Monte Carlo (MC) simulations have been used to determine the near- and super-critical fluid phase behavior of hexane, and the adsorption thermodynamics of hexane on cobalt. Canonical ensemble MC simulations have been used to compute the equation of state and chemical potentials, and Gibbs ensemble MC simulations (GEMC) have been used to compute fluid phase coexistence. Overall, MC results with the TraPPE-UA force field were successful in representing near- and super-critical hexane thermodynamic behavior, and represent an improvement over treatments with classical cubic equations of state.

In addition, GEMC

simulations predict a critical point for TraPPE-UA hexane of 509 K, 35 bar, and 0.22 g/ml, in excellent agreement with critical constants for real hexane. With regard to physisorption, hexane physisorption onto a (0001) cobalt surface at 523.15 K was modeled using grand canonical Monte Carlo (GCMC). GCMC excess adsorption results show crossover from adsorption to depletion around a bulk density of 0.1 g/ml, and a global minimum in depletion at 0.44 g/ml. Energetically, GCMC results indicate that the energy of adsorption decreases with increasing bulk density. Both molecule loading and adsorption energy data suggest three distinct adsorption regions: 1.) a low density (vapor-like) region, 2.) a near-critical region with little to no change in loading or energy of adsorption, and 3.) a high density (liquid-like) region.


Page 3 of 29

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


I. Introduction Adsorption and catalysis in supercritical fluids (SCFs) has been a subject of investigation for almost thirty years. Previous SCF adsorption studies have focused either on general physisorption characteristics or separation in porous media, often for permanent gases such as CH4 and N2.1,2 In the area of catalysis, SCFs have been employed as solvent environments for a number of chemistries, including oxidation, hydrothermolysis, and amination.3,4 Particularly, supercritical hexane has been found to be an excellent solvent media for Fischer-Tropsch (FT) synthesis, as part of an overall gas-to-liquids conversion technology.5-7 Elbashir et al. report a number of advantages to SCF FT synthesis over conventional gas- or liquid-phase FT synthesis, including improved heat transfer, reduced coking/deactivation of catalyst, and improved product distribution.7 While much experimental work has been conducted on catalytic SCF systems, very little attention has been given to computational studies. In particular, there have been no studies directed at SCF adsorption of substances which act as solvent environments in heterogeneous catalytic systems. In this regard, computational investigations, especially atomistic level modeling studies, hold much potential to provide molecular-level details on catalytic SCF systems, given their inherent microheterogeneity.8 In this work, we examine the fluid phase thermodynamics of near- and super-critical hexane, and the adsorption thermodynamics of SCF hexane on a model cobalt catalyst surface. This work is the first of its kind to examine physisorption of a realistic SCF solvent molecule (hexane) on a realistic catalyst (cobalt). The results presented herein contribute to the general understanding of physical chemistry of surface events in the presence of SCFs, with particular relevance to adsorption and catalysis in SCFs, especially catalytic SCF FT synthesis.



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

II. Background Near- and super-critical fluids are fluids that are heated and compressed near or beyond their critical temperature and pressure. In the vicinity of the critical point, the properties of a fluid vary greatly. Among these “tunable” properties are density, heat capacity, and isothermal compressibility.9 At the molecular level, the general phenomena of clustering and microheterogeneity in SCF environments likewise has been known for almost thirty years.8,10-13 However, some recent simulation studies have shown that the magnitude of clustering depends on the nature of the solute/solvent systems.14,15 Scientific and engineering applications involving SCFs often look to exploit this tunable property feature and unique solvent microstructure to optimize system performance, often with regard to solubility, mass transfer, or reactivity.13,16 The influence of SCF conditions on the physisorption of relatively permanent gases (such as CH4 and N2) has been fairly well studied. These studies include experimental, computational, and theoretical studies.1,2,17,18 In particular, Donohue proposed a new classification of IUPAC adsorption isotherms for SCF systems, based on the nature of the fluid-solid interaction (weak or strong).2 In this regard, most published work to date involves cases of relatively strong fluidsolid interaction, such as the case of CH4 or N2 adsorption onto activated carbon. In all of these cases of strong (that is, similar) fluid-solid interactions, investigators noted enhanced (or positive excess) adsorption under SCF conditions. At the molecular level, this likely relates to preferential adsorbate-adsorbent interactions, which leads to clustering of adsorbate molecules at the surface, with adsorbed layer densities greater than the bulk fluid density limit. In the most related to work to catalysis in SCFs, Asiaee and Benjamin have computed the potential of mean force for CO chemisorption onto a cobalt catalyst in the presence of SCF hexane, and determined the SCF solvent effects on adsorption energy and free energy.19 Specifically, these authors have 4 ACS Paragon Plus Environment

Page 4 of 29

Page 5 of 29

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


shown that SCF hexane tends to stabilize CO on the cobalt catalyst surface, relative to the ideal gas adsorption limit. With regard to physisorption, there have been a few recent reports of depletion (that is, negative excess adsorption) under SCF fluid conditions.20,21 For the simulated case of a SCF Lennard-Jones fluid onto a surface with a weak fluid-solid interaction, Oleinikova and Brovchenko noted a crossover from enhanced (positive excess) adsorption to depletion (negative excess adsorption) as one traverses from the ideal gas limit to SCF conditions in the bulk fluid phase.21 Moreover, this work highlighted the relative lack of information on physisorption of SCFs on surfaces with relatively weak fluid-solid interaction. In the case of heterogeneous catalysis in the presence of a SCF solvent, the extent of adsorption (or depletion) and characteristics within the adsorbed layer are extremely important, and likely determine the magnitude of the solvent effect on surface chemistry and physics. Experimental studies of Fischer-Tropsch synthesis in SCF hexane have suggested just such solvent effects on elementary catalytic kinetics and mechanisms.7

III. Computational Details III.A. Models Hexane is modeled with the united-atom transferable potentials for phase equilibria (TraPPE-UA) force field.22 This pairwise, united atom force field treats the CH3 and CH2 groups of hexane as united atom entities. Moreover, this is a semi-rigid model, meaning bond lengths are maintained at their equilibrium distances, while the model is flexible with regard to bond



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 6 of 29

angle bending and torsion. This is a reasonable approximation for fluid and adsorption phase thermodynamic investigations, as they are dominated more by intermolecular forces, rather than intramolecular forces (such as bond stretching). Intermolecular forces described by the TraPPE family of force fields include both van der Waals and electrostatic forces. Van der Waals interactions are modeled with a Lennard-Jones (LJ) potential:

(1) where σ and ε are the collision diameter and well depth, respectively. For the TraPPE-UA family of force fields, electrostatics are modeled with Coulomb’s law and point charges:

(2)

However, for the case of nonpolar hexane, there are no point charges (q values of zero) in the TraPPE-UA model. Therefore, there are no electrostatic forces between the fluid phase hexane molecules and the cobalt atoms of the catalyst surface. Intramolecular bond angle bending is described with a harmonic potential: ⁄

(3)

Torsional interactions are modeled with an OPLS functional form23:


(4)

Page 7 of 29

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


Cobalt atoms are modeled with the ReaxFF potential created by Labrosse et al.24 The ReaxFF functional form is designed to capture both the physics and chemistry existing between interacting atoms, and its complete development and form has been reported elsewhere.25 While ReaxFF is able to capture varying extents of chemical bonding between atoms, our simulations inherently are non-reactive, and strictly demonstrate the physical adsorption of hexane to a cobalt surface. The ReaxFF functional type was not supported in the Monte Carlo code used in this study. (More details on the Monte Carlo methods can be found in the following subsection.) Therefore, the ReaxFF potential for cobalt was re-cast into a Lennard-Jones form, to make it compatible with the TraPPE-UA potential used for hexane used within our MC simulations (described above). The LJ functional form was fitted to ReaxFF generated van der Waals potential energy data, and the deviation of the fitted LJ potential energy curve was minimized to determine the appropriate LJ σ and ε values. Figure 1 shows a comparison and the goodness of fit (R2 = 0.985) between the original ReaxFF potential and the fitted LJ potential for cobalt, and Table 1 reports the LJ parameters for Co (and united atom CH3 and CH2 for hexane) used in this study.



0.90 E-ReaxFF

0.70

E-LJ

0.50

Energy (kJ/mole)

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

0.30 0.10 -0.10 0

1

2

3

4

5

6

-0.30 -0.50

Distance (Angstroms)

-0.70 -0.90

Figure 1: ReaxFF vs. LJ potentials for Co-Co interaction. The solid line is the original van der Walls portion of the ReaxFF potential from Labrosse24, and the dashed line is the Lennard-Jones best-fit curve for implementation within MC simulations in this study.

atom

σ (Å)

ε (kJ/mol)

Co

1.564

0.7640

CH3

3.750

0.8146

CH2

3.950

0.3824

Table 1: Lennard-Jones potential parameters for atom types used in MC simulations. Cobalt LJ values regressed from the ReaxFF potential of Labrosse24, TraPPE-UA values for hexane from Martin and Siepmann22.

A potential cutoff of 2.5σ was used for all interactions, and analytical tail corrections were applied. Lorentz-Berthelot combining rules are used for all cross-interactions.26

III.B. Methods Monte Carlo (MC) simulations are conducted with the Towhee Monte Carlo program.27 A series of different statistical mechanical ensembles were used to compute properties of single


Page 8 of 29

Page 9 of 29

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


phases SCF hexane, two-phase high pressure hexane, and hexane physisorption on cobalt. These different ensembles and simulations are represented in Figure 2 and described below. Adsor pt i on Ph ase Equ i l i br i u m

Bulk Fluid Behavior (PVT + µ) Canonical and Grand Canonical Ensemble Monte Carlo (NVT, µVT)

Fluid Phase Hexane

Phase Equilibrium Grand Canonical Ensemble Monte Carlo (µVT)

Adsorbed Phase Hexane

(Determine loading)

Catalyst Surface

Figure 2: Simulation approach for adsorption of hexane under SCF conditions.

III.B.i. Single Phase SCF Behavior Canonical (NVT) and grand canonical (µVT) ensemble MC simulations were conducted to determine the equation of state and chemical potentials of SCF hexane, at a supercritical temperature of 523 K (above the critical point of real and TraPPE-UA hexane; see sections III.A and IV.B for more information). NVT simulations contained 64 hexane molecules, with volumes varying to span a density range from 0.15 to 0.60 g/ml. (Simulations were also conducted with 128 and 256 hexane molecule systems, and confirmed the absence of finite size effects.) Fluid pressures were computed with the standard virial equation.28 MC moves utilized included center-of-mass translation, center-of-mass rotation, atom translation (to allow bond angle bending and torsional sampling), and particle regrowth, based on the configurational bias MC method.29 NVT simulations were run for 6,250 MC cycles for equilibration and over 71,000 MC cycles for production. Block averaging was used to quantify uncertainty in computed values. 9 ACS Paragon Plus Environment


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 29

Chemical potentials of SCF hexane thermodynamic states were determined from a series of µVT ensemble MC simulations. A cubic simulation box with a side length of 30 Å was used to set the simulation volume, and the temperature was fixed at 523 K. In addition to the MC moves listed above (for the NVT simulations), configurational biased particle insertion/deletion moves were added. A range of chemical potentials from -500 to -1400 K were fixed as constraints during these simulations, and fluid phase densities and pressures were computed. Equilibration and production were conducted for at least 2,100 and 40,000 MC cycles, respectively. The PVT behavior determined from these µVT simulations agreed with the behavior from the original NVT simulations. Moreover, these simulations provided confirmation of the correct chemical potential at each bulk density to be used in subsequent physisorption simulations.

III.B. ii. High Pressure Phase Behavior Proper analysis and interpretation of near- and super-critical fluid phenomena requires knowledge of the critical point of the material. In this regard, the critical point of TraPPE-UA hexane is unknown.

Therefore, Gibbs ensemble Monte Carlo (GEMC) were conducted to

compute the TraPPE-UA hexane critical point. More information on the GEMC method and its use for phase equilibrium calculations can be found elsewhere.30 GEMC simulations were run as two-box NVT simulations, at temperatures between 450-480 K, and with a total of 200 molecules. In addition to the MC moves used for the single-phase NVT simulations (see section III.B.i), configurational biased particle insertion and deletion MC moves were added. In addition to coexistence densities in the separate boxes, saturation pressures were determined from the


Page 11 of 29

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


virial equation, as applied to the vapor phase box. Equilibration and production were run for 5,000 and 35,000 MC cycles, respectively. It is well-known that simulating too close to the critical point is difficult and prone to error due to extreme fluctuations, as represented by the molecular clusters and voids discussed previously. Recently, Dinpajooh et al. has provided several insights into the proper use of GEMC simulations to predict very accurate and precise critical properties.31 One of the major revelations of this work was the realization that finite size effects are relatively small in GEMC simulations, as compared to other techniques, and that smaller systems with shorter cutoff distances still provide reliable critical point data. Given this, rather than simulating much larger systems with larger cutoff distances, we have chosen to use scaling laws and associated relationships to estimate the critical point. These scaling laws will be presented and discussed during the corresponding results section.

III.B.iii. Hexane Physisorption Hexane physisorption onto Co(0001) was studied using grand canonical (µVT) Monte Carlo (GCMC) methods. Bulk chemical potentials computed from the methods outlined in section III.B.i are used as constraints for the physisorption GCMC simulations. A 10x10 Co(0001) slab with four layers is used for the catalyst surface. A picture of the catalyst surface is shown in Figure 3. The surface has an overlayer area of (20Å x 20Å)R120°. An adsorbed layer (vacuum) height of 20 Å was used, with a hard reflective wall placed at the top of the simulation box. Further, the height of the hard wall was set to 300 Å to prevent any interaction between neighboring catalyst surfaces, thereby further representing adsorption to a free surface 11 ACS Paragon Plus Environment


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

(versus a nano-confined arrangement). (Select points were also simulated with a height of 40 Å, to confirm the validity of the 20 Å adsorption loading results.) During these simulations, the cobalt atoms were held at fixed positions, and not permitted to move. The initial coordination of cobalt atoms was determined from a density functional theory (DFT) calculation in order to find the minimum level of energy. DFT calculations were carried out by the DMol3 method32 of BIOVIA Materials Studio via the Perdew-Wang generalized-gradient approximation (GGAPW91).33 Additional details of the DFT calculation can be found in the literature.19 MC moves sampled included center-of-mass translation, atom translation (to sample angle bending and torsion), center-of-mass rotation, and configurational biased insertion/deletion (to allow satisfaction of the chemical potential constraint). Equilibration and production were run for at least 8,300 and 200,000 MC cycles, respectively. (Simulations were also conducted with a system four times larger in catalyst area than the one described here, and confirmed the absence of finite size effects on computed loadings and energetics.)

Figure 3: Top view of Co(0001) surface.


Page 12 of 29

Page 13 of 29

IV. Results and Discussion IV.A. High Pressure Phase Equilibrium Successful interpretation of SCF hexane physisorption results requires accurate prediction of high pressure hexane properties, including critical properties. As mentioned in section III.B.ii, GEMC simulations were used to compute high pressure hexane phase equilibrium and to estimate a critical point for TraPPE-UA hexane. Figure 4 compares GEMC results with TraPPE-UA to experimental data from NIST.34 The uncertainty in computed points (error bars not shown) is between 4-6% for saturated liquid densities and between 14-18% for saturated vapor densities. As the figure shows, the agreement between simulation and experiment is excellent. This result is to be expected, given the ability of TraPPE-UA to successfully reproduce phase diagrams for other hydrocarbons.22 480

Temperature (K)

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


475 470 465 460 455 450 0.00

0.10

0.20

0.30

0.40

0.50

Density (g/ml)

Figure 4: Coexistence curve for high pressure hexane. Open circles are GEMC simulation data and filled circles are data from NIST.34

Due to excessive fluctuations approaching the critical point, we were unable to successfully simulate phase coexistence above 480 K. Therefore, we have used a scaling law to estimate the critical temperature for TraPPE-UA hexane35: !"# $ % % &

(5) 13 ACS Paragon Plus Environment


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 29

In equation (5) above, a critical exponent of β = 0.325 is used, which is common for hydrocarbon fluids. The critical density is computed from the law of rectilinear diameters36

' !"# ( ) % %

(6)

Finally, the critical pressure is determined from the Clausius-Clapeyron equation37: *+ ," -

-.

(7)

%

From this post-analysis of the GEMC data, the critical constants for TraPPE-UA hexane are shown in Table 2. The critical temperature of 509 K calculated here is in good agreement with the one provided previously using the same forcefield.38 As one notes, the TraPPE-UA force field does an excellent job of reproducing the critical temperature and density (within 6%), and has reasonable agreement with the experimental critical pressure (within 18%). This type of over-prediction of critical pressures is common for pairwise potential models, such as TraPPEUA. NIST34

GEMC/TraPPE

% error

(this study) Tc (K)

507.82

50923

0.23

Pc (bar)

30.34

352

18

ρc (g/ml)

0.2332

0.21911

5.9

Table 2: Critical properties of hexane – real hexane34 versus GEMC/TraPPE-UA values from this study. Numbers in subscript represent uncertainties (+/-) in computed GEMC properties.


Page 15 of 29

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


Beyond general property information, the confirmation of TraPPE-UA critical temperature provides reassurance that temperatures above 509 K are truly supercritical, with regard to simulated TraPPE-UA hexane bulk and adsorption behavior.

IV.B. SCF Single Phase Results Given the confirmation regarding high pressure phase behavior and the critical point for TraPPE-UA hexane, we now turn our attention to single phase supercritical TraPPE-UA hexane bulk fluid behavior. In particular, we seek to explore behavior along an isotherm of 523 K, which is a common supercritical reaction temperature for SCF hexane FT synthesis. Figure 5 displays pressure versus density data for hexane along a 523 K isotherm, from both NVT-MC simulations and NIST.34 In this figure, the uncertainties in the computed pressures range from 316% over the specified density range. As the figure shows, the agreement here for single phase PVT properties is also excellent. For further comparison, Figure 5 also shows PVT predictions from the Peng-Robinson equation of state (EOS), a commonly used cubic EOS in engineering calculations.39 From the figure, it is clear that NVT-MC predictions with TraPPE-UA hexane represent a modest improvement over such empirical cubic EOS for the prediction of SCF PVT behavior. (Please note that the R2 values for the goodness of fit for the PREOS and MC results to the NIST data are 0.98 and 0.99, respectively.)


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

Pressure (kPa)


90000 80000 70000 60000 50000 40000 30000 20000 10000 0 0.15

MC NIST PREOS

0.25

0.35

0.45

0.55

0.65

Density (g/ml)

Figure 5: Hexane isotherm at 523 K. Open circles are MC simulation data, filled circles are data from NIST34, and the line is from the Peng-Robinson equation of state39.

According to the simulation approach outlined in Figure 2, the chemical potentials of bulk SCF hexane are required, as constraints for subsequent µVT physisorption simulations. Therefore, GCMC simulations of bulk hexane at 523 K were conducted, as outlined in section III.B.iii. For these simulations, chemical potentials were set as constraints, and bulk densities were computed from the ensemble averaged number of molecules. These µVT results are presented in Figure 6, in which they are simultaneously compared to NIST hexane results.34 Once again, one sees the suitability of TraPPE-UA hexane to reproduce the chemical potential behavior of real hexane as a function of density at 523 K. These results now allow us to map the subsequent physisorption behavior to the bulk fluid thermodynamic state as a function of bulk temperature and density (or pressure), rather than temperature and chemical potential. This is more convenient, both from an experimental and industrial operations standpoint.


Page 16 of 29

Page 17 of 29

0.6 0.5

Density (g/ml)

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


0.4 0.3 0.2 0.1 0 -1600

-1400

-1200

-1000

-800

-600

-400

Chemical Potential (K)

Figure 6: Density versus chemical potential for hexane at 523 K. Filled circles are GCMC simulation data and open circles are data from NIST.34

IV.C. Physisorption Loading In examining physisorption, it is often common to analyze adsorption loading results as a function of bulk density and/or pressure, rather than chemical potential. Here, we present and discuss both cases. The results for physisorption loading versus bulk density at 523 K is shown in Figure 7. Here, one notes that loading increases with bulk density from the ideal gas limit (zero density) up to the critical density of around 0.2 g/ml. This type of adsorption behavior from the gas phase is to be expected. However, between 0.2-0.4 g/ml bulk density, there is essentially no change in the adsorption loading. This behavior is very much consistent with the limited change in chemical potential as a function of density near the critical point, as shown previously in Figure 6. On a molecular level, this limited change in adsorption loading is likely evidence of preferential clustering between hexane molecules in the near-critical fluid state, as opposed to hexane-cobalt clustering at the catalyst surface. At densities larger than 0.4 g/ml and up to 0.6 g/ml, hexane loading again increases with increasing bulk density. Again, this type of adsorption behavior would be expected in a larger density, liquid-like regime. 17 ACS Paragon Plus Environment


0.05

Loading (molecules/A^2)

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

0.04 0.03 0.02 0.01 0 0

0.1

0.2

0.3

0.4

0.5

0.6

Bulk Density (g/ml)

Figure 7: Hexane physisorption loading versus bulk hexane density at 523 K.

The same GCMC physisorption data is re-plotted in Figure 8, but this time as a function of bulk pressure. In this case, bulk pressures were mapped to bulk densities from the NVT-MC results presented in section IV.B. When analyzed in this manner, we again see two distinct adsorption regions, gas-like and liquid-like, separated by the computed TraPPE-UA hexane critical pressure of 35 bar. While loading increases with increasing bulk pressure in both of these regions, the rates of increase (slope) are different in the two regions, with a larger rate/slope in the gas phase than in the liquid phase. Generally, for systems of gas adsorption on porous solids with strong fluid-solid interaction, a maximum appears in the plot of loading versus pressure (due to the divergent behavior of the fluid isothermal compressibility).18 However, the free cobalt surface (not nano-confined) modeled in the current study is more similar to the mesoporous adsorbent categories, whose corresponding plots have no maximum point.2 It should be noted that the uncertainties in the computed loadings presented in Figures 7 and 8 range between 3-6%.


Page 18 of 29

Page 19 of 29

0.045


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


0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0.000 0

50

100

150

200

250

Bulk Pressure (bar)

Figure 8: Hexane physisorption loading versus bulk hexane pressure at 523 K.

Beyond straightforward adsorption isotherms, another manner of interpreting adsorption loadings is through the concept of excess adsorption. Excess adsorption is defined as follows: /0 )# 2

" 34

(8)

4

where ρadsorbed and ρbulk are the densities of the adsorbed and bulk environments, respectively. Figure 9 presents excess adsorption as a function of bulk density. Here we see that at densities less than 0.1 g/ml, excess adsorption is positive and decreases with increasing bulk density. At densities between 0.1-0.4 g/ml, excess adsorption becomes negative and decreases with increasing bulk density. This phenomenon of negative excess adsorption is called “depletion” and has been observed before both experimentally and through molecular simulation.20,21 In particular, Oleinikova and Brovchenko observed a similar crossover from adsorption to depletion for a model Lennard-Jones fluid on a weak interacting surface.21 Such a large negative excess adsorption in this near-critical region again suggests that hexane molecules prefer interaction with themselves over interaction with the cobalt catalyst atoms, relatively. This same preference 19 ACS Paragon Plus Environment


for adsorbate-adsorbate clustering/interaction, and simultaneously weak adsorbate-adsorbent interactions, in the supercritical state was observed by Rother et al. as well, for propane physisorption onto a silica aerogel.20 Finally, at bulk densities greater than 0.4 g/ml, excess adsorption remains negative, but increases with increasing densities. 40%

Excess Adsorption

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 20 of 29

20%

Bulk Density (g/ml)

0% 0

0.1

0.2

0.3

0.4

0.5

0.6

-20% -40% -60% -80%

Figure 9: Excess hexane adsorption versus bulk hexane density at 523 K.

IV.D. Physisorption Energetics In addition to the effects of SCF conditions on hexane loading, it is also interesting to explore the energetic aspects of adsorption in SCF conditions. To begin this discussion, we present the total (internal) energy of adsorption as a function of bulk density, as determined from ensemble-averaged energies of our MC simulations. Here, the total energy of adsorption is calculated by: ∆/%" /"674 /" /74

(9)

The dependence of total internal energy of adsorption on density is shown in Figure 10, where the uncertainties in calculated energies range from 5-18%. Presented in this manner, it suggests that the energy of adsorption increases from exothermic to endothermic with increasing bulk 20 ACS Paragon Plus Environment

Page 21 of 29

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


density, which seems counterintuitive and incorrect. To examine these results deeper, we consider the thermodynamic cycle shown in Figure 11. This cycle presents an alternative path to “direct” adsorption, represented by the quantity ∆Etotal. In this alternative path, an individual hexane molecule is “removed” from a bulk hexane environment in the first step, and then subsequently adsorbed on the catalyst in the second step. This cycle is represented by the expression: ∆/%" ∆/89! ∆/)

(10)

By this alternative thermodynamic cycle, one can compute the true, separate energy of adsorption (∆Eadsorb) as follows: ∆/) ∆/%" ∆/89!

(11)

The quantities ∆Etotal are already known, as presented in Figure 10. The quantity ∆Eremove is computed from: ∆/9! ∆/: 94 ∆/4 74

(12)

Here, ∆Ebulk fluid is the ensemble averaged energy from the NVT-MC bulk simulations (sections III.B.i and IV.B), and ∆Esingle molecule is the ensemble-averaged energy of a single-molecule (ideal gas limit).


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

Energy of Adsorption (kJ/mol)


5 4 3 2 1 0 -1 0 -2 -3 -4 -5 -6

0.1

0.2

0.3

0.4

0.5

0.6

Bulk Density (g/ml)

Figure 10: Total energy of adsorption versus bulk hexane density at 523 K.

∆EAdsorb

∆ERemove

∆ETotal Cobalt Catalyst

Bulk SCF Hexane

Figure 11: Thermodynamic cycle for adsorption in supercritical fluid conditions.

These separate energetic contributions to adsorption are plotted in Figure 12 (specifically, the left hand graph). Here, one notes that the total energy of adsorption results from the counteracting effects of “removal” and SCF “adsorption”. As one contribution, ∆Eremove increases continuously with increasing bulk density, as one would expect, given its essential equivalence to the change in energy upon vaporization (molecule going from liquid to gas environment). However, for the specific contribution of adsorption, (∆Eadsorb), one notes that 22 ACS Paragon Plus Environment

Page 22 of 29

Page 23 of 29

overall this quantity is negative (exothermic) and decreases with increasing bulk density. Upon finer examination (from the right hand graph of Figure 12), we again see three distinct regions here. In the low density (< 0.2 g/ml) gas-like region, the energy of adsorption decreases with increasing bulk density. In the intermediate, near-critical density region (0.2-0.4 g/ml), the energy of adsorption remains nearly constant. This constant energy change is directly tied to the constant loading over this region (Figure 7), and is a direct consequence of the fact that no additional hexane molecules adsorb to the surface over this density range, likely due to their preferential interaction with themselves in the bulk over the cobalt catalyst atoms. At larger densities (> 0.4 g/ml), the energy of adsorption again decreases with increasing bulk density. In both the gas- and liquid-like regimes, the decrease in adsorption energy (∆Eadsorb) with increasing density suggests that the larger hexane bulk fluid densities help stabilize hexane molecules adsorbed to the surface, thereby lowering the energy.

18 16 14 12 10 8 6 4 2 0 -2 0 -4 -6 -8 -10 -12 -14

-6

DE(remove) DE(ads) DE(Total)

0.1

0.2

0.3

0.4

0.5

0.6

Adsorption Energy (kJ/mol)

Energy (kJ/mol)

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


-7 -8 -9 -10 -11 -12 -13 0

0.1

Bulk Density (g/ml)

0.2

0.3

0.4

0.5

Bulk Density (g/ml)

Figure 12: Energetic components of SCF adsorption versus bulk density at 523 K. Left.) All three steps in thermodynamic cycle, Right.) Enlarged view of ∆E(adsorption) behavior.


0.6


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

V. Conclusions Monte Carlo simulations demonstrate the suitability of the TraPPE-UA force field for reproducing the fluid phase thermodynamic properties of near- and super-critical hexane, including densities, chemical potentials, and phase equilibrium. In particular, GEMC simulations predict a critical point for TraPPE-UA hexane of 509 K, 35 bar, and 0.22 g/ml, in excellent agreement with critical constants for real hexane. With regard to hexane physisorption onto a (0001) cobalt surface at 523.15 K, GCMC excess adsorption results show crossover from adsorption to depletion around a bulk density of 0.1 g/ml, and a global minimum in depletion at 0.44 g/ml. Energetically, GCMC results indicate that the energy of adsorption decreases with increasing bulk density. Both molecule loading and energy data suggest three distinct adsorption regions: 1.) a low density (vapor-like) region, 2.) a near-critical region with little to no change in loading or energy of adsorption, and 3.) a high density (liquid-like) region. Upon consideration of the specific application of catalytic SCF hexane FT synthesis on a cobalt catalyst, it is now apparent that adsorbed layer hexane densities will be smaller than corresponding bulk densities (set by bulk operating temperature and pressure). This will result in fewer hexane molecules available in the adsorbed layer to solvate adsorbates and transition states during surface reactions, and will likely change reaction energetics (heats of reaction and activation energy barriers). Therefore, any future investigations of elementary surface kinetics within catalytic SCF FT, or any catalytic system involving relatively weak adsorbate-surface


Page 24 of 29

Page 25 of 29

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


interactions, must account for the extent of SCF solvent physisorption to properly capture solvent effects on heterogeneous catalytic reaction kinetics and pathways.

Supporting Information Original Monte Carlo simulation numerical data for Figures 1, 4-10, and 12.

Acknowledgements The authors would like to recognize Kyle Benjamin for assistance in running physisorption GCMC simulations and analyzing physisorption data.

Author Information * To whom correspondence should be addressed: [email protected]; 605-394-2636



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) Kelley, F. D.; Chimnowitz, E. H. Near-Critical Phenomena and Resolution in Supercritical Fluid Chromatography. AIChE J. 1990, 36, 1163-1175. (2) Donohue, M. D.; Aranovich, G. L. A New Classification of Isotherms for Gibbs Adsorption of Gases on Solids. Fluid Phase Equil. 1999, 158-160, 557-563. (3) Baiker, A. Supercritical Fluids in Heterogeneous Catalysis. Chem. Rev. 1999, 99, 453-474. (4) Savage, P. E. A Perspective on Catalysis in Sub- and Supercritical Water. J. Supercrit. Fluids 2009, 47, 407-414. (5) Yokota, K.; Hanakata, Y.; Fujimoto, K. Supercritical Phase Fischer-Tropsch Synthesis. Chem. Eng. Sci. 1990, 45, 2743-2749. (6) Bochniak, D. J.; Subramaniam, B. Fischer-Tropsch Synthesis in Near-Critical n-Hexane: Pressure-Tuning Effects. AIChE J. 1998, 44, 1889-1896. (7) Elbashir, N. O.; Bukur, D. B.; Durham, E.; Roberts, C. B. Advancement of Fischer-Tropsch Synthesis via Utilization of Supercritical Fluid Reaction Media, AIChE J. 2010, 56, 997-1015. (8) Tucker, S. C. Solvent Density Inhomogeneities in Supercritical Fluids. Chem. Rev. 1999, 99, 391-418. (9) Levelt Sengers, J. M. H. Thermodynamics of solutions near the solvent's critical point. In Supercritical Fluid Technology: Reviews in Modern Theory and Applications; Bruno, T. J. Ely, J. F., Eds.; CRC Press: Boca Raton, 1991; pp. 1-56.   (10) Debenedetti, P. G.; Mohamed, R. S. Attractive, Weakly Attractive, and Repulsive NearCritical Systems. J. Chem. Phys. 1989, 90, 4528. (11) Petsche, I. B.; Debenedetti, P. G. Solute-Solvent Interactions in Infinitely Dilute Supercritical Mixtures: A Molecular Dynamics Investigation. J. Chem. Phys. 1989, 91, 7075. (12) Debenedetti, P. G.; Petsche, I. B.; Mohamed, R. S. Clustering in Supercritical Mixtures: Theory, Applications and Simulations. Fluid Phase Equil. 1989, 52, 347. (13) Eckert, C. A.; Knutson, B. L.; Debenedetti, P. G. Supercritical Fluids as Solvents for Chemical and Materials Processing. Nature 1996, 383, 313-318. (14) Stubbs, J. M.; Drake-Wilhelm, D. D.; Siepmann, J. I. Partial Molar Volume and Solvation Structure of Naphthalene in Supercritical Carbon Dioxide: A Monte Carlo Simulation Study. J. Phys. Chem. B 2005, 109, 19885-19892. (15) Anderson, K. E.; Siepmann, J. I. Solubility in Supercritical Carbon Dioxide: Importance of the Poynting Correction and Entrainer Effects. J. Phys. Chem. B 2008, 112, 11374. (16) Savage, P. E.; Gopalan, S.; Mizan, T. I.; Martino, C. J.; Brock, E. E. Reactions at Supercritical Conditions: Applications and Fundamentals. AIChE J. 1995, 41, 1723-1778. 26 ACS Paragon Plus Environment

Page 26 of 29

Page 27 of 29

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


(17) Tan, Z.; Gubbins, K. E. Adsorption in carbon micropores at supercritical temperatures. J. Phys. Chem. 1990, 94, 6061-6069. (18) Findenegg, G. H.; Löring, R. Fluid Adsorption up to the Critical Point. Experimental Study of a Wetting Fluid/Solid Interface. J. Chem. Phys. 1984, 81, 3270-3276. (19) Asiaee, A.; Benjamin, K. M. Molecular Simulation of CO Chemisorption on Co(0001) in Presence of Supercritical Fluid Solvent: A Potential of Mean Force Study. J. Chem. Phys., in press. (20) Rother, G. Melnichenko, Y. B. Cole, D. R. Frielinghaus, H. Wignall, G. D. Microstructural Characterization of Adsorption and Depletion Regimes of Supercritical Fluids in Nanopores. J. Phys. Chem. C 2007, 111, 15736. (21) Oleinikova, A.; Brovchenko, I. Nonmonotonic Crossover from Adsorption to Desorption in Supercritical Fluid Near a Weakly Attractive Surface. Phys. Rev. E 2008, 78, 061601. (22) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569-2577. (23) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638-6646. (24) LaBrosse, M. Statistical Mechanical and Quantum Mechanical Modeling of Condensed Phase Systems; Ph.D. Thesis: University of Pittsburgh, 2009. (25) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396-9409. (26) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (27) Martin, M. G. MCCCS Towhee: A Tool for Monte Carlo Molecular Simulation. Mol. Sim. 2013, 39, 1212-1222; http://towhee.sourceforge.net (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press; Oxford; 1987.  (29) Siepmann, J. I.; Frenkel, D. Configurational Bias Monte Carlo: A New Sampling Scheme for Flexible Chains. Mol. Phys. 1992, 75, 59-70. (30) Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble. Mol. Phys. 1987, 61, 813-826. (31) Dinpajooh, M.; Bai, P.; Allan, D. A.; Siepmann, J. I. Accurate and Precise Determination of Critical Properties from Gibbs Ensemble Monte Carlo Simulations. J. Chem. Phys. 2015, 143, 114113. (32) Delley, B. From Molecules to Solids with the DMOL3 Approach. J. Chem. Phys. 2000, 113, 7756.



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

(33) Perdew, J. P.; Wang, Y. Accurate and Simple Analytic Representation of the Electron-Gas Correlation Energy. Phys. Rev. B 1992, 45, 13244. (34) Span, R. Multiparameter Equations of State – An Accurate Source of Thermodynamic Property Data; Springer: Berlin, 2000. (35) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Oxford University Press: New York, 1989.  (36) Rowlinson, J. S.; Swinton, F. L. Liquids and Liquid Mixtures; 3rd ed.; Butterworth: London, 1982.  (37) Atkins, P. W. Physical Chemistry; 4th ed.; W.H. Freeman and Company: New York, 1990. (38) Potoff, J. J.; Siepmann, J. I. Vapor-Liquid Phase Equilibria for Linear and Branched Alkane Monolayers Physisorbed on Au(111). Langmuir 2002, 18, 6088-6095.   (39) Peng, D. Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fund. 1976, 15, 59-64.


Page 28 of 29

Page 29 of 29

Table of Contents (TOC) Image

-6

Adsorption Energy (kJ/mol)

0.05


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


0.04

0.03

0.02

0.01

-7 -8 -9 -10 -11 -12

0 0

0.1

0.2

0.3

0.4

0.5

-13

0.6

0

0.1

Bulk Density (g/ml)

0.2

0.3

0.4

Bulk Density (g/ml)


0.5

0.6

Monte Carlo Simulation of Near- and Supercritical ... - ACS Publications

Recommend Documents