Development of the Transferable Potentials for Phase Equilibria Model

May 18, 2015 - •S Supporting Information. ABSTRACT: The transferable potentials for phase equilibria force field is extended to hydrogen sulfide...
0 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCB

Development of the Transferable Potentials for Phase Equilibria Model for Hydrogen Sulfide Mansi S. Shah,†,§ Michael Tsapatsis,† and J. Ilja Siepmann*,†,‡,§ †

Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Avenue SE, Minneapolis, Minnesota 55455-0132, United States ‡ Department of Chemistry and §Chemical Theory Center, University of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *

ABSTRACT: The transferable potentials for phase equilibria force field is extended to hydrogen sulfide. The pure-component and binary vapor−liquid equilibria with methane and carbon dioxide and the liquidphase relative permittivity are used for the parametrization of the Lennard−Jones (LJ) and Coulomb interactions, and models with three and four interaction sites are considered. For the three-site models, partial point charges are placed on the sites representing the three atoms, while the negative partial charge is moved to an off-atom site for the four-site models. The effect of molecular shape is probed using either only a single LJ interaction site on the sulfur atom or adding sites also on the hydrogen atoms. This procedure results in four distinct models, but only those with three LJ sites can accurately reproduce all properties considered for the parametrization. These two are further assessed for predictions of the liquid-phase structure, the lattice parameters and relative permittivity for the face-centered-cubic solid, and the triple point. An effective balance between LJ interactions and the dipolar and quadrupolar terms of the first-order electrostatic interactions is struck in order to obtain a four-site model that describes the condensed-phase properties and the phase equilibria with high accuracy.

■

level of theory. However, this force field results in a significant underestimation of the vapor pressure.5 Subsequently in 1997, Kristóf and Liszi5 reparametrized the four-site model of Forester to accurately reproduce the vapor−liquid coexistence curve (VLCC) of H2S. The liquid-phase RP for this model at T = 212 K is 16.1 compared to the experimental value of 8.04.8 Also, this model overestimates the strength of the interactions between H2S and CO2 when binary phase equilibria are simulated. In 2000, Delhommelle et al.9 stressed the need for a polarizable force field for H2S to better predict binary phase equilibria for polar/nonpolar mixtures, and proposed a five-site polarizable model for H2S to improve pressure−composition curves for binary mixtures of H2S and n-pentane. In an attempt to develop an atom-based force field, Nath developed a threesite model.10 This model predicts with good accuracy the phase equilibria of H2S with higher alkanes, but it underestimates the boiling point by 7 K. Kamath et al.11 proposed four models in 2005, each of them being a three-site model with varying values of the partial charges. RP values at T = 194.6 K for their models A, B, C, and D are 8.2, 20.4, 30, and 37, respectively, while the experimental value is 8.998 at this temperature. In a subsequent paper by Kamath and Potoff,12 the deficiency of their H2S model to correctly estimate the interactions between CO2 and

INTRODUCTION Hydrogen sulfide (H2S) is a very hazardous compound; its safe and efficient handling poses a tremendous challenge to oil and gas industries. For large-scale sweetening of natural gas, highly energy-intensive amine-based absorption processes are employed.1 The resulting acid gas stream is either reinjected into reservoirs or used for sulfur recovery using the Claus process.2 In order to meet the emission regulations for the off gases from the Claus unit, expensive tail gas treatment such as subdewpoint processes or amine-based absorptive separation is employed.2 Innovative solutions for these separations are of immense environmental and economic interest. Considering the health hazards of H2S, with concentrations as low as 500 ppm being fatal,3 predictive molecular modeling may prove instrumental for providing molecular-level insights and for guiding the development of improved separation processes, but the success of such modeling studies hinges on the availability of accurate force fields to describe the systems of interest. Molecular simulations of systems containing H2S have been carried out for almost three decades. In 1986, Jorgensen4 proposed the first force field to describe H2S: a three-site model developed to reproduce its enthalpy of vaporization and liquid density at the normal boiling point. However, this model grossly overpredicts the vapor pressure5 and the liquid-phase relative permittivity (RP)6 of H2S. In 1989, Forester et al.7 introduced a four-site model by fitting to the dimer structure predicted by a distributed-multipole analysis at the 6-31 G* © 2015 American Chemical Society

