Molecular simulation of the solubility and diffusion of carbon dioxide and hydrogen sulfide in polyethylene melts

Available online at www.sciencedirect.com

Fluid Phase Equilibria 261 (2007) 168–175

Molecular simulation of the solubility and diffusion of carbon dioxide and hydrogen sulfide in polyethylene melts Franc¸ois Faure a,b , Bernard Rousseau a,∗ , V´eronique Lachet b , Philippe Ungerer b a

b

Laboratoire de Chimie Physique, Universit´e Paris Sud 11, CNRS, 91405 Orsay Cedex, France Institut Fran¸cais du P´etrole (IFP), 1&4 Avenue de Bois-Pr´eau, 92852 Rueil-Malmaison Cedex, France Received 30 April 2007; received in revised form 9 July 2007; accepted 9 July 2007 Available online 17 July 2007

Abstract We present molecular simulation results of solubility and self-diffusion of carbon dioxide and hydrogen sulfide in linear polyethylene melts. Simulations were carried out at 433 K and pressures in the range 0.1–10 MPa. Solubilities were computed using Monte Carlo simulations in the osmotic ensemble, with constant number of polymer chains, pressure, temperature and gases chemical potential. Diffusion coefficients were obtained from equilibrium molecular dynamics. It is found that hydrogen sulfide has a higher solubility in polyethylene than carbon dioxide at a given pressure. Self-diffusion coefficients are found to slightly increase with increasing penetrant weight fraction. © 2007 Elsevier B.V. All rights reserved. Keywords: Solubility; Henry constant; Diffusion; Polymer; Gas; Monte Carlo; Osmotic ensemble; Molecular dynamics

1. Introduction The permeation of gases in polymers is a major issue in the oil and gas industry. Flexible pipes convey gases from the sea bottom to the surface in the exploitation of offshore oil pools. Besides the metallic shields which provide mechanical resistance to the pipe, polymeric layers are inserted into the metallic structure to ensure the impermeability of the pipe, both from the sea water outside and from the fluids conveyed inside. In order to prevent gases from damaging the structure of the pipe, it is necessary to acquire a good understanding of the permeation phenomenon and to develop reliable and quantitative models to predict permeation behaviours in the polymeric layers. The permeability P of a penetrant through a material is given by P = S × D, where S is the solubility of the penetrant and D is the diffusion coefficient. The approach chosen in this work is molecular simulation, using molecular dynamics for the diffusion coefficient and Monte Carlo simulations for solubility. Since full equilibration of dense phases of long chains polymer melts is a difficult task to achieve in molecular simulations, special Monte Carlo moves were added to the code to improve the configurational space sampling. ∗

Corresponding author. E-mail address: [email protected] (B. Rousseau).

0378-3812/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2007.07.032

The classical way to study liquid–vapor systems is to use the Gibbs ensemble, introduced by Panagiotopoulos [1]. But this statistical ensemble is less appropriate in the case of gas mixtures of known composition. Moreover, due to their very low vapor pressure, the presence of polymer molecules in the gas phase is unrealistic. The transfer of such molecules from the polymer phase to the gas phase attempted in the Gibbs ensemble is therefore useless. Thus, an other statistical ensemble, the osmotic ensemble, has been put forward [2–4] to study solubilization of gases in polymers. In this one phase ensemble, the polymer-rich phase is constrained to be in equilibrium with a virtual gas phase. The equilibrium between the polymer-rich phase and the pure gas phase is ensured by insertion/deletion Monte Carlo moves where the chemical potential of the gas species is imposed. The number of polymer chains remains unchanged during the simulation. The polymer studied in this work is the linear polyethylene (PE), represented by chains containing 100 or 200 carbon atoms. Carbon dioxide and hydrogen sulfide gases were studied primarily, due to their relevance in the oil and gas industry, particularly in the context of the increase of natural gas production with a high acid gases content. In the next section, the osmotic ensemble and the simulation details are described. In the third section, we present the validation tests for the simulation methods and discuss several results

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175

on the solubility and diffusion of carbon dioxide and hydrogen sulfide at 433 K and different pressures.

Table 1 Molecular potential parameters Model

The Gibbs ensemble [1] is the reference ensemble to study liquid–vapor phase equilibrium. Each phase is simulated by a box and particles are transfered between the boxes to equalize chemical potentials. However, it does not allow to specify the chemical composition of a gas mixture solubilizing in a polymer matrix. The aim of this work is to simulate systems where the chemical composition, i.e. the fugacity, of each species is specified. We consider a system with p gas species in equilibrium with q polymer species (1 ≤ p ≤ p + q). The number of molecules of the q species (polymers) is constant whereas the number of molecules of the p species (gases) is variable. We note Ni is the number of molecules of the gas species i, V the volume of the simulation box, μi the chemical potential of the gas species i and U is the potential energy of the system. The configurational partition function Θ is: Θ=

 V

Ni

⎡ exp ⎣−βU +

p 

⎤ βμj Nj − βPV ⎦

(1)

j=1

U

where β = 1/kB T with kB the Boltzmann constant and T the temperature. This expression is consistent with those found in previous works [2,3]. The simulation of the polymer chains has been handled with the implementation of several Monte Carlo moves specially designed for long chain molecules, by Theodorou and coworkers: Concerted Rotation (ConRot) [5], Double Bridging (DB) and Internal Double Rebridging (IDR) [6–9]. Beside these three moves, reptation [10], regrowth [11] and flip have been used to sample the configurational space of the polymer chains. Classical translations and rotations were used for gas molecules as well as insertions/deletions. Since the pressure is imposed on the system, volume changes have also been used. 2.2. Infinite dilution method To evaluate the Henry constant of a penetrant, we used the infinite dilution method [12,3]. This method is only valid for low pressures, where the gas can be considered as ideal, i.e. there are no interactions between solubilized gas particles. Particle insertion tests of a phantom gas molecule are carried out in a pure polymer box. In the NpT ensemble, the weight-average Henry constant H W of the gas is obtained by the following expression:  HW =

MWPE MWgas



ρPE V NPE pT β VWtest NPE pT

Lennard–Jones

Electrostatic charge