Received: March 16, 2015 Revised: May 15, 2015 Published: May 18, 2015 7041

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B

the standard Lorentz−Berthelot combining rules are used for all unlike interactions:16

H2S has been highlighted. These authors mention that induced polarization effects are not expected to significantly improve the predicted phase behavior for the H2S/CO2 system because of the relatively small dipole moment of H2S. More recently, Drude-polarizable force fields for H2S have been proposed.6,13 In this work, we seek an H2S model for the Transferable Potentials for Phase Equilibria (TraPPE) force field that overcomes prior deficiencies and more accurately describes the interactions of H2S with CO2 and CH4 that are pivotal for modeling of sour gas mixtures. Targets for our force field include the following: reproducing the experimental critical temperature of H2S within 1%, its saturated liquid density at the normal boiling point within 1%, its saturated vapor pressures over the complete VLCC range within 10%, and the compositions in binary vapor−liquid equilibria (VLE) with CO2 and CH4 within an average deviation of less than 5% using standard combining rules for unlike interactions. To our knowledge, none of the available force fields in the literature meet all of these constraints. Furthermore, the performance of the resulting model is also assessed via other condensed-phase properties and the prediction of the triple point temperature.

σij =

FORCE FIELDS The H2S models developed and assessed in this work and also the TraPPE all-atom models14,15 for CH4 and CO2 are nonpolarizable and have a rigid geometry. The internal geometry of the H2S model is the same as that of the model by Jorgensen.4 Bond lengths and bond angles for molecules investigated here are listed in Table 1. Nonbonded interactions Table 1. Bond Lengths and Bond Angles for TraPPE Models bond type

length [Å]

angle type

angle [deg]

H2S CH4 CO2

rH−S rC−H rC−O

1.34 1.10 1.16

∠HSH ∠HCH ∠OCO

92 109.4712 180