˚ ) δ (A

/kB (K)

˚ ) σ (A

CH2 CH3

0.38 0.22

86.29 120.15

3.461 3.607

CO2 C O

0 0

H2 S S H Charge

0 – –

2. Simulation method 2.1. Monte Carlo

169

28.129 80.507 250 – –

– –

2.757 3.033

+0.652 −0.326

3.73 – –

+0.4 +0.25 −0.9

δ is the distance between the Lennard–Jones force center and the carbon atom in the anisotropic united atoms model.  and σ are, respectively, the energetic factor and the distance parameter of the Lennard–Jones potential. The C O ˚ long and the S H 1.340 A. ˚ The negative charge in H2 S is set bond is 1.149 A ˚ of the sulfur atom, along the symmetry axis of the molecule. The at 0.1862 A H–S–H angle value is 92◦ .

where MW is the molecular weight, ρPE the average density, NPE the number of polymer molecules and Wtest is the Rosenbluth weight for the insertion-test of the solute [13]. This method provides a way of evaluating solubilities on a range of pressure where Henry’s law is valid. The weight fraction w of the penetrant in the polymer at the pressure P is given by w = P/H W . 2.3. Molecular dynamics Molecular dynamics simulations were realized in the canonical ensemble (constant number of particles, volume and temperature) using a Nos´e-Hoover thermostat and an explicit time reversible integrator [14] with a timestep of 2 fs. The carbon–carbon bonds in polyethylene were constrained to their equilibrium value (see Table 1) using the RATTLE algorithm [15]. CO2 and H2 S were treated as rigid molecules. The same force field parameters, cutoff radius and long range corrections were used in molecular dynamics and Monte Carlo simulations (see below). Electrostatic interactions were computed using Ewald summation. Initial MD configurations were obtained from the output of the osmotic ensemble Monte Carlo simulations. For each penetrant type and each pressure, a MC configuration was selected with a number of penetrant molecules close to the average number computed from the MC run. The volume of the selected box was slightly modified to obtain a box at a density equal to the average density computed from MC. Initial velocities of centre of force were obtained from a Gaussian distribution and the configuration was then thermalized for 100 ps. The self-diffusion coefficient was computed from the meansquare displacement, msd(t), of the N penetrant molecules center of mass position, rcm i (t): N

(2)

msd(t) =

1  cm 2 [ri (t) − rcm i (0)] , N i=1

(3)

170

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175 Table 2 Intramolecular potentials

using the Einstein formulation [16]: D=

1 d lim msd(t) t→∞ 6 dt

(4)

For each system, the production run last for 10 ns and the msd(t) was computed over a time interval of 0.4 and 0.5 ns for CO2 and H2 S, respectively. During this period, the average mean square displacement was larger than tens of molecular diameter and we checked that the penetrant molecules reached the diffusive mode (i.e. the linear regime of msd(t) is observed). 2.4. Force field

k (cos(θ) − cos(θ0 ))2 2

(5)

where θ is the angle between two consecutive bonds, θ0 the equilibrium bending angle and k is the force constant. The torsion potential is described by Toxvaerd’s model [21]: Vtors (φ) =

8 

ai cosi (φ)

Bending CH3 –CH2 –CH2

CH2 –CH2 ˚ 1.535 A CH2 –CH2 –CH2 θ0 = 114◦ , k = 74900

Torsion CH3 –CH2 –CH2 –CH2 CH2 –CH2 –CH2 –CH2 a0 / . . . /a8 = 1001.36/2129.52/ −303.06/ −3612.27/2226.71/1965.93/ −4489.34/ −1736.22/2817.37 Constants k and ai for bending and torsion potentials are expressed in K.

The AUA4 Anisotropic United Atoms force field [17] has been used to model the polymer chains. The hydrogen sulfide molecule is described by Kristof and Liszi’s United Atom model [18], the carbon dioxide molecule by Harris and Young’s All Atoms model [19]. The details of these models are given in Table 1. Intermolecular interactions were represented by a Lennard–Jones potential, for dispersive–repulsive interactions, and by a coulombic potential, for electrostatic interactions. The molecular Lennard–Jones parameters are combined with ˚ has been used Kong’s mixing rules [20]. A cutoff radius of 10 A for Lennard–Jones interactions and it has been compensated with tail correction for long-range interactions. Since gas molecules are considered as rigid, intramolecular interactions only concern the polymer chains. Bending and torsion have been used for intramolecular interactions as well as Lennard–Jones potential for non-bonded intramolecular interactions. The bending potential is given by: Vbend (θ) =

Bond lengths CH3 –CH2

in the phase appears for any given pressure. Osmotic and Gibbs ensemble simulations show identical results. As far as these “single penetrant gas/polymer” systems are concerned, it seems that both methods are equivalent. The infinite dilution simulation reveals that the calculated weight-average Henry constant fits well with the beginning of the experimental curve. The Henry’s regime is valid in a rather small range of pressure ( 1 MPa) for this system. Results are in very good agreement with the experimental data. However, the global overestimate in the gas concentration is independant from the method used, especially at low pressures. 3.2. Solubilities at 433 K Longer chains have then been studied in order to be able to compare our results with experimental data on polyethylene. Two different systems have been studied. The first system is composed of 15 chains of C100 and the second is made of 8 chains of C200 . Simulations are held at 433 K, well above the melting point of linear polyethylene. Carbon dioxide has already been the target of experiments [23,24], but no data have been found for hydrogen sulfide, which is highly toxic. The concentration–pressure curve is shown in Fig. 2.

(6)

i=0

where φ is the torsional angle of three consecutive bonds. The trans configuration corresponds to the zero value of φ. The numerical values of the parameters for these potentials are given in Table 2. No stretching potential is used, bond lengths are constrained to their equilibrium value. 3. Results 3.1. Validation The solubility of CO2 in short PE chains was first studied to test the different Monte Carlo simulation methods. The box contained 40 chains of n-C44 at 423 K. The pressure–composition curve of the system is plotted in Fig. 1. Simulations are in good agreement with the experimental points, though a slight overestimate of the concentration of gas

Fig. 1. Pressure–composition curve of CO2 in n-C44 at 423 K. Filled symbols are experimental points [22]. Empty circles are osmotic simulations and empty triangles are Gibbs ensemble simulations. Statistical errors are the size of the symbols. The line corresponds to the infinite dilution simulation.

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175

Fig. 2. Concentration–pressure curve of CO2 and H2 S in linear PE at 433 K. Filled symbols are experimental data of CO2 solubility (diamonds are Sato’s [23], triangles are Chaudhary’s [24]). Empty circles are osmotic simulations of CO2 in C200 . Stars are osmotic simulations of CO2 in C100 Squares stand for osmotic ensemble simulations of H2 S in C200 . The lines are infinite dilution simulations in C100 ; dashed line: CO2 , dotted line: H2 S.

Concerning CO2 , we have reported solubilities computed in C100 and C200 using the osmotic ensemble along with results at infinite dilution in C100 . It is observed that the solubility increases linearly with increasing pressure for both chain lengths. For the pressure range considered here, the adsorption isotherm is caracteristic of Henry’s regime. The Henry constants evaluated for these systems are shown in Table 3. For osmotic simulations, the Henry constants have been evaluated from a linear regression of the data. A reasonable agreement for H W is obtained from the different data sets with uncertainties of the order of 15–20%. This result is typical of such simulations [3], for which a low convergence is observed, even with a large amount of Monte Carlo steps. Comparison of our simulation data with CO2 experiments on low and high density polyethylene demonstrates that our simulations slightly overestimate solubility in polyethylene. This trend was already observed for CO2 /n-C44 and could be attributed to the force field used in this work. We checked that the volumetric properties of pure PE with the AUA4 force field are correctly represented by our model. An illustration of such result is given in Fig. 3 where the specific volumes of linear C100 calculated at 0.1 MPa and different temperatures in the range 400–500 K are compared with experimental data of Zoller [25] and Dee [26]. Furthermore, a structural property obtained in our simulations, the end-to-end distance, was compared with a correlation

171

Fig. 3. Full circles: Specific volume of linear C100 calculated at 0.1 MPa and temperatures in the range 400–500 K. Simulations are compared with experimental data of Zoller [25](open square) and Dee [26](open triangle) on linear C71 (Zoller) and C78 (Dee) (dotted lines) and C150 (full lines).

˚ 2 mol g−1 . based on experimental data [27]: R2ee /M= 1.25 A For C100 , we obtained an average end-to-end distance equal ˚ roughly 3% below the value given by the correlato 40.5 A, tion. Therefore, we tend to believe that the slight discrepancies observed between experimental and calculated CO2 solubilities originate in the description of gas–polymer interactions and/or gas–gas interactions at high gas concentration. Concerning H2 S simulations, a good agreement is observed between osmotic and infinite dilution simulations. We have again assumed that we are in the Henry’s regime and reported the H W in Table 3. Clearly, H2 S solubility is higher than CO2 solubility. It confirms what has been observed in short linear alkanes (n-C15 and n-C16 ) [28,29]. We have reported in Fig. 4 the pressure–composition curve for H2 S in n-C15 obtained from experiments [28] and simulations [30]. The very good agreement between experiments and simulations is a strong validation of the force fields used in this work to describe the H2 S/PE mixture. A single experimental data point [29] of phase equilibrium for the CO2 /n-C16 system is also reported. As can be seen, for a

Table 3 Henry constants (H W (MPa)) for CO2 and H2 S in PE melts at 433 K System CO2 /C200

CO2 /C100

Osmotica

Osmotica

Infinite

97.2

82.8

86.4

H2 S/C100 dilutiona

Infinite 52.6

dilutiona

H2 S/C200 Osmotica 37.8

For osmotic simulations, the Henry constants have been evaluated from a linear regression of the data. a Methods.

Fig. 4. Pressure–composition curve of H2 S in n-C15 at 426.6 K. Circles are experimental data from [28], squares are simulation data in the Gibbs ensemble. We report, for comparison, a single experimental data point (triangle) [29] for CO2 in n-C16 .

172

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175

pressure of ∼2.5 MPa, the concentration of CO2 is half that of H2 S, indicating that hydrogen sulfide is indeed more soluble in linear alkanes than carbon dioxide. 3.3. Diffusion results at 433 K During the past two decades, diffusion processes in polymers have been investigated by molecular dynamics. The method has been successful in providing a detailed description of the microscopic processes involved in the permeation of penetrants in polymer. Contributions from Takeushi [31], M¨uller-Plathe [32,33] and Boyd and coll. [34,35] have given a microscopic picture of the mobility of small gases like hydrogen, oxygen, methane in various hydrocarbon polymers. At low temperature (below Tg ), the penetrant is trapped for long periods (several tens-hundred of ps) in a cage of surrounding polymer and makes unfrequent jumps. As temperature increases, the size of jumps increases and the diffusion process evolves from a solid-like behaviour to a liquid-like behaviour, with a broad spectrum of frequency jumps through the polymer. This change in diffusion mechanism seems to be accompanied by a change of the activation energy [36] but is not especially connected with the glass transition [37]. For the systems and the temperatures investigated in the study, the diffusion behaviour is typical of a liquid-like mechanism, as can be seen from Fig. 5. We present in Fig. 6 the self-diffusion data for CO2 and H2 S in polyethylene (8 chains of C200 ) as a function of the gas weight fraction. The total simulation run was splitted into 10 independent blocks. For each instant t, the standard deviation ES(t) of msd(t) over the 10 blocks was computed. From this, a “msd min(t) = msd(t) − ES(t)”, “msd max(t) = msd(t) + ES(t)”, and msd(t) were plotted and the slope of each was computed, giving the statistical errors on D. The data are also given in Table 4 for completeness. It is found that the self-diffusion of gases slightly increases with increasing gas pressure, or equivalently, with the gas weight fraction. It is known that diffusion coefficients are independent of the penetrant gas pressure or concentration for supercritical gases with very low critical temperatures (such as oxygen or methane) [38]. Subcritical gases are much more soluble in polymers and consequently, the above behaviour is observed only at low gas concentration in polymers. As the penetrant gas pressure is raised, the polymer is plasticized by the penetrant gas and the Table 4 Self-diffusion and swelling values for the studied systems for different pressure and corresponding concentrations Gas

P (MPa)

wgas (%)

Dgas (10−9 m2 s−1 )

V/V (%)

CO2

0.1 2.0 3.0 10.0

0.21 2.5 3.4 15.7

5.6 ± 0.5 6.89 ± 0.3 6.92 ± 0.2 10.04 ± 0.1

3.1 4.7 5.1 20.4

H2 S

0.1 1.0 2.0

4.5 ± 0.3 6.2 ± 0.5 7.9 ± 0.5

1.7 5.3 14.0

0.152 1.97 5.61

The reference density is the pure gas density at 433 K and 0.1 MPa. It is assumed that polyethylene density is almost constant in the considered pressure range.

Fig. 5. (a) Displacement of three different CO2 molecules in 8 polyethylene (C200 ) chains versus simulation time. A large spectrum of jump frequencies can be observed. (b) Average mean squared displacement computed on a time interval of 400 ps for a simulation run of 10 ns. Results are reported for three different pressures at 433 K. The dotted lines are linear function of time, indicating the diffusive regime. The diffusion coefficients were computed from the slope in the diffusing regime, i.e. for t > 300 ps.

diffusion coefficient may increase rapidly, in some cases, exponentially with increasing pressure. The critical coordinates of CO2 and H2 S are, respectively, Tc = 304 K, Pc = 7.38 MPa and Tc = 373 K, Pc = 9.01 MPa. Therefore, gases studied are in a supercritical state, although H2 S is close from its critical point. As can be seen in Fig. 6, diffusion coefficients of both gases increase with concentration with a more pronounced trend for H2 S, in agreement with expected trends. We have also indicated in Table 4 the evolution of the swelling ( V/V ) with pressure. These swelling values should be considered with caution as we assumed that the density of the pure polyethylene is almost constant with pressure in the range considered here. However, they clearly indicate that an important swelling occurs in H2 S samples at lower pressures than in CO2 . Accordingly, self-diffusion increases more rapidly for H2 S than for CO2 . We have included in Fig. 6 the experimental diffusion data from Sato et al. of CO2 in high density polyethylene. The diffusion coefficients were determined by measuring pressure versus time during gas dissolution and fit the results with a diffusion equation. The experimental data were obtained at a temperature

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175

Fig. 6. Self-diffusion coefficients for CO2 and H2 S in polyethylene at 433 K versus gas weight fraction. Triangles down are simulations for CO2 and triangles up simulations for H2 S. Circles are experimental data from Sato et al. [23]] at 453.2 K.

of 453.2 K. We are not aware of any experimental value for the activation energy of CO2 in polyethylene in this temperature range. If we consider an activation energy value in the range 35–45 kJ/mol (based on the activation energy for methane in polyethylene, see [39] and references therein), we can roughly estimate an increase of the diffusion coefficient when temperature is raised from 433 to 453 K by a factor of 1.5–1.7, i.e. of the order of magnitude of the ratio between experimental and simulation data. In agreement with Sato et al. [23], we observe an increase of diffusion with gas concentration. This behaviour is usually interpreted by a plasticization of the polymer and, at a microscopic level, by an increase of the segmental mobility of the polymer chains. In order to test this hypothesis, we have studied the local dynamical processes in polyethylene for different gases concentrations. Our analysis is based on the study of the Rouse modes associated with the atomic coordinates of the chain atoms ri (t). For a chain containing n carbon atoms, the mode Xk (t) is written:    n 1 kπ 1 ri (t) cos i+ ; k = 1, . . . , n − 1 Xk (t) = n n 2 i=1 (7) Xk (t) represents the motion of the chain segment including n/k bonds. The autocorrelation function of theses modes can be written:   t Xk (t) · Xk (0) = exp − (8) Ck (t) = τk X2k (0) where τk is the relaxation time associated with mode k. For a realistic polymer chain, whose behaviour is different from the Gaussian chain introduced in the Rouse model, the Rouse mode autocorrelation functions can better be described by a stretched exponential form:  βk