are modeled using a pairwise-additive potential consisting of Lennard−Jones (LJ) 12−6 and Coulomb terms: ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj σij σij U (rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ + r ⎝ rij ⎠ ⎦ 4πε0rij ⎣⎝ ij ⎠

2

and εij =

εiiεjj

(2)

Model Development. This subsection describes four different models for H2S that are explored in our effort to extend the TraPPE force field to H2S. Broadly, two variations and their combinations are explored: (i) the charge distribution and (ii) the molecular shape as defined by the LJ interactions. Shifting the partial negative charge from the location of the S atom (where it is placed for models with three interaction sites) to an additional off-atom X-site on the H−S−H angle bisector and placed toward the H atoms (for models with four interaction sites) allows one to explore the effects of changing the charge distribution and adjusting the ratio of dipole to quadrupole moments of H2S. As far as changing the molecular shape through the LJ interactions is concerned, additional LJ sites can be placed on the H atoms in addition to the S atom to account for the nonspherical shape of the H2S molecule. Exploring the two variations results in four different model types, which are illustrated in Figure 1. Here, type A-B corresponds to a model with A interaction sites and B LJ interaction centers. For instance, type 4-3 implies four interaction sites with three LJ interaction sites placed on the atomic positions and three partial charges placed on the locations of the H atoms and the off-atom X site. This 4-3 model requires fitting of four LJ parameters (for the sites placed on H and S atoms), one X-site displacement (for the distance between S atom and X site), and the partial negative charge on the X site (where the constraint of molecular neutrality requires qX = −2 qH). In principle, one could also construct four-site models with partial charges and/or LJ sites placed on all four sites. This would significantly expand the parameter space because it would allow for two adjustable parameters for partial charges and two additional LJ parameters, but such models are not considered in order to limit the already high dimensional parameter space (6-dimensional for the 4-3 model) to be explored for a wide range of properties. Model Type 3-1. Type 3-1 is the simplest of all models and consists of three sites with partial charges and a single LJ site located at the sulfur atom. In this case, there are three parameters (εSS, σSS, and qS) that need to be optimized. Screening models for the liquid-phase RP reveals that this property is not very sensitive to the LJ parameters (as long as the density is close to the experimental value), but is almost entirely determined by the partial charges. A qS value between −0.25 and −0.26 |e| results in a reasonable RP for the model. Having settled on a value for qS, LJ parameters are optimized to

■

molecule

σii + σjj

(1)

where rij, εij, σij, qi, and qj are the site−site separation, LJ well depth, LJ diameter, and partial charges for beads i and j, respectively. In order to allow applicability over a wide spectrum of potential compositions of natural gas streams and also for the sake of consistency with other TraPPE models,

Figure 1. Schematic representation of the types of H2S models. Type A-B implies an A-site model with B LJ interaction sites. Yellow, black, and cyan spheres represent the S atom, H atoms, and the off-atom X site, respectively. 7042

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B

Table 2. Force Field Parameters: Partial Charges (qi), LJ Parameters (εii and σii), and Displacement (δS−X) of Off-Atom X Site molecule

reference

model

qS [|e|]

H2S H2S H2S H2S H2S H2S H2S

Kamath et al.11 this work Kristóf and Liszi5 this work Nath10 this work this work

3-1 3-1 4*-1 4-1 3-3 3-3 4-3

−0.252 −0.25 +0.40

CH4

TraPPE−EH14

CO2

TraPPE15

qX [|e|]

εSS/kB [K]

σSS [Å]

278 278 250 248 250 125 122 εCC/kB 0.01 εCC/kB 27

3.71 3.71 3.73 3.75 3.72 3.60 3.60 σCC 3.31 σCC 2.80

−0.90 −0.64

−0.248 −0.28 −0.42

qC +0.70

fit to the critical temperature, saturated vapor pressures, and liquid densities. The parameters given in Table 2 are only a representative set of the many type 3-1 models possible with −0.26 |e| ≤qS ≤ −0.25 |e|. These parameters are close to one of the four H2S force fields proposed by Kamath et al.11 Model Type 4-1. To study the effect of changing the ratio of dipole to quadrupole moments of H2S, the negative charge can be shifted from S by a distance δS−X toward the H atoms along the H−S−H angle bisector. The type 4-1 model requires optimization of four parameters (εSS, σSS, qX, and δS−X). The van der Waals interactions continue to be modeled using a single LJ site on the S atom. Once again, the liquid-phase RP is found to be fairly insensitive to the LJ parameters (as long as the liquid density is well reproduced), but mostly determined by qX and δS−X. Multiple combinations of qX and δS−X yield an accurate liquid-phase RP. For each combination, the LJ parameters were obtained by fitting to the critical temperature and liquid density, and the four-parameter combination resulting in the correct slope of logarithmic pressure versus inverse temperature line is used to obtain the optimized set of parameters listed in Table 2. Model Type 3-3. Model type 3-3 differs from type 3-1 by additional LJ sites on the H atoms. Distributing the van der Waals interactions helps to effectively capture the shape of the H2S molecule. The 3-3 model requires optimization of five parameters: εSS, σSS, εHH, σHH, and qS. Given a reasonable set of LJ parameters, the partial charges can again be adjusted by fitting to the liquid-phase RP. Since the binary VLE for the CO2/H2S mixture is quite sensitive to σHH and εHH, this property dominates the choice of LJ parameters for the H atom. Finally, the LJ parameters for S are fitted to reproduce the single-component VLE of H2S. Multiple iterations of the two sets of LJ parameters are required to reach a complete set of parameters that can describe all the properties of interest. The optimized set of parameters for the proposed 3-3 model are listed in Table 2. This type is similar to the atom-based force field of Nath,10 but smaller values of the LJ well depth and diameter on the H atom could not yield a satisfactory normal boiling point for three-site models. Larger LJ parameters used in this work are mainly the result of fitting to the binary VLE with CO2. Although an extensive parametrization is performed, it is important to note that relatively large steps are used in obtaining the LJ parameters for hydrogen (ΔσHH = 0.5 Å and ΔεHH/kB = 10 K) because five parameters need to be optimized simultaneously and the properties involved in fitting are expensive to compute. Model Type 4-3. With six parameters (εSS, σSS, εHH, σHH, qX, and δS−X) to be optimized, the 4-3 model is the most complex

εHH/kB [K]

σHH [Å]

δS−X [Å]

0.1862 0.5 3.9 50 50 εC−H/kB 15.3 εOO/kB 79

0.98 2.5 2.5 σC−H 3.31 σOO 3.05

0.3

of the four models investigated in this work. Although the 4-1 model performs very well for pure-component properties, it shows large deviations for the H2S/CO2 binary mixture. Again, as also for the 3-3 model, distributing the LJ interactions over both S and H atoms allows one to adjust the strength of the H2S/CO2 interactions. In combination with the adjustments of the dipole moment/quadrupole moment ratio through the δS−X parameter and the corresponding partial charges, the 4-3 model can be parametrized to satisfy the accuracy requirements for all of the target properties. The parameters for this model are reported in Table 2. Methane and Carbon Dioxide. The explicit-hydrogen version of the TraPPE force field14 is used for the CH4 molecules. This is a rigid five-site model with LJ interaction sites representing the valence electrons placed on each C−H bond center and an additional LJ site with a very shallow well for the core electrons on the C atom. CO2 is modeled as a rigid, three-site molecule with LJ interaction sites and partial charges on each of the three atoms with parameters taken from the existing TraPPE force field.15 Force field parameters for both CH4 and CO2 are listed in Table 2.

■

SIMULATION DETAILS Gibbs Ensemble Monte Carlo (GEMC) simulations17−19 in the canonical (NVT) ensemble are used to simulate pure, binary (H2S/CH4 and H2S/CO2), and ternary H2S/CO2/CH4 VLE. For binary systems, one has the freedom to choose between the NpT or NVT ensemble due to the additional degree of freedom. However, due to its easier setup, the NVT ensemble is used in this work for the binary VLE simulations. A system size of 1000 total molecules is used for all GEMC simulations of unary, binary, and ternary systems. Following the standard for the TraPPE force field, LJ interactions are truncated at 14 Å and analytical tail corrections are applied (this rcut value is slightly larger than those used for the development of the CH4 and CO2 models). The Ewald summation method20 with a screening parameter of κ = 3.2/rcut and Kmax = κLbox + 1 for the upper bound of the reciprocal space summation is used for the calculation of the Coulomb energy. The total system volume is adjusted for each state point to yield a vapor phase containing about 20% of the molecules.21 Since this can lead to rather large volumes (and, correspondingly, linear dimensions) of the vapor-phase box, rcut for the vapor phase is set to approximately 40% of the box length to reduce the cost of the Ewald sum calculations. Four kinds of Monte Carlo moves, including translational, rotational, volume exchange, and particle transfer moves, are used to sample the phase space. The coupled− decoupled configurational-bias Monte Carlo algorithm22 is used 7043

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B

remainder divided equally between translations and rotations. Due to the heterogeneous nature of the condensed-phase box, LJ tail corrections are not applied but a larger cutoff at 19 Å is used together with the Ewald summation for the Coulomb interactions. After equilibration for at least 50 000 MCCs, data are collected from production periods ranging between 400 000 to 600 000 MCCs for the different simulation temperatures. For all systems investigated in this work, eight independent Monte Carlo simulations are carried out and the statistical uncertainties reported in the following sections are the standard errors of the mean calculated from these independent simulations. For the calculation of the self-diffusion coefficient, the liquid density is first obtained using Monte Carlo simulations in the NpT ensemble at T = 206.5 K and p = 1 atm. Using this density as input for a 1000-particle system, a molecular dynamics run (with the GROMACS 4.6.3 software29,30 and the SHAKE algorithm for the rigid-body constraints31) in the canonical ensemble of 1 ns in duration is used for equilibration. This is followed by a 10 ns trajectory in the microcanonical ensemble. The self-diffusion coefficient is computed from a linear fit to the mean square displacement versus time (using multiple time origins) to only the region from 0.5 to 3.5 ns of the trajectory to avoid contamination from the initial ballistic region and insufficient statistics toward the end of the simulation. The error is estimated by dividing the trajectory into two parts.

to enhance the acceptance rate for particle transfer moves. The probabilities for volume and transfer moves are set to yield approximately one accepted move of each type per Monte Carlo cycle (MCC),21 where an MCC consists of N = 1000 randomly selected moves. The remaining moves are divided equally between translations and rotations. An equilibration period of 50 000 to 100 000 MCCs is used for all simulations, and between 100 000 and 200 000 MCCs are used for the production phase. The static RP values of the liquid and solid phases are calculated from fluctuations of the system dipole moment sampled during a canonical-ensemble simulation:23 εD = 1 +

1 (⟨M2⟩ − ⟨M ⟩2 ) 3ε0kBVT

(3)

where M, ε0, and T are the total system dipole moment, the permittivity of free space, and the absolute temperature, respectively; V is the volume of the simulation box that is determined from a prior simulation. For the liquid phase, the equilibrium density for the model system is obtained by performing an isotropic NpT simulation24 for 500 molecules at T = 194.6 K and p = 1 atm. For the solid phase, the simulation is initiated with 500 molecules placed on Fm3̅m lattice sites with a 5 × 5 × 5 supercell, and an anisotropic (orthogonal) NpT simulation, 25 allowing the cell lengths to relax independently, but maintaining the cell angles orthogonal, is carried out at T = 157.9 K and p = 1 atm. After equilibration for 100 000 MCCs, the solid is found to maintain its face-centered cubic structure and the average density is obtained from a production period of 100 000 MCCs. The average densities obtained for the liquid and solid phases are then used for the subsequent canonical-ensemble simulations to determine the corresponding static RP via eq 3. These simulations consist of equilibration and production periods of 50 000 and 100 000 MCCs, respectively. For the isobaric−isothermal simulations, the probability of volume moves is set to result in approximately one accepted move per MCC and the remaining moves are divided equally between rotations and translations, as is also done for the canonical-ensemble simulations. In addition, anisotropic (orthogonal) NpT simulations are carried out to obtain the lattice parameter at T = 142 K and p = 1 atm. Isobaric−isothermal ensemble simulations are performed with 500 molecules at T = 298 K and p = 31 bar to obtain the intermolecular radial distribution functions at the same conditions as the experimental data.26 The system is equilibrated for 100 000 MCCs, followed by a production run of 200 000 MCCs. To determine the triple point, Gibbs ensemble Monte Carlo simulations17−19,27 are extended to cover solid−vapor equilibria and metastable states on the VLCC using a slab setup introduced by Chen et al.28 In this setup, the condensed-phase simulation box contains a slab of material surrounded by vapor and the other box contains a homogeneous vapor phase; the condensed-phase box is elongated along the z-axis and includes two surfaces parallel to the xy-plane. These simulations are started using a solid slab consisting of 2744 molecules (7 × 7 × 14 unit cells) and an empty vapor box. The vapor box volume is adjusted for each temperature to yield about 100 molecules in the vapor phase at equilibrium. For the condensed-phase box, volume changes are carried out by randomly changing one of the cell lengths while keeping the orthogonal shape. Again, the fraction of moves is adjusted to yield about one accepted volume and one accepted particle swap move per cycle with the

■

RESULTS AND DISCUSSION Liquid-Phase Relative Permittivity. The liquid-phase RP is a convenient measure of the polarity of a solvent. Due to induced polarization effects, the (average) charge distribution of a molecule in a condensed-phase environment differs from the distribution found for an isolated molecule in the gas phase. Thus, nonpolarizable models usually employ a larger dipole moment than present in the gas phase. It should be noted, however, that the liquid-phase RP not only depends on the molecular dipole moment but also reflects the degree of orientational ordering in the liquid. For weakly dipolar molecules, such as H2S, the computation of the liquid-phase RP is relatively inexpensive compared to the computation of binary vapor−liquid equilibria with a nonpolar molecule. Thus, evaluation of the liquid-phase RP is used here to quickly narrow the range of partial charges for a given type of H2S model. The static dipole moments for the four types of models optimized in this work are listed in Table 3. As can be seen, the μD values for three of the models are about 30% larger than the experimentally determined gas-phase value of 0.98 D.32 The exception is the 3-1 model for which compromises are necessary to achieve a reasonable VLCC (see below). The Table 3. Dipole Moment, Liquid-Phase Density, and Relative Permittivity at 194.6 K and 1 atma ρliq [g/mL]

εD

3-1 4-1 3-3 4-3

1.12 1.32 1.25 1.27

0.9681 0.9741 0.9751 0.9771

8.01 9.21 9.31 8.71

expt

0.9832

0.98133

8.998

model

μD [D]

a

The subscripts denote the statistical uncertainty of the last digit (standard error of the mean).

7044

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B liquid-phase RP values at T = 194.6 K and p = 1 atm (i.e., close to the triple point) for the other three models are found to deviate by less than 4% from the experimental value of 8.99,8 whereas the 3-1 model yields a value of only 8.0. Overall, it is clear that computation of the liquid-phase RP helps to constrain the range of suitable force field parameters, but it is not sufficient to discard specific model types. Unary Vapor−Liquid Equilibria. Vapor−liquid coexistence curves and Clausius−Clapeyron plots for all four models are shown in Figures 2 and 3 (and numerical data are provided

Table 4. Normal Boiling Point, Critical Properties, and Accentric Factor for Hydrogen Sufidea model

Tb [K]

Tc [K]

ρc [g/cm3]

pc [bar]

ω

3-1 4-1 3-3 4-3

206.91 213.32 212.62 212.62

378.85 377.06 375.75 374.56

0.3542 0.3602 0.3602 0.3612

953 934 943 933

0.0094 0.0814 0.0873 0.0923

expt

212.8537 212.8733

373.0740 373.3238 373.433 373.5441

0.34833 0.34941

89.4640 89.737 90.0538 90.0741

0.09633

a

The subscripts denote the statistical uncertainty of the last digit (standard error of the mean).

properties, saturated liquid and vapor densities are fitted to the scaling law for the critical temperature34 ρliq (T ) − ρvap (T ) = A(T − Tc)β

and the law of rectilinear diameters

(4)

35

ρliq (T ) + ρvap (T )

= ρc + B(T − Tc) (5) 2 where β = 0.326 is the universal critical exponent for threedimensional systems. The critical pressure is determined by extrapolation of the saturated vapor pressures to the computed critical temperature using the Antoine equation. Only VLCC data at T ≥ 340 K (i.e, T ≥ 0.9 Tc) are used for the determination of the critical point. The normal boiling point is obtained by interpolation of the saturated vapor pressures using the Clausius−Clapeyron equation for only the two data points closest to Tb. With the exception of the type 3-1 model, the other three models allow for a very accurate description of the unary VLE of H2S. The 3-1 model is highly simplified and, with the constraint of yielding an acceptable liquid-phase RP, is incapable of capturing the correct dependence of p on T (see Figures 3 and 4), and the best parametrization yields a 3-1 model that underestimates Tb but overestimates Tc and pc. As a result, the 3-1 model qualitatively fails for the accentric factor that is underpredicted by an order of magnitude. The 3-1 model also yields larger deviations for the orthobaric liquid density and significantly under predicts the enthalpy of vaporization at T ≤ 320 K (see Figure S1 and Table S3 in the Supporting Information). It should be noted here that the parameters of the 3-1 model found in this work are nearly the same as those for model A developed by Kamath et al.11 However, the normal boiling point and critical properties reported here are somewhat different, likely due to the use of a larger cutoff in this work (14 versus 10 Å) and the inclusion of only near-critical data for the determination of the critical point in this work (T ≥ 0.9 Tc versus T ≥ 0.7 Tc); the applicability of the Ising scaling law far away from Tc is questionable particularly for polar compounds. Based on the deficiencies of the 3-1 model for the computation of the unary VLE, this model is abandoned and not considered for the calculation of other properties. For the 4-1, 3-3, and 4-3 model types, parameters can be found that satisfy the target accuracy for the unary VLE of H2S. The normal boiling point33,37 is reproduced to within 0.2% (just outside of the statistical uncertainties) by all three models,

Figure 2. Vapor−liquid coexistence curves for hydrogen sulfide. Experimental data by Beaton et al.,36 Cubitt et al.,37 Goodwin,33 and Reamer et al.39 are represented by black plusses, diamonds, squares, and crosses. The orange up triangles, green down triangles, magenta left triangles, and cyan right triangles show the simulation results for the 3-1, 4-1, 3-3, and 4-3 models, respectively. Statistical uncertainties for the simulation data are smaller than the symbol size.

Figure 3. Saturated vapor pressure versus inverse temperature for hydrogen sulfide. Symbol styles and colors are the same as those in Figure 2. Statistical uncertainties for the simulation data are smaller than the symbol size.

in Tables S1−S2 in the Supporting Information), respectively, and numerical values of the normal boiling points, critical properties, and accentric factor for these models are summarized in Table 4. For the calculation of the critical 7045

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B

Sabocinski and Kurata43 covering a wider range of 225 K ≤ T ≤ 364 K. Recently, Chapoy et al.44 reported binary data for 258 K ≤ T ≤ 313 K. Considering the practical importance of this system and also for improving the transferability of the H2S force field, binary VLE for the H2S/CO2 mixture are considered as an additional criterion in the force field development. The critical temperatures of CO2 and H2S are 304 and 373 K, respectively. Thus, T = 293.16 and 333.16 K were selected for the calculation of the binary VLE to include data below and above the critical temperature of the low-boiling component (CO2) and because two experimental data sets are available at these temperatures. Figures 5 and 6 show the pressure−

Figure 4. Relative deviations in orthobaric liquid density and saturated vapor pressure from the recommended experimental values33 as a function of temperature. Symbol styles and colors are the same as those in Figure 2.

and the saturated vapor pressures for all three models fall within 5% of the experimental values over the entire stability range of the liquid phase (see Figure 4). All three models slightly underpredict the orthobaric liquid density for T < 300 K, whereas it is overpredicted in the near-critical region. The 43 model is slightly more accurate for ρliq over the full temperature range than the 4-1 and 3-3 models. The three models also slightly overpredict the enthalpy of vaporization (see Figure S1 and Table S3 in the Supporting Information). The critical temperature is overestimated by 1%, 0.7%, and 0.3% for models 4-1, 3-3, and 4-3, respectively. The critical pressures determined for these three models are very close to each other but fall about 4% above the average of the experimental values.37,38,40,41 Similarly, the critical densities for the three models are very close to each other but are about 3% higher than the experimental values.33,41 The 4-3 model yields a more accurate accentric factor than the 4-1 and 3-3 models. The ability of the 4-1, 3-3, and 4-3 types of models to reproduce the unary VLE and the liquid-phase RP of H2S indicates that either adjusting the dipole/quadrupole ratio or capturing the nonspherical shape of the H2S molecule provides sufficient flexibility for the model parametrization when only considering these properties. Thus, additional experimental data are needed to select the best model. Here, it is found that the binary VLE of H2S/CO2 can help to resolve this issue. Binary Mixture with Carbon Dioxide. Carbon dioxide and hydrogen sulfide constitute a large fraction of the waste streams resulting from oil and gas sweetening. These compounds can be reinjected into reservoirs either as compressed vapors or as liquified acidic streams through compression followed by cooling. These processes have aroused significant interest in the determination of the binary CO2/H2S VLE. Already in 1953, Bierlein and Kay42 reported on measurements of the binary VLE in the temperature range from 273 to 333 K. This was followed in 1959 by the work of

Figure 5. Pressure−composition diagram (top) and separation factor (bottom) for the H2S/CO2 mixture at T = 293.16 K. The experimental data of Bierlein and Kay42 and of Chapoy et al.44 (at T = 293.47 K) are shown as crosses and stars, respectively. The green up triangles, magenta left triangles, and cyan right triangles depict the computed data for the 4-1, 3-3, and 4-3 models, respectively.

composition diagrams and separation factors for the CO2/H2S mixture at 293.16 and 333.16 K, respectively. The pxy data for the 4-1 model are shifted significantly to the right (i.e., higher CO2 mole fractions for both phases at a given pressure) compared to the experimental data, whereas the data for the 3-3 and 4-3 models closely follow the experimental data. This behavior implies that the 4-1 model somewhat overpredicts the strength of the H2S/CO2 interactions. Actually, the data for the 4-1 model closely match the predictions from the Peng− Robinson equation of state45 with the kij parameter set to zero, whereas kij = 0.0960 is required to correlate the experimental 7046

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B

Figure 7. [H](H2S)−[O](CO2) radial distribution function and the corresponding number integrals at T = 293.16 K and xCO2 ≈ 0.45 for the 4-1, 3-3, and 4-3 models shown as green, magenta, and cyan lines, respectively.

where the liquid phase is close to equimolar. The RDF for the 4-1 model shows a shoulder for separations smaller than 3 Å, whereas the corresponding shoulder is shifted to larger separations by about 0.5 Å for the 3-3 and 4-3 models. As can be seen from the number integrals, these shoulders account for the nearest neighbor oxygen atom around a given hydrogen atom. This peak represents a weak hydrogen bond that is further weakened by the LJ site on the hydrogen donor atom for models 3-3 and 4-3. Based on the difference in the LJ diameter between the sulfur atom in the present models and the oxygen atom in CO2 (but similar values are common for water and alcohol models), one can expect that the S−H distance cutoff for a hydrogen bond should be larger by ∼0.3 Å than the commonly used value of 2.6 Å for O−H distances.46 The weakening of the hydrogen bond for the 3-3 and 4-3 models decreases the solubility of carbon dioxide in the H2S liquid phase and, hence, allows one to reproduce the positive deviations from regular mixing. As mentioned above, the Lorentz−Berthelot combining rules are used here consistently to compute the LJ parameters for unlike interactions belonging to the same type of molecule or to two different types of molecules. Other combining rules are available that yield a larger value for the unlike LJ diameter and a smaller value for the unlike well depth than the Lorentz− Berthelot combining rules when applied to sites with different LJ diameters.16 Use of these different combining rules for the S(H2S)−C(CO2), S(H2S)−O(CO2), H(H2S)−C(CO2), and H(H2S)−O(CO2) interactions would indeed significantly weaken the interactions between H2S and CO2 molecules: e.g., by more than 29% and 17% for the S(H2S)−C(CO2) and S(H2S)−O(CO2) well depth, respectively, with the 4-1 model using the Waldman−Hagler rules, i.e., much more than the ∼10% adjustment required for the Peng−Robinson equation of state.44 However, an even bigger effect would be observed for the S(H2S)−H(H2S) unlike interaction that would require significant reparameterization of the 3-3 and 4-3 H2S models (and also of the TraPPE CH4 model) and to a lesser extent of the TraPPE CO2 model. Such a complete reparameterization is beyond the scope of the present work. Based on its failure to predict the H2S/CO2 VLE with satisfactory accuracy, the 4-1 model is abandoned at this point. However, neither the binary VLE nor the structural data for this

Figure 6. Pressure−composition diagram (top) and separation factor (bottom) for the H2S/CO2 mixture at T = 333.16 K. The experimental data of Kay and Bierlein42 and of Sobocinski and Kurata43 are shown as crosses and stars, respectively. Symbols and colors for the simulation data are the same as those in Figure 5.

data at T = 293 K.44 Apparently, the 3-3 and 4-3 models are able to capture this positive deviation from regular mixing that requires a weakening of the unlike interactions. The deviations for the separation factor are significantly smaller with the 4-1 model underpredicting the separation factor at lower CO2 mole fraction (xCO2 < 0.5 at T = 293.16 K and xCO2 < 0.25 at T = 333.16 K) and overpredicting at higher mole fractions, whereas the 3-3 and 4-3 models perform better at low and high CO2 mole fractions (but of course not at the point where the data for the 4-1 model intersect the experimental data). In order to correctly predict the binary VLE for the H2S/ CO2 mixture with nonpolarizable models, the key ingredient that is missing in most models in the literature is accounting for the nonspherical repulsive shape of the triatomic H2S molecule which requires additional LJ sites on the hydrogen atoms. The H2S force field by Nath includes these additional LJ sites, but both the LJ diameter and well depth are too small (see Table 2) to yield a satisfactory representation of the unlike interactions.10 It is found here that relatively large values for the σHH and εHH parameters are essential to fit the binary VLE with CO2. Figure 7 illustrates the [H](H2S)−[O](CO2) radial distribution functions (RDF) and the corresponding number integrals (average number of [O] atoms surrounding a given [H] atom) for the 4-1, 3-3, and 4-3 models at a state point 7047

DOI: 10.1021/acs.jpcb.5b02536 J. Phys. Chem. B 2015, 119, 7041−7052

Article

The Journal of Physical Chemistry B

due to polarization in polar environments) and LJ contributions to the cohesive energy. Ternary Mixture with Methane and Carbon Dioxide. Ternary mixtures involving CH4, CO2, and H2S are of considerable importance in designing new processes for natural gas sweetening. To further assess the predictive capabilities of the TraPPE force field, the ternary phase diagram for the H2S/ CO2/CH4 mixture at T = 238.76 K and p = 20.68, 34.47, and 4 48.26 bar (also see SI) are computed. In this case, TCH