t (9) Ck (t) = exp − ∗ τk

173

Fig. 7. Evolution of the relaxation times τk associated with the Rouse mode k versus k for different pressures for CO2 in polyethylene chains with n = 200. Circles are simulations at 0.1 MPa, triangles at 3 MPa and diamonds at 10 MPa. The mode number k describes the relaxation of a chain segment comprising n/k carbon–carbon bonds. The statistical error for the slowest mode is about 25%. It decreases with increasing k.

where the relaxation times τk∗ and the stretching parameters βk are associated with the local mode k. The relaxation times τk are recomputed from τk∗ and k using the following relation: τk = τk∗ /βk Γ (1/βk ) where Γ (x) is the Euler function. We present in Fig. 7 the evolution of the relaxation times τk for CO2 in polyethylene at three different pressures. The βk values are not represented here but no significant difference between modes and pressure was observed, all βk lying between 0.72 and 0.78. From the Fig. 7, it can be seen that the relaxation time associated with a sub-chain of length n/k decreases when k increased, i.e. short sub-chains relax more rapidly than long sub-chains. A remarkable feature is that not all modes seem to be affected in the same way when pressure is raised. Indeed, the dynamical properties of sub-chains of length ≈ 30 (k = 6 and 7) seem to be more strongly modified suggesting that this characteristic subchain length is more influenced by the plasticizing effect. The fact that self-diffusion increases with pressure seems to indicate that the dynamics of CO2 molecules is coupled with the relaxation of these particular modes. We believe that this kind of approach can bring interesting and new informations on the local dynamics of penetrants in polymers and a more systematic study of this behaviour is scheduled. 4. Conclusions We have studied the solubility and diffusion properties of carbon dioxide and hydrogen sulfide in linear polyethylene using molecular simulation approaches. Simulation were done at a temperature of 433 K and pressures in the range 0.1–10 MPa. Solubility was computed using Monte Carlo simulations in the osmotic ensemble. Infinite dilution simulations were used to compute Henry constant and estimate the solubility in the limit of low pressures. For carbon dioxide, a reasonable agreement with experimental data from Sato [23] and Chaudhary [24] was obtained. As expected from results on short alkane chains in comparable conditions, H2 S solubility is higher than CO2 solu-

174

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175

bility in PE by a factor of ∼2 in the Henry’s regime. This first simulation result on the hydrogen sulfide/polyethylene mixture would benefit from experimental validation. Starting from equilibrated configurations obtained in osmotic ensemble simulations, diffusion coefficients of the penetrants were computed using classical equilibrium molecular dynamics. A slight dependence of the diffusion coefficient with concentration has been observed for both gases. In agreement with Sato [23] experiments, diffusion increases with increasing pressure and/or gas weight fraction. This trend is more pronounced for hydrogen sulfide. The local dynamics of the polymer chains was studied as a function of gas concentration using the Rouse mode formalism. It was observed that, with increasing gas weight fraction, the relaxation times of Rouse modes were not uniformly modified. Modes 6 and 7 (k − values), corresponding to subchains of 30 carbon atoms, relax more rapidly than others with increasing gas weight fraction. This suggests that the intermediate modes are more sensitive to a plasticization effect than short and long modes (respectively, long and short sub-chains). The increased mobility observed indicates a coupling between segmental motion on a length scale of 30 carbons and carbon dioxide mobility. From the following results on solubility and self-diffusion, it can be inferred that hydrogen sulfide permeability at pressures larger than 1 MPa is higher than carbon dioxide’s. Acknowledgements The authors would like to thank IFP and ANRT for financial support through a PhD grant (F.F.). Special thanks to Marie-H´el`ene Klopffer and Bruno Flaconn`eche (IFP) for fruitful discussions, and Jean-Marie Teuler (LCP) for support in code development. References [1] A. Panagiotopoulos, Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble, Mol. Phys. 61 (1987) 813–826. [2] F. Escobedo, Novel pseudoensembles for simulation of multicomponent phase equilibria, J. Chem. Phys. 108 (21) (1998) 8761–8772. [3] B. Banaszak, R. Faller, J. de Pablo, Simulation of the effects of chain architecture on the sorption of ethylene in polyethylene, J. Chem. Phys. 120 (23) (2004) 11304–11315. [4] T. Spyriouni, I. Economou, D. Theodorou, Phase equilibria of mixtures containing chain molecules predicted through a novel simulation scheme, Phys. Rev. Lett. 80 (20) (1998) 4466–4469. [5] L. Dodd, T. Boone, D. Theodorou, A concerted rotation algorithm for atomistic Monte Carlo simulation of polymer melts and glasses, Mol. Phys. 78 (4) (1993) 961–996. [6] V. Mavrantzas, T. Boone, E. Zervopoulou, D. Theodorou, End-bridging Monte Carlo: a fast algorithm for atomistic simulation of condensed phases of long polymer chains, Macromolecules 32 (1999) 5072–5096. [7] N. Karayiannis, A. Giannousaki, V. Mavrantzas, D. Theodorou, Atomistic Monte Carlo simulation of strictly monodisperse long polyethylene melts through a generalized chain bridging algorithm, J. Chem. Phys. 117 (11) (2002) 5465–5479. [8] N. Karayiannis, V. Mavrantzas, D. Theodorou, A novel Monte Carlo scheme for the rapid equilibration of atomistic model polymer systems of precisely defined molecular architecture, Phys. Rev. Lett. 88 (10) (2002) 105503-1–105503-4.

[9] D. Theodorou, Understanding and predicting structure–property relations in polymeric materials through molecular simulations, Mol. Phys. 102 (2) (2004) 147–166. [10] R. Boyd, An off-lattice constant-pressure simulation of polymethylene, Macromolecules 22 (5) (1989) 2477–2481. [11] J. de Pablo, M. Laso, U. Suter, Simulation of polyethylene above and below the melting point, J. Chem. Phys. 96 (1992) 2395–2403. [12] F. M¨uller-Plathe, Calculation of the free energy for gas absorption in amorphous polypropylene, Macromolecules 24 (1991) 6475–6479. [13] P. Ungerer, B. Tavitian, A. Boutin, Application of Molecular Simulation in Oil and Gas Technology—Monte Carlo Methods, Editions Technip, Paris, 2005. [14] G.J. Martyna, M.E. Tuckerman, D.J. Tobias, M.L. Klein, Explicit reversible integrators for extended systems dynamics, Mol. Phys. 87 (5) (1996) 1117–1157. [15] H.C. Andersen, Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations, J. Comp. Phys. 52 (1983) 24–34. [16] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science publications, 1989. [17] P. Ungerer, C. Beauvais, J. Delhommelle, A. Boutin, B. Rousseau, A. Fuchs, Optimization of the anisotropic united atoms intramolecular potential for n-alkanes, J. Chem. Phys. 112 (2000) 5499–5510. [18] T. Kristof, J. Liszi, Effective intermolecular potential for fluid hydrogen sulfide, J. Phys. Chem. B 101 (1997) 5480–5483. [19] J. Harris, K.H. Yung, Carbon dioxide’s liquid–vapor coexistence curve and critical properties as predicted by a simple molecular model, J. Phys. Chem. 99 (1995) 12021–12024. [20] C. Kong, Combining rules for intermolecular potential parameters. II: Rules for the Lennard–Jones (12–6) potential and the Morse potential, J. Chem. Phys. 59 (1973) 2464–2467. [21] S. Toxvaerd, Equation of state of alkanes II, J. Chem. Phys. 107 (1997) 5197–5204. [22] K. Gasem, R. Robinson Jr., Solubilities of carbon dioxide in heavy normal paraffins (C20 – C44 ) at pressures to 9.6 MPa and temperatures from 323 to 423 K, J. Chem. Eng. Data 30 (1) (1985) 53–56. [23] Y. Sato, K. Fujiwara, T. Takikawa, S. Sumarno, H. Takishima, Masuoka, Solubilities and diffusion coefficients of carbon dioxide and nitrogen in polypropylene, high-density polyethylene and polystyrene under high pressures and temperatures, Fluid Phase Equilib. 162 (1999) 261–276. [24] B. Chaudhary, A. Johns, Solubilities of nitrogen, isobutane and carbon dioxide in polyethylene, J. Cell. Plastics 34 (1998) 312–328. [25] P. Zoller, The pressure–volume-temperature properties of three wellcharacterized low-density polyethylenes, J. Appl. Polym. Sci. 23 (1979) 1051–1056. [26] G.T. Dee, T. Ougizawa, D.J. Walsh, The pressure–volume–temperature properties of polyethylene, poly(dimethylsiloxane), poly(ethylene gycol) and poly(propylene glycol) as a function of molecular weight, Polymer 33 (16) (1979) 3462–3469. [27] L.J. Fetters, D.J. Lohse, D. Richter, T.A. Witten, A. Zirkel, Connection between polymer molecular weight, density, chain dimensions, and melt viscoelastic properties, Macromolecules 27 (17) (1994) 4639– 4647. [28] S. Laugier, D. Richon, Vapor–liquid equilibria for hydrogen sulfide + hexane, + cyclohexane, + benzene, + pentadecane, and (hexane + pentadecane), J. Chem. Eng. Data 40 (1995) 153–159. [29] B. Breman, A. Beenackers, E. Rietjens, R. Stege, Gas–Liquid solubilities of carbon monoxide, carbon dioxide, hydrogen, water, 1-alcohols and n-paraffins in hexadecane, octacosane, 1-hexadecanol, phenanthrene, and tetraethylene glycol at pressures up to 5.5 MPa and temperatures from 293 K to 553 K, J. Chem. Eng. Data 39 (4) (1994) 647–666. [30] P. Ungerer, V. Lachet, B. Tavitian, Application of molecular simulation in oil and gas production and processing, Oil & Gas Science and Technology 61 (3) (2006) 387–403. [31] H. Takeuchi, K. Okazaki, Molecular dynamics simulation of diffusion of simple gas molecules in a short chain polymer, J. Chem. Phys. 92 (9) (1990) 5643–5652. [32] F. M¨uller-Plathe, Diffusion of penetrants in amorphous polymers: a molecular dynamics study, J. Chem. Phys. 94 (4) (1991) 3192–3199.

F. Faure et al. / Fluid Phase Equilibria 261 (2007) 168–175 [33] F. M¨uller-Plathe, Molecular dynamics simulation of gas transport in amorphous polypropylene, J. Chem. Phys. 96 (4) (1992) 3200–3205. [34] P.V.K. Pant, R.H. Boyd, Molecular dynamics simulation of diffusion of small penetrants in polymers, Macromolecules 26 (4) (1993) 679–686. [35] R.H. Gee, R.H. Boyd, Small penetrant diffusion in polybutadiene: A molecular dynamics simulation study, Polymer 36 (7) (1995) 1435–1440. [36] M. Meunier, Diffusion coefficients of small gas molecules in amorphous cis-1,4-polybutadiene estimated by molecular dynamics simulations, J. Chem. Phys. 123 (2005) 134906–134907.

175

[37] J. Han, R.H. Boyd, Small-molecule penetrant diffusion in hydrocarbon polymers as studied by molecular dynamics simulation, Macromolecules 27 (19) (1994) 5365–5370. [38] S.A. Stern, B. Krishnakumar, S.M. Nadakatti, Permeability of polymers to gases and vapors, in: J.E. Mark (Ed.), Physical Properties of Polymers Handbook, AIP, Woodbury NY, 1996 (Chapter 50). [39] P.V.K. Pant, R.H. Boyd, Simulation of diffusion of small-molecule penetrants in polymers, Macromolecules 25 (1) (1992) 494–495.