Molecular Simulation of the Effect of Temperature ... - ACS Publications

5646

J. Phys. Chem. B 2008, 112, 5646-5660

Molecular Simulation of the Effect of Temperature and Architecture on Polyethylene Barrier Properties† Patricia Gestoso*,‡,£ and Nikos Ch. Karayiannis§,¶,⊥ Accelrys Ltd., 334 Cambridge Science Park, Cambridge, CB4 OWN, U.K., Rhodia Recherches, Centre de Recherches de Lyon, BP 62, 69192, Saint-Fons Cedex, France, Institute of Chemical Engineering and High Temperature Chemical Processes (ICE/HT-FORTH), 26500 Patras, Greece, and Department of Chemical Engineering, UniVersity of Patras, 26500 Patras, Greece ReceiVed: May 14, 2007; In Final Form: January 29, 2008

We present a multiscale approach for calculating the low-concentration solubility, diffusivity, and selectivity of small molecules through polymer matrixes. The proposed modeling scheme consists of two main stages; first, thoroughly equilibrated and representative poly(ethylene) (PE) atomistic melt configurations were obtained through the application of a Monte Carlo (MC) scheme based on advanced chain-connectivity altering moves (linear architectures) or the combination of localized MC moves followed by molecular dynamics. In the second phase, transition-state theory (TST), as proposed by Gusev and Suter [Gusev, A. A.; Suter, U. W. J. Chem. Phys. 1993, 99, 2228], was invoked in a coarser level of description to calculate the barrier properties of the studied macromolecules to small gas molecules at infinite dilution. The multiscale methodology was successfully applied on PE melts characterized by various molecular weights (MW) (from C78 up to C1000) and polydispersity indices at a wide range of temperature conditions. The effect of molecular architecture on the barrier properties was examined through the comparison between linear and short-chain branched structures bearing the same total number of carbon atoms. Simulation results were found to be in very good agreement with available experimental data. Additionally, the new scheme has been further validated by comparing the qualitative behavior of solubility, diffusivity, and selectivity with previously reported trends in the literature based on both experimental and simulation studies. The present study concludes that density plays a dominant role that determines the behavior of the polymer as a barrier material, especially in terms of diffusivity. Additionally, it is evidenced that short-chain branching has a small effect on the barrier properties of PE when the comparison is performed on purely amorphous samples. The hierarchical method presented here not only is faster when compared against conventional molecular dynamics simulations, but in some cases, like the vicinity of the glass transition temperature or for long polymer chain melts, it opens the way to the calculation of the barrier properties at realistic simulation times.

I. Introduction Permeation of small molecules through polymeric materials is of technological importance in a variety of areas. Specific industrial and commercial applications demand barrier properties of the end product that span a very broad range of characteristics; a polymer designed for use in food packaging has to be practically impermeable to certain gases such as oxygen (O2), whereas fast diffusion is desirable in the case of purification from residual monomers in polymer synthesis. Other examples constitute the gas separation in biosensors where the goal is to develop a barrier selective to one molecular species or the drug delivery systems in which the speed of the drug release from a polymer matrix has to be fully controlled.1 † Originally submitted for the “Keith E. Gubbins Festschrift”, published as the November 1, 2007, issue of J. Phys. Chem. C (Vol. 111, No. 43). * To whom correspondence should be addressed. Address: Accelrys Ltd., 334 Cambridge Science Park, Cambridge CB4 OWN, U.K. E-mail: [email protected]. ‡ Accelrys Ltd. £ Rhodia Recherches. § Institute of Chemical Engineering and High Temperature Chemical Processes. ¶ University of Patras. ⊥ E-mail: [email protected].

The increasing demand for tailored characteristics and the wide range of applications of the barrier materials requires the execution of a long series of exploratory experiments, focused on the design of superior materials,2 that may proven, in practice, to be either too expensive or time-consuming. In this context, novel computational methods may assort experimental research by providing reliable predictions of the sorption and diffusion of small molecules in dense polymer matrixes at ideal conditions.3 Nowadays, the continuous increase in computing power and the affordable budget in obtaining farms of workstations reduces the required computational time and lowers the cost of the modeling investment. Even more, computer simulations, by resorting to the atomistic level of description, may reveal relaxation mechanisms that play a prominent role in the barrier properties that otherwise could not be traced by conventional permeability experiments.4,5 From the theoretical point of view, sorption and diffusion in and through macromolecular systems are explained and predicted mainly through the concepts of free volume (FV),6-8 the dual-mode transport model,9-12 and molecular theories.13-16 The packing of chains and their constituent atoms, even in the case of very dense polymers, leaves holes (or microvoids),17 the sum of which forms the free volume of the system influenced by the polymer matrix and the applied conditions. In principle, these

10.1021/jp073676q CCC: $40.75 © 2008 American Chemical Society Published on Web 04/12/2008

Polyethylene Barrier Properties clusters of free space can be occupied by penetrant molecules, depending on the size and shape ratio between the molecule and the cavity; in other words, the penetrant must fit somehow in the free space. Additionally, the penetrant-polymer interaction in relation to the penetrant-penetrant and polymerpolymer interactions also plays an important role.3 The collection of all free cavities that can host the penetrants comprises the accessible volume (AV) of the system, which depends on the polymer matrix and the penetrant molecule. In the melt state, as chains experience increased mobility, the free volume is continuously redistributed in the system, and the penetrant molecules diffuse through its fluctuations. According to this concept, transport proceeds by the following mechanism: through the segmental motion of the polymer matrix, a void opens up in the close proximity of the penetrant molecule. The penetrant, drifted by the diffusive motion of the ambient polymer atoms, performs a quick move from its original position into a new void, leaving a cavity of similar size behind, which is filled by other polymer segments. On the basis of this qualitative picture, in a semicrystalline material, it is expected that solubility and diffusivity are much higher in the amorphous than through the crystalline phase since amorphous regimes are less dense and their free volume is redistributed more easily and frequently. Quantitative methods, theoretical approaches, and semiempirical relationships, based mainly on the concept of free volume, have been developed in the past years for the barrier properties of molten or glassy polymers. These relationships have tried to correlate bulk properties and molecular characteristics (such as volume and thermal expansion coefficient) to permeation.18 While they have been very successful in relating permeation to the aforementioned quantities, they cannot provide a detailed picture of the diffusion mechanism or fail to incorporate additional mechanisms that may have a profound effect on the barrier properties. As mentioned above, besides theory and experiments, computational modeling, a rather new but very promising approach, can assist in the study of the barrier properties of macromolecular systems by combining and integrating various simulation techniques applied on different time and length scales. In the present work, we demonstrate a new hierarchical and multiscale approach to calculate the infinite dilution solubility, diffusivity, and selectivity of small gas molecules (oxygen (O2), nitrogen (N2), argon (Ar), carbon dioxide (CO2), and methane (CH4)) in polyethylene (PE) matrixes characterized by various chain lengths, polydispersity indices, and molecular architectures (linear and short-chain branched) at different temperatures. The proposed methodology consists of two main stages. Initially, fast and robust equilibration of the PE systems of a wide range of length scales is achieved through extensive application of advanced chain-connectivity altering MC moves based on the end-bridging19,20 and double-bridging1-23 algorithms. In the second phase, transition-state theory (TST)24-27 is applied on a large ensemble of representative, fully relaxed, and statistically uncorrelated PE configurations to extract the barrier properties of the simulated macromolecules to O2, N2, Ar, CO2, and CH4 molecules at infinite dilution. The paper is organized as follows. Section II presents a concise literature survey of simulation approaches for the prediction of gas solubility and diffusivity in amorphous (molten or glassy) polymers and of past simulation efforts on PE and related materials. Section III describes the multiscale methodology followed to generate well-equilibrated structures of linear and short-chain branched PE melts and the approach implemented to convert the united-atom (UA) representation to an

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5647 explicit-atom (EA) description of the PE structures. The presence of hydrogen atoms is crucial for calculating the barrier properties as the UA representation overestimates diffusion coefficients.27 The barrier properties were calculated using the gsnet and gsdiff programs of the commercial package InsightII 400P+.28 Results concerning the barrier properties of polyethylene as a function of molecular weight (MW), degree of polydispersity, and molecular architecture (linear vs short-chain branched structures) at various temperatures are presented in section IV. Finally, section V summarizes our findings and conclusions and discusses future plans and extensions. II. Modeling of Gas Permeability in Polymers State-of-the-art simulation studies of the barrier properties of macromolecules use either molecular dynamics (MD)29 or transition-state theory (TST)24-27 methods in order to understand and, to a certain extent, predict the transport behavior of small gas molecules in polymers. Both approaches are complementary, with MD being the less coarse-grained. The strength of this method lies in the relatively few approximations made, mainly the assumptions made in the applied force field. Its major drawback is the computational cost, which limits its applicability to deal with penetrant/polymer systems exhibiting high diffusion coefficients at temperatures well above the glass transition temperature (Tg). Even in the melt state, small penetrants in complex macromolecular environments may exhibit anomalous diffusive behavior25,26,30 for time scales on the order of dozens of nanoseconds before Fickian diffusion settles in, rendering conventional MD practically inefficient to capture the true transport phenomena. In the most common implementation, MD simulations track the dynamics of both the polymer matrix and the penetrant molecules, which either coexist from the very beginning (i.e., the polymer chain is originally grown and relaxed in the presence of the gas molecules) or have the penetrants inserted a posteriori in the simulation cell in an energetically biased way to avoid overlaps with the existing polymer segments.30-47 The TST method, as introduced by Arrizi, Gusev, and Suter for glassy polymers,24-27 is well suited to extend the time scale of the observation as it does not require classical dynamics but involves a number of assumptions that must be justified. This approach is especially effective for the study of gas transport evolving via a succession of infrequent events where diffusion is too slow to be tracked by conventional MD.4,48-61 Additionally, to this date, the Gusev-Suter implementation of the TST method is restricted to polymer-penentrant systems where the matrix response to the guest molecule is elastic. This is a valid assumption for light gas molecules at infinite dilution but breaks down in systems where the penetrants induce massive deformation in the polymer matrix, as with ions in complex polymers or in the case of macromolecules being swollen by penetrants at high concentration. In the majority of the simulated polymerpenetrant systems, TST was found to provide very good to excellent predictions of the solubility and diffusivity (and thus permeability) of small spherical gas molecules (H2, O2, N2, and CH4).4,46-61 As the penetrant size increases and its shape becomes anisotropic (like the CO2 molecule), conventional TST fails to capture the corresponding transport behavior. Up to a certain degree, a solution to the aforementioned problems can be provided by the implementation and application of multidimensional TST as proposed by Greenfield and Theodorou.62-66 Gusev-Suter TST, as implemented in InsightII,28 is initiated as a 3-D fine-resolution grid is laid on the fully equilibrated polymer configuration covering the entire simulation box. Next,

5648 J. Phys. Chem. B, Vol. 112, No. 18, 2008

Gestoso and Karayiannis

a spherical probe, with a radius identical to that of the gas penetrant, is inserted in all grid points (x, y, z), and the resulting nonbonded energy, Eins(x, y, z), is calculated between the test probe and all atoms of the polymer matrix. Then, the amorphous cell is divided in low-energy regions (clusters of accessible volume containing local minima of the energy) and high-energy domains (excluded volume characterized by high intermolecular energy). Under the condition of low concentration (i.e., there are no potential interactions between solved molecules at infinite dilution) according to the Widom method,67,68 the excess chemical potential, µex, is given by

〈 ( )〉

Eins µex ) RT ln exp kBT

(1)

where R is the universal gas constant, T is the temperature, kB the Boltzmann constant, and the brackets, 〈〉, denote the average over all grid points and probe insertions. The low-concentration solubility coefficient, S, is related to the excess chemical potential, µex, through

( ) µex RT

S ) exp -

(2)

The identified sorption sites (microstates) of local energy minima are separated by high-energy barrier surfaces. Penetrant diffusion can be envisioned as a series of infrequent transitions between adjacent microstates controlled by first-order rate constants given by24-27

kifj )

kBT Qij h Qi

(3)

where kifj denotes the rate constant for the transition from microstate i to microstate j, h is Plank’s constant, and Qi and Qij are the partition functions of the penetrant in microstate i and on the separating surface between microstates i and j, respectively. The collection of all microstates form a 3-D coarsegrained spatial network on which penetrant molecules hop between adjacent microstates with rate constants defined in eq 3. The process of displacing penetrant molecules in the coarsegrained lattice is repeated over a large number of time steps and a large population of ghost walkers through a kinetic MC (kMC) scheme.69,70 The mean square displacement (MSD), 〈[rp(t) - rp(0)]2〉, is calculated as a function of time, t, from the trajectories of all penetrant walkers. At long time scales when normal (Fickian) diffusion is established, the low-concentration diffusivity coefficient, D, can be readily calculated from

{

}

〈[rp(t) - rp(0)]2〉 D ) lim tf∞ 6t

(4)

where the brackets, 〈〉, denote the average over all trajectories. Finally, the infinite-dilution permeability coefficient, P, is calculated simply as the product P ) DS. In the original formulation of TST for glassy systems, the polymer matrix was considered as frozen, a simplification which led to significant underestimation of the solubility and diffusivity coefficients. In the improved approach used in this work, thermal fluctuations of the constituent atoms are taken into account under the simplification that the matrix segments execute independent “elastic” motions, tethered to their equilibrium positions through harmonic springs. The displacements around the equilibrium positions follow a Gaussian probability density function

(

W(∆r1,∆r2,...,∆rN) ∝ exp -

N

∑ i)1

3∆ri2 2〈∆ri2〉

)

(5)

where N is the total number of polymer atoms, ∆ri ) ri - 〈ri〉 is the displacement from the equilibrium position, and 〈ri〉 and 〈∆ri2〉 is the corresponding mean square displacement, usually assigned an averaged value of 3∆2, with ∆2 being the “smearing factor” for all atoms, regardless of the type. By incorporating this rather simplistic form of the smearing factor, TST is limited to track the transport and sorption behavior of small molecules (H2, O2, N2) since their presence in the matrix does not affect the polymer environment. Larger molecules (CO2 or n-alkanes) may force the surrounding macromolecular segments to undergo significant rearrangements to accommodate their presence. In addition, complex segmental relaxation mechanisms like phenyl ring vibrations and flips that may have a profound effect on the “dynamic flexibility”4,71 of the system are not fully captured by the smearing factor of eq 5. Whether one chooses to adopt either MD or TST approaches, a key factor for the successful prediction of solubility and diffusivity coefficients is the generation and use of a large number of fully equilibrated, statistically uncorrelated, and realistic polymer configurations. In order to fulfill these requirements, it is necessary to apply accurate, finely tuned potential force fields and, in parallel, to implement modeling techniques that achieve robust and fast equilibration at all length scales. The equilibration of long, entangled polymer melts is a longstanding problem in the simulation community. The large spectrum of characteristic time and length scales renders most conventional simulation methods like MD practically inefficient to provide, within reasonable computational time, robust and rapid relaxation of the studied macromolecular systems. Even the state-of-the-art parallel implementation of MD is not capable of accessing simulation times exceeding a couple of microseconds, lying well below the longest relaxation time of a typical medium-MW polyethylene (PE) system at standard conditions. Recently, the introduction of advanced chain-connectivity altering Monte Carlo (MC) moves made possible the full-scale equilibration of truly long, entangled melt configurations characterized by various molecular weight distributions. Through the end-bridging (EB)19,20 and double-bridging (DB)21,22 MC moves, as well as their intramolecular analogues, polydisperse and monodisperse PE samples can be fully equilibrated within modest computational time. Very recently, the core EB and DB algorithms have been further modified and expanded to handle complex, nonlinear molecular architectures bearing short or long branches.23 Nowadays, the generation of vast MC-generated trajectories, consisting of many uncorrelated and realistic configurations, constitutes the imperative groundwork for hierarchical schemes exploring and predicting the static,21,22,72 dynamic, and topological73 properties of the simulated macromolecular systems in atomistic detail. In the present work, a multiscale modeling approach, based on chain-connectivity altering algorithms, is applied for the prediction of the barrier properties (solubility and diffusivity) of PE melts, characterized by either a linear or short-chain branched architecture, to small gas penetrant molecules (O2, N2, Ar, CO2, and CH4) at low concentration. The proposed methodology, being reliable and robust, outperforms conventional MD simulations in the long-range equilibration of dense polymeric systems by many orders of magnitude, especially at temperatures in the vicinity of the glass transition. While the atomistic chain-connectivity altering moves have been success-

Polyethylene Barrier Properties

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5649

TABLE 1: Details of the Simulated Linear PE Systems number of chains

number of atoms

polydispersity index

C78

40

3120

1.083

C142 C142C0 C500 C1000

22 22 6 3

3124 3124 3000 3000

1.083 1.000 1.053 1.053

TABLE 2: Details of the Simulated Short-Chain Branched PE Systems number of number number branches polydispersity temperature index (K) of chains of atoms per chain

temperature (K) 250, 300, 350, 400, 450, 500, 550, 600 300, 450 300, 450 300, 450 250, 300, 350, 400, 450, 500, 550, 600

fully applied on macromolecular systems with a chemical constitution more complex than that of PE (e.g., PI, PIB, PP, PEO, and PEG),74-78 the major disadvantage of the current multiscale modeling approach is its inability to treat even more complex macromolecular entities, for example, chains bearing phenyl rings either in the backbone or as substitutes or monomer units consisting of many different atoms. The relaxation of chemically complex macromolecular systems can be achieved, to a certain degree, by resorting to other methods based on static structure minimization and successive steps of short MD runs consisting of cooling-annealing compressing-decompressing,48-58 requiring a much larger amount of computational time and resources and producing a rather limited sample of relaxed configurations per cycle.4 Alternatively, one could transform the atomistic structure to a coarser representation through pretabulated potentials,79 then apply the chain-connectivity altering moves19-22 on the coarsegrained medium, and finally, return to the atomistic level of description using fine-graining techniques. III. Multiscale Modeling Approach The proposed hierarchical methodology consists of three main steps: (1) MC equilibration of linear and short-branched polyethylene structures at various temperatures is achieved through chain-connectivity altering algorithms and molecular dynamics. (2) Conversion of selected, representative, and fully equilibrated united-atom (UA) PE MC-generated frames to InsightIIcompatible explicit-atom (EA) configurations by adjusting the corresponding hydrogen atoms. Because the addition of hydrogen atoms may create energetic overlaps (and consequently unrealistic local packing), extra segmental relaxation is performed through a succession of short minimization/NVT MD simulations. (3) TST is invoked to calculate the barrier properties of the PE matrixes to small penetrant molecules. III.1. Systems StudiedsMolecular Model. All linear, polydisperse, and monodisperse PE systems simulated in the course of the present work are listed in Table 1. As in all previous works for polydisperse PE melt configurations, the equilibrium chain length distribution is fully controlled through the use of an appropriately selected spectrum of relative chemical potentials µ* cast in a semigrand canonical statistical ensemble.19-23 A uniform distribution of chain lengths was adopted in all polydisperse cases in the closed interval [Nav (1 - ∆), Nav (1 + ∆)], where Nav denotes the number-average chain length and ∆ is the half-width of the chain length distribution reduced by Nav. The polydispersity index, I, is related to ∆ through I ) 1 + ∆2/3. The MC scheme, used for the equilibration of all polydisperse linear systems, consists of the following moves: flip (or internal libration) (5%),80 rotation (3%), intermolecular reptation (15%),19,23,81,82 concerted rotation (20%),83 end-bridging (15%),19,20 double-bridging (20%),21,22 intramolecular double-rebridging

C132C1 C126C2 C112C5 C110C8

22 22 22 22

3124 3124 3124 3124

10 8 6 4

1.000 1.000 1.000 1.000

300, 450 300, 450 300, 450 300, 450

(10%),21,22 intramolecular end-bridging (10%),82 and volume fluctuation (2%). In the case of strictly monodisperse systems, the condition of constant molecular weight throughout the whole MC simulation prohibits the use of algorithms that, by altering the lengths of the participating chains, impose even a small degree of polydispersity. As a consequence, the corresponding MC scheme for the monodisperse C142C0 configurations is reduced to the following mix: flip (5%), rotation (3%), reptation (15%), concerted rotation (20%), double-bridging (35%), intramolecular double-rebridging (10%), intramolecular endbridging (10%), and volume fluctuation (2%). Short-chain branched systems require the use of a wider sample of advanced local or chain-connectivity altering moves due to the presence of multiple branch points. Consequently, a number of special moves (like the double concerted rotation and the H-BR (H-shaped branch rebridging))23 was developed and implemented recently, allowing the efficient equilibration of the branch point and its adjacent atoms as well as the longrange conformational characteristics of long-chain branched structures (H-shaped).23 Details about the simulated shortbranched configurations can be found in Table 2. The notation used throughout the whole manuscript to describe these systems is CxCy, where x and y stand for the total number of carbon atoms in the backbone and in each branch, respectively. To render the comparison of the barrier properties as a function of the molecular architectures fair, all combinations of branch lengths and frequencies were selected so that the MW of the branched chains corresponded to that of a linear C142C0 chain (i.e., 1990 g/mol). In order to keep the branch length and frequency fixed (to avoid the introduction of additional parameters that may affect the barrier properties), only localized MC moves were adopted for the simulation of short-chain branched systems, significantly reducing the efficiency of the MC protocol in equilibrating the long-range characteristics of the macromolecules. For this reason, selected MC-generated configurations were subjected to an additional relaxation provided through NPT MD simulations on the order of 100 and 400 ns at T ) 450 and 300 K, respectively. All NPT MD simulations were conducted using the LAMMPS (large-scale atomic/molecular massively parallel simulator) software.84,85 For the integration of the equations of motion, the multiple time step algorithm rRESPA86,87 was employed with small and large integration time steps set to 1 and 5 fs, respectively. In all NPT MD simulations of shortchain branched systems, the Nose´-Hoover thermostat88,89 and the Andersen barostat90 were employed. The molecular model employed in the MC simulations for the initial phase of the full-scale equilibration of the samples is identical to that used in our previous studies of linear and nonlinear PE melts, which provided excellent estimates of the volumetric,21-23 conformational,21-23 and dynamic91 properties for a variety of temperature conditions. This scheme adopts the united-atom representation, according to which all CH, CH2, and CH3 units are treated as single, spherically symmetric interacting sites. More details about the employed potential force field can be found elsewhere.91



Figure 1. Density dependence on the average molecular length as obtained from MC simulations on polydisperse PE systems at T ) 300 and 450 K (P ) 1 atm). Experimental data at the corresponding temperatures and best fits of hyperbolic functions on the simulation predictions are shown for comparison purposes.

III.2. Monte Carlo Equilibration: Efficiency and Volumetric Properties. The decay of the autocorrelation function of the end-to-end unit vector, , as a function of CPU time was used to verify the equilibration of the melts, as reported previously for similar systems.20-22 The generation of MC trajectories containing vast amounts of independent and fully equilibrated conformations allows the calculation of the density, F, and the mean-squared end-to-end distance, , of the simulated PE systems as average values over all recorded frames after the initial equilibration period. On the basis of the calculated values of F and , representative configurations from the trajectory with instantaneous values of F and that match exactly the average properties can be selected. This procedure is repeated over all linear systems and simulation conditions. As the temperature decreases, the number of chosen frames increases in an effort to reduce the statistical error in the successive calculation of the barrier properties by sampling over an extended set of different configurations. In the case of the short-chain branched analogues, the corresponding frames are selected after the application of the final NPT MD simulation. Figure 1 presents the variation of density, F, with the average molecular length, Nav, as obtained from MC simulations on polydisperse samples at T ) 450 and 300 K (P ) 1atm). It can be seen that, regardless of the applied temperature, the density reaches a plateau value at around N ∼ C224; no appreciable effect on the density is observed for longer chain lengths. This is indicative of true polymeric behavior, and it is in agreement with experimental qualitative trends. Also shown in Figure 1 are experimental results found in the literature regarding the volumetric properties of PE systems under the same conditions of temperature and pressure as those for our MC simulations.92,93 It should be noted that since at room temperature (T ) 300 K) PE is a semicrystalline material, the density value reported in ref 93 is corrected so as to correspond to purely amorphous PE, removing the effect of crystallinity on density. Simulation data in Figure 1 are fitted by hyperbolic expressions of the form

F(Nav,T,P) )

F∞(T,P)Nav a0(T,P) + Nav

(6)

where F∞(T,P) is the density value at infinite chain length and a0(T,P) is a dimensionless constant, describing the rate at which

Figure 2. Density as a function of temperature, T, as obtained from MC simulations on polydisperse C78 and C1000 linear systems. Fits of eq 7 to the simulation data are shown for comparison purposes.

F increases with Nav at the given conditions of temperature and pressure. Clearly, the hyperbolic fits of eq 6 manage to capture the qualitative trends of the simulation findings in Figure 1 at both temperatures. The corresponding pairs of fitting parameters are F∞(300K,1atm) ) 0.861 g/cm3, a0(300K,1atm)) 2.40 and F∞(450K,1atm) ) 0.775 g/cm3, a0(450K,1atm)) 3.20. The calculated density values at infinite molecular weight are in excellent agreement with the experimental measurements of Fexp(300K,1atm) ) 0.855 g/cm3 and Fexp(450K,1atm) ) 0.766 g/cm3,92,93 leading to relative errors of 0.6 and 1.1% at T ) 300 and 450 K, respectively. As pointed out in a previous publication,4 simulated densities should be as close as possible to the experimental observables; otherwise, a deviation of even 2-3% may prove critical to the correct prediction of the transport properties of the corresponding macromolecular systems. For the PE systems studied here, the application of advanced MC algorithms and the adaptation of a very accurate UA force field ensured the creation of fully equilibrated and realistic configurations mimicking the experimental amorphous samples. The dependence of density on temperature is shown in Figure 2 as obtained from MC simulations on linear, polydisperse C78 and C1000 PE samples (P ) 1atm). Assuming that the thermal expansion coefficient (β) is independent of the temperature, the corresponding density can be fit to the expression

F(T) ) F(0)e-β0T

(7)

for the whole range of applied temperatures. The resulting fitting curves are also shown in the figure, the corresponding parameters being F(0) ) 1.0482 g/cm3 and β0 ) -0.000748 K-1 for C78 and F(0)) 1.0606 g/cm3 and β0 ) -0.000692 K-1 for C1000. The square regression coefficient was over 0.998 for both fittings, pointing out the validity of the assumption that the thermal expansion coefficient (β) is independent of the temperature. Applying the definition of the volumetric thermal expansion coefficient, β

1 ∂F β)F ∂T P

( )

(8)

it follows that for the systems studied β ) β0. Therefore, according to the simulation data of Figure 2, β ) 0.000748 and 0.000692 K-1 for the C78 and C1000 PE systems, respectively. The dependence of density on molecular architecture (shortchain branching) is plotted in Figure 3. The same trend is


Figure 3. Effect of molecular architecture (short-chain branching) on density as obtained from combined MC/MD simulations on monodisperse PE systems at T ) 300 and 450 K (P ) 1 atm).

observed at T ) 300 and 450 K; by setting the linear C142 system as the reference (bearing the same MW as the short-branched counterparts), it appears that the addition of branches with a single CH3 unit (i.e., C1) causes only a slight increase in density (0.5% for T ) 300 K and 0.4% for T ) 450 K), leading to a maximum value. The C126C2 system exhibits almost identical density as the linear C142C0 analogue, whereas the addition of longer branches (C5 and C8) results in a minimal reduction in density (0.5% for C8 at T ) 300 K and 0.8% for C5 at T ) 450 K). As can be seen, even in the most extreme cases, the relative difference in the density values does not exceed 1%, lying marginally within the statistical error of the simulation method. As a consequence, our simulation data clearly suggest that shortchain branching has a very small effect on the density of the system, at least for the range of chain lengths, branch lengths and frequencies studied here. In addition, the aforementioned simulation findings suggest that the difference between the densities of short-branched (LLDPE) and linear (HDPE) PE end products arises mainly due to the different degrees of crystallinity; linear semicrystalline materials possess a larger number of crystalline domains than the short-branched analogues. III.3. Addition the Hydrogen AtomssLocal Relaxation of the Explicit-Atom Configurations. As mentioned in section III.2, the initial phase of equilibration of the studied systems was undertaken by complex chain-connectivity altering MC algorithms, providing a large set of fully relaxed and representative configurations in a united-atom (UA) description. The UA representation in the corresponding MC simulations saves a remarkable amount of computational time while this level of coarse-graining does not cause a reduction in the resolution or in the quality of the results. The reasons are the chemically simple constitution of PE and the application of a finely tuned potential force field that provides very accurate prediction of the conformational and volumetric properties, as demonstrated in section III.2. On the other hand, it has been reported in the literature that the use of an explicit-atom (EA) over UA representation improves, notably, the values of the barrier properties obtained.27 To handle this discrepancy, a specific strategy was adopted to transform the MC-generated UA structures into EA configurations compatible to the existing file formats in InsightII. Besides the technical side of implementing the existing fullequilibrated UA configurations to be simulated through the commercial packages, the physical aspect of the addition of hydrogen atoms is a key factor for the accurate prediction of

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5651 the barrier properties of the studied macromolecular materials. Hydrogen atoms in the polymer matrix exhibit enhanced mobility due to their continuous displacements and rotations. Additionally, they contribute to the dynamic flexibility of the system by providing local rearrangements of the accessible volume that play a significant role, especially in the melt phase where penetrants are “carried along” by density fluctuations caused by the thermal motion of the polymer segments.3 The MC-generated structures and the environment of the commercial software were integrated through a developed homemade algorithm, where the UA structures were modified by inserting the corresponding hydrogen atoms to the sp3 carbon atoms at which they were bonded. The charges where automatically assigned by the atom typing rules when the COMPASS force field (condensed-phase optimized molecular potentials for atomistic simulation studies) was loaded. This force field was chosen as it is considered nowadays as one of the most accurate, general-purpose force fields.94,95 It is expected that the addition of hydrogen atoms may cause some interatomic overlaps between an existing carbon atom and an introduced hydrogen atom or between added hydrogens. These overlaps should be removed before the final stage of TST as these topological configurations result in unrealistic local packing of the polymer segments. Therefore, an additional local relaxation procedure is demanded to energetically “clean” the newly created EA structure. A protocol consisting of a succession of static energy minimizations at constant density (molecular mechanics, MM)96 and short NVT MD simulations were implemented using the Discover module of the commercial package Materials Studio.97 The static energy minimization was initiated by a steepest descent method until derivatives were less than 100 kcal mol-1. Afterward, a Polak-Ribiere conjugate gradient algorithm was applied, allowing a deviation of the energy derivatives of 1 kcal mol-1. The total duration of the NVT MD cycle was limited to 50 ps, the integration time step was set equal to 1 fs, and the temperature was controlled by the Berendsen thermostat. In all cases, the NVT statistical ensemble was selected as the density, which resulted from the preceding MC simulations, was in excellent agreement with available experimental data, a critical requirement for the successful prediction of the barrier properties. The number of the required cycles (i.e., iterations of the same procedure of static energy minimizations and NVT MD) was determined by the evolution of the potential energy during preliminary execution of a large number of minimization/MD pairs. Convergence of the total energy to a plateau value indicates that the segments of the EA system reconfigure locally so that their configuration is compatible with the applied force field (COMPASS). Figure 4 shows the time evolution of the potential energy of an EA system, C1000 and T ) 300 K, as obtained from the application of 6 NVT MD simulations with corresponding static energy minimizations executed before the initiation of each MD step, as indicated by the change in line style in the potential energy every 50 ps. As seen in Figure 4, potential energy stabilizes and oscillates around an average value during the third MD step; therefore, further local relaxation of the system can be considered as unnecessary. As mentioned above, it is essential to verify the impact of the changes in the force field and representation (from UA to EA) on the quality of the resulting configurations and how the new relaxation cycle, by combining MM and short MD steps, rectifies any unrealistic local packing caused by this transformation. The carbon-carbon intermolecular and intramolecular pair distribution functions, g(r) and w(r), for the initial configuration


Gestoso and Karayiannis TABLE 3: Collision Diameter, σ, and Well Depth, E, of the Lennard-Jones Potential for the Gas Penetrant Molecules in Spherical, United-Atom Representations as Implemented in Gusev-Suter TST and the van der Waals Molecular Surface Area (WSA) and the Effective Well Depth (E/kB)eff Parameter Used in the Yampolskii’s101 and Teplyakov’s100 Correlation Formulas, Respectively

Figure 4. Potential energy as a function of time as obtained from successive 50 ps NVT MD simulations for linear C1000 melts at T ) 300 K using the COMPASS force field. Before the application of each MD simulation, the system undergoes an additional step of static energy minimization.

(after the addition of hydrogen atoms and before any minimization) and those obtained at the end of each successive MM/ MD cycle were compared, the results indicating that the intermolecular and intramolecular C-C packing, which resulted from the application of the MM/MD cycles, is very similar to that obtained from the MC equilibration. This suggests that there is no appreciable perturbation of the configuration due to the application of a different force field and the addition of the hydrogen atoms. From these results, it is concluded that after the application of two combined MM/MD cycles, the EA configuration is fully relaxed through the application of the COMPASS force field and that the barrier properties can be safely calculated by the application of the Gusev-Suter TST after the end of the third relaxation cycle. For each polymer structure selected, the barrier properties were calculated after the third and fourth minimization/NVT MD cycles. III.4. TST Application. TST method, as formulated by Arrizi24 and Gusev et al.,25-27 was invoked through the gsnet and gsdiff subroutines in InsightII. The sorption sites (i.e., clusters of accessible volume for a given penetrant), their connectivity, and the rate constants for the transitions between them were identified by the gsnet module. Initially, a 3-D grid was laid on the fully equilibrated simulation box, covering the whole volume of it. The spacing of the grid was set at 0.5 Å, corresponding to a resolution slightly lower than that reported for past TST calculations (around 0.3 Å).4,48-58,61 The simulation boxes, created in the course of the MC-based equilibration, were much larger (40-45 Å) than the cell dimensions (25-30 Å) of a typical detailed atomistic simulation for the study of the barrier properties.4 Current implementation of the TST in InsightII’s gsnet and gsdiff subroutines prevents the use of cells larger than 40-45 Å since there is a limit in the maximum number of allocated sites comprising the 3-D grid. Several test simulations were carried out to ensure that the selected grid spacing of 0.5 did not affect the smearing factor and the successive prediction of solubility and diffusivity coefficients. It was further concluded that a grid spacing of 0.5 Å did not introduce artificial estimates of the measured quantities for the PE matrix/penetrant molecule combinations simulated here. As mentioned in section II, in the current implementation of TST, all penetrant molecules are represented as single, spherical,

penetrant

σ (Å)

(kcal/mol)

WSA (Å2/molecule)

(/kB)eff (K)

O2 N2 CH4 Ar CO2

3.46 3.70 3.82 3.84 4.00

0.2344 0.1889 0.2945 0.2464 0.4500

39.6 40.9 48.0 44.4 52.4

112.7 83.0 154.7 122.3 213.4

united-atom sites whose short-range interactions with the polymer atoms are described through a 9-6 Lennard-Jones (LJ) potential. The collision diameters, σ, and the well depths, , are listed in Table 3. In all cases, the nonbonded interactions were truncated at 9.5 Å as test runs with larger cutoff distances revealed no change in the calculated barrier properties. After applying the 3-D grid, the smearing factor, ∆2, was calculated for each penetrant through a self-consistent (SC) scheme involving information about the mean square displacement (MSD) of the polymer atoms from their equilibrium positions. The MSD curve was obtained from an additional 50 ps NVT MD simulation executed after the completion of the final relaxation cycle. The two factors entering the self-consistent scheme were the smearing factor, ∆2, and the most probable residence time of the penetrant in the sorption sites. The SC method converged when the relative difference between two successive ∆2 values was less than 2.5%. The resulting value depends not only on the polymer matrix and the applied temperature but also on the penetrant molecule. More details about the adopted self-consistent scheme for the calculation of the smearing factor can be found elsewhere.4 On the basis of the above procedure, the gsnet program allows the identification of the free and accessible volume, which is used to calculate the rate constants of all of the transitions between sorption sites and subsequently allows the prediction of the low-concentration solubility coefficient, S, based on the Widom method67,68 (eqs 1 and 2). This information is subsequently used by the gsdiff program, which uses a kinetic MC (kMC) scheme to track the MSD of the penetrant as a function of time and to calculate the low-concentration diffusion coefficient, D, in the Fickian limit (eq 4). The total duration of the kMC was 10-4 s, and the MSD was averaged over the trajectories of 1000 penetrant walkers. As pointed out by previous simulations efforts,4 for the accurate prediction of the barrier properties, it is important to sample over as many different configurations as computational time and resources permit. This allows the calculation of the S and D values as statistical averages over the whole available set of structures. The present multiscale modeling methodology fits perfectly to this approach. For example, the proposed MC technique provides a fully equilibrated frame for the linear, polydisperse C1000 PE system at T ) 300 K within less than 42 h of CPU time (the required time drops to just 2 h if the same system is simulated at T ) 450 K), the additional relaxation procedure (to integrate MC-generated structures with the InsightII software and transform from UA to EA representations) consumes approximately 34.9 CPU hours (Pentium IV, 3.0 Hz) for the same system, and finally, Gusev-Suter TST requires an average of 0.2 CPU hours (Origin 200 R12000 360 MHz) to calculate S and D. In total, by combining the ultrafast MCbased equilibration method and the successive TST calculation, the barrier properties of an entangled C1000 configuration to small



TABLE 4: Simulation Results for the Low-Concentration Solubility, S, for Different Penetrant Molecules in a Linear C1000 PE Matrix at T ) 300 K and Experimental Data, Corrected for Crystallinity, from Reference 98

TABLE 5: Simulation Results for the Low-Concentration Diffusivity, D, for Different Penetrant Molecules in a Linear C1000 PE Matrix at T ) 300 K and Experimental Data, Corrected for Crystallinity, from Reference 98

penetrant

S(simulation) (10-6 cm3 (STP)/cm3 Pa)

S(experimental)a (10-6 cm3 (STP)/cm3 Pa)

penetrant

D(simulation) (10-7 cm2/s)

D(experimental) (10-7 cm2/s)

O2 N2 CH4 Ar CO2

1.14 ( 0.08 0.418 ( 0.08 1.18 ( 0.16 0.570 ( 0.1 3.20 ( 0.3

0.828 0.400 1.97 1.00 4.47

O2 N2 CH4 Ar CO2

13.2 ( 1.2 4.84 ( 0.8 1.60 ( 0.05 1.80 ( 0.1 0.41 ( 0.1

12.1 8.42 5.08 9.47 9.79

a Solubility calculated as S ) P/D from available diffusivity and permeability data.98

gas molecules at T ) 300 K can be predicted within less than 77.1 CPU hours. This extraordinary feature of the hierarchical approach allowed the extensive analysis of a large number of totally different and realistic PE configurations at different temperatures, resulting in a significant reduction in the statistical error of the barrier properties’ calculations. Ideally, the whole equilibrated NPT/MC trajectory can be used as an expanded ensemble for TST predictions, but due to CPU limitations, a significantly smaller set of representative configurations was used, depending on the applied temperature. According to preliminary simulations, sampling over 4 structures was adequate at T g 500 K, over 5 in the range of 500 K > T g 350 K, and over 8-10 at lower temperatures (350 K > T g 250 K). Reported errors bars correspond to standard deviations of the distributions of the measured quantities. IV. TST Results IV.1. Comparison with Experiments and Correlation Functions. Before attempting to explore and identify the effect of temperature, molecular weight, and architecture on the barrier properties of model PE systems, it was important to validate the simulation results obtained with the proposed multiscale methodology (MC combined with TST) against available experimental data. It should be noted that at standard conditions (T ) 300 K and P ) 1 atm), the vast majority of experimental data correspond to semicrystalline PE materials, whereas our computer-generated PE structures are purely amorphous. On the basis of the assumptions that (a) the penetrant molecules cannot be absorbed or diffused through crystalline regimes and (b) the presence of a crystalline phase does not affect the density of the amorphous regimes (i.e., no densification of the amorphous volume occurs),2 a lower bound estimate for the “corrected” experimental diffusivity, D′exp, so as to correspond to purely amorphous (ideal) systems, is provided by Pant and Boyd32 through

D′exp )

3 Dexp 2 (1 - Vφ)

(9)

where Dexp is the reported experimental diffusion coefficient of the penetrant in the semicrystalline material and Vφ is the crystalline volume fraction. Tables 4 and 5 present S and D, respectively, as obtained by the present simulation methodology for the molecules listed in Table 3 in a polydisperse, linear C1000 system at T ) 300 K. The experimental data in the tables correspond to the values published by Michaels and Bixler98 for C26000 (T ) 298 K, Vφ ) 0.43) corrected by the crystallinity percentage.32 The reported experimental solubility coefficient for the gas molecules in PE was calculated indirectly as S ) P/D based on diffusivity and permeability data of ref 98. In all cases, simulated coefficients

are obtained as statistical mean values averaged over eight different configurations. In general, there is a good agreement between experimental data and simulation findings for solubility and diffusivity, especially for the relatively small penetrants (O2, N2, and CH4). For the low-concentration solubility, our simulation predictions are in very good agreement with the experiments even for the larger gases (Ar, CO2) that do not fulfill the requirements and simplifications invoked by the Gusev-Suter TST method. For all samples, simulated solubilities lay within a factor of 2 (in most cases even better) from experimental data. In addition, TST-based simulations predict correctly the experimental trends; in the extreme behaviors, CO2 and N2 show the highest and lowest solubility in the PE matrix, respectively. In general, the sorption coefficient of gases and vapors increases with increasing condensability of the penetrant. Thus, sorption selectivity favors large, more condensable molecules. Several expressions have been proposed trying to correlate condensability with the boiling temperature (Tb), the critical temperature (Tc),99 the well depth parameter (/kΒ) of the L-J potential,100 and the molecular surface area101 of the dissolved species. According to Yampolskii and co-workers,101 the solubility coefficient is correlated with the van der Waals surface area (WSA) of the penetrant molecule through

ln(S) ) a0 + a1WSA

(10)

where a0 and a1 are the regression coefficients. The values for the van der Waals molecular surface areas (WSA) used for the gas molecules studied here can be found in Table 3. Equation 10 was fitted to both simulation findings and experimental data for solubility coefficients (Table 4) as a function of the type of penetrant quantified through WSA. As seen in Figure 5, Yampolskii’s correlation functions capture very accurately both the qualitative and quantitative solubility trends exhibited in PE, except in the case of oxygen (O2), where significant deviation from linearity in the ln(S) versus WSA scaling is observed for both simulation and experimental data. This is a rather unexpected result since oxygen’s WSA parameter was proven quite accurate in the prediction of solubility in a large set of different polymeric systems.101 From fitting the experimental data to Yampolskii’s formula, we obtained a0(experimental) ) -23.1 and a1(experimental) ) 0.206, whereas values of a0(simulation) ) -22.3 and a1(simulation) ) 0.181 were calculated using our simulation results. The fitting was obtained excluding oxygen from the analysis of the solubility correlation for the amorphous PE at T ) 300 K. Correlations for the solubility coefficients of gases can be also established by determining an effective Lennard-Jones well depth, (/kB)eff, for each penetrant as proposed by Teplyakov and co-workers100

log(S) ) K3 + K4

() kB

eff

(11)



Figure 5. Correlation of ln(S) (in cm3 (STP)/cm3 Pa) with the van der Waals molecular surface (WSA) using simulation and experimental data for the low-concentration solubility, S, of gas penetrants in amorphous PE at T ) 300 K. Also shown are corresponding linear regressions (excluding the predictions for oxygen) based on Yampolskii’s correlation function (eq 10).

Figure 7. Low-concentration solubility, S, of O2, N2, and CH4 as a function of molecular length as obtained from present TST simulations on linear polydisperse PE samples at (a) T ) 450 K and (b) T ) 300K (also shown are hyperbolic fits of the simulation data).

Figure 6. Correlation of log(S) (in cm3 (STP)/cm3 Pa) with the well depth parameter of the Lennard-Jones potential, /kB, for the lowconcentration solubility, S, as obtained from present TST-based simulations and experimental data for purely amorphous PE (T ) 300 K). Also shown is the best linear fit of the experimental predictions according to the Teplyakov correlation function (eq 11).

where K3 and K4 are the regression coefficients. This correlation (eq 11) is reminiscent of the linear dependence on the penetrant’s heat of vaporization as suggested by regular solution theory.18 Table 3 lists the effective well depth parameters (/kB)eff for the gas molecules studied here. Figure 6 presents the dependence of the logarithm of the solubility coefficient on the L-J well depth of the gas penetrant molecule as obtained from our simulation and the available experimental data listed in Table 4. As in the case of the Yampolskii correlation function, the best linear fit of the Teplyakov analogue (eq 11) follows closely the experimental solubility data for all existing gas penetrants providing the values for the regression coefficients K3 ) -6.99 and K4 ) 0.788 × 10-2. TST-based predictions seem to scatter around the regression line in a rather reasonable agreement. Previous simulation efforts suggested that solubility coefficients, especially for large gases (CH4 and CO2), are strongly influenced by system size effects, and whenever small simulation boxes are utilized, solubility coefficients deviate significantly from experimental data.61,102,103 In the present simulations, the dimensions of the MC-generated cells are almost twice as large

as the ones used in previous TST studies,4 ruling out the possibility of simulation artifacts. Table 5 summarizes the diffusion coefficient results, D, for C1000 and T ) 300 K. TST simulations capture both the qualitative and quantitative diffusion trends for gas penetrants as large as methane (CH4). As for O2 and N2, the agreement between simulation and experiment is very good. Clearly, because of its simplifying assumption about the spherical representation of the gas penetrants (especially for CO2) and the rather coarse implementation through the smearing factor ∆2 of the structural relaxation of the polymer atoms in the presence of large dissolved entities, TST is unable to reliably predict the transport properties of larger gas molecules. This is a well-established observation regarding the present formulation of the Gusev-Suter TST, with CH4 being widely considered as the gas penetrant molecule that lies marginally within its effective range of application. Further analysis of the diffusion data for the molecules that are tractable by TST (i.e., O2, N2, and CH4) suggests that the smaller the molecule size, the larger the diffusion coefficient of the penetrant, in agreement with past experimental data98 and simulation trends42,55,58,61,104,105 for the transport of small gases through polymer matrixes. IV.2. Effect of Molecular Weight. Figure 7a and b presents the low-concentration solubility, S, of O2, N2, and CH4 as a function of molecular weight (MW) of linear polyethylene at T ) 450 and 300 K, respectively. Also shown in Figure 7b are the corresponding experimental solubility data for linear PE (corrected for crystallinity, as listed in Table 4).39

Polyethylene Barrier Properties At T ) 300 K, S appears to be independent of the chain length, within the standard error, and follows the trend S(CH4) > S(O2) > S(N2), in agreement with the experimental data. Conversely, at T ) 450 K, for all gas penetrants, S decreases hyperbolically with MW in a fashion reminiscent of the density hyperbolic increment as shown in Figure 1; as density increases with MW, the portion of the free volume (i.e., the volume of the simulation box not occupied by polymer segments) decreases, and consequently, the clusters of accessible volume (i.e., sorption sites) are reduced. Nevertheless, it is important to point out that the increase in density for chain lengths between C78 and C1000 is less pronounced (3.2%) than the increase of solubility for the different gases (12.8% for O2, 16.5% for N2, and 19.3% for CH4). The observed solubility trend for T ) 450 K suggests that S(O2) > S(CH4) > S(N2). The change of behavior with temperature is not exceptional; experimental solubility data for H2 and N2 in polybutadiene106 showed that, whereas at low temperatures (270 K) the solubility of H2 is higher than that for N2 (following the correlation functions as proposed by Yampolskii101 and Teplyakov100), for higher temperatures (e.g., 330 K), the trend is reversed. This point will be further discussed in section IV.4. These results suggest that solubility behavior with MW obeys different mechanisms depending on the temperature. At high temperatures (melt state), MW has a profound effect on solubility, whereas at standard conditions, this impact diminishes. Additionally, temperature appears to affect the solubility selectivity (S(i)/S(j)) of PE regarding the sorption of gas molecules. This will be discussed further in section IV.4. Figure 8a and b shows the dependence of low-concentration diffusivity, D, on MW at T ) 450 and 300 K. At standard conditions, TST-based diffusivity results lie very close to the experimental values, especially for oxygen. In addition, excellent agreement is also observed at T ) 450 K for nitrogen (no experimental data were available for O2 and CH4) in the melt state,107 with the relative error being less than 10%. GusevSuter TST was developed primarily to model permeation in dense, glassy polymeric systems, where the characteristic time frame of the related transport phenomena is out of reach for the all-conventional MD methods. In contrast to the thermodynamically driven phenomenon of sorption, kinetically driven diffusion is strongly correlated to the density (and consequently free volume) of the polymer matrix, as clearly indicated by TST results at both temperatures, in agreement with the observations of Charati and Stern102 and Lim et al.108 As in the case of solubility, variation in the diffusion coefficient is more significant than that in density for chain lengths between C78 and C1000; whereas at T ) 300 and 450 K the increase in density is 2.3 and 3.2%, respectively, the decrease in the diffusion coefficient is 64.5 (O2), 72.0 (N2), and 79.4% (CH4) at 300 K and 28.2 (O2), 28.4 (N2), and 29.1% (CH4) at 450 K. Besides the fraction and distribution of free volume, additional mechanisms affect gas transport in melt and rubber PE; as demonstrated in Figure 1, linear C500 and C1000 PE systems exhibit practically identical densities, whereas Figure 8 suggests that for all penetrants in the temperature range studied here, low-concentration D in C1000 is systematically smaller than that in C500. This reveals the effect of local relaxation of the polymer atoms (segmental dynamics) on the gas diffusion process. The effect of MW on diffusion is amplified as temperature decreases; at T ) 450 K, the diffusivity ratio (DC78(i)/DC1000(i)) is similar for all penetrants around a value of 1.4, whereas at T ) 300 K, this ratio increases up to 2.8, 3.6, and


Figure 8. Low-concentration diffusivity, D, of O2, N2, and CH4 as a function of molecular length as obtained from present TST simulations on linear polydisperse PE samples at (a) T ) 450 K and (b) T ) 300 K.

4.9 for O2, N2, and CH4, respectively. Even through the incorporation of the simplified smearing factor, a rather coarse quantification of the segmental polymer motion, the TST method is capable of distinguishing between the different dynamical flexibility of the polymer host as a consequence of the different molecular weight. IV.3. Effect of Polydispersity. The effect of polydispersity is explored through the simulation of the barrier properties of two linear C142 systems characterized by the same average molecular weight, one being strictly monodisperse (I ) 1.000) and the other with uniformly distributed chain lengths ranging from C71 up to C213 corresponding to a polydispersity index of I ) 1.083. MC simulations, used for the full-scale equilibration of both samples, proved that polydispersity has no appreciable effect on density as long as the mean chain length remains constant. Table 6 summarizes the calculated D and S coefficients for all tested penetrants (O2, N2, and CH4) through the monoand polydisperse C142 PE systems at T ) 300 and 450 K. The results suggest that polydispersity does not influence the barrier properties of the polymer matrix, at least for the range of molecular length distributions studied here. IV.4. Effect of Temperature. The effect of temperature on the barrier properties of linear, polydisperse PE systems resulting from our TST calculations is depicted in Figures 9 and 10. Figure 9a and b presents the logarithm of solubility, log(S), of O2, N2, and CH4 as a function of reciprocal temperature for the C78 and C1000 melts, respectively, in the range 250 K e T e 600 K. Penetrant solubility decreases with temperature in all



TABLE 6: Solubility (S) and Diffusivity (D) for O2, N2, and CH4 as Obtained from TST Simulations on Monodisperse (I ) 1.000) and Polydisperse (I ) 1.083) Linear C142 Systems at T ) 300 and 450 K S (10-7 cm3 (STP)/cm3 Pa)

D (10-6 cm2/s)

system

O2

N2

CH4

O2

N2

CH4

I ) 1.083, T ) 300 K I ) 1.000, T ) 300 K I ) 1.083, T ) 450 K I ) 1.000, T ) 450 K

11.4 ( 0.8

4.30 ( 0.08

12 ( 2

2.5 ( 0.1

1.08 ( 0.08

0.42 ( 0.01

10.8 ( 0.9

4.0 ( 0.1

11.2 ( 1

2.5 ( 0.1

0.98 ( 0.06

0.37 ( 0.03

1.7 ( 0.1

0.83 ( 0.05

1.25 ( 0.08

113 ( 10

93 ( 8

72 ( 8

1.80 ( 0.09

0.87 ( 0.06

1.33 ( 0.08

114 ( 10

92 ( 6

68 ( 9

PE matrixes as the gas molecule experiences more difficulties in condensing the polymer when the temperature is raised. Although some polymer systems show an increase of gas solubility with temperature,109 the common rule suggests that solubility decreases with temperature, as confirmed by the experimental data of Mi et al.110 on the solubility of CO2 in PET and several simulation works (CO2 and He in PE,111 CO2 and CH4 in polyetherimide,108 CH4 and CO2 in HDPE,112 and n-alkanes in poly(dimethylsilamethylene)).43 It is also interesting to notice that whereas O2 and N2 behavior with temperature is very similar, in the case of methane, the trend is different. At high temperatures, solubility values for methane values are very close to N2. As temperature decreases, the solubility values increase, approaching O2 values. Finally, for temperatures below 300 K, the solubility is higher than that for O2. This difference

for the methane, which is the largest molecule studied, is in agreement with previous reports in the literature, pointing to a more pronounced solubility decrease with temperature for large gas molecules.43,113 Figure 10a and b hosts the dependence of log(D) on the inverse of temperature for C78 and C1000, respectively. For the entangled polymeric system (C1000), diffusivity results for methane are compared against available experimental data both in the melt114 and at standard conditions;98 in all cases, the quantitative agreement is very good. The combination of exhaustive pre-equilibration using chain-connectivity altering MC moves and the application of Gusev-Suter TST is proven to be relatively accurate even at high temperatures where penetrant diffusion is highly coupled with the relaxation and mobility exhibited by the atoms of the host polymer. In order

Figure 9. Logarithm of solubility, S, (in cm3 (STP)/cm3 Pa) of O2, N2, and CH4 as a function of reciprocal temperature as obtained from TST simulation for (a) C78 and (b) C1000 linear PE systems.

Figure 10. Logarithm of diffusivity, D, (in cm2/s) of O2, N2, and CH4 as a function of reciprocal temperature as obtained from TST simulation for (a) C78 and (b) C1000 linear PE systems.



TABLE 7: Solubility (S(i)/S(j)), Diffusivity (D(i)/D(j)), and Permeability (P(i)/P(j)) Selectivities between the Pairs N2-O2 and N2-CH4 as a Function of Temperature, T, for Linear C1000 PE T (K)

S(O2)/ S(N2)

S(O2)/ S(CH4)

D(O2)/ D(N2)

D(O2)/ D(CH4)

P(O2)/ P(N2)

P(O2)/ P(CH4)

250 300 350 400 450 500 550 600

2.98 2.73 2.44 2.28 2.14 2.00 1.91 1.82

0.684 0.966 1.12 1.28 1.42 1.47 1.52 1.55

2.71 1.64 1.37 1.27 1.20 1.12 1.10

8.25 3.05 1.98 1.65 1.49 1.37 1.35

7.44 3.99 3.11 2.70 2.41 2.14 2.01

7.97 3.41 2.53 2.34 2.18 2.08 2.10

to explain the differences between C78 and C1000, it is useful to see the dependence of density on temperature for both chain lengths (Figure 1). From the comparison of both figures, we can see that the decrease in diffusivity with the increase of chain length is directly related to the decrease in density, which was already mentioned in previous sections. As shown in Figure 10, the dependence of log(D) on reciprocal temperature for all penetrants (O2, N2, and CH4) follows a non-Arrhenius behavior. This distinct departure from Arrhenius behavior was originally identified by Pant and Boyd for CH4 diffusion in PE matrixes.32 Present TST data suggest that this temperature dependence holds also for oxygen and nitrogen transport in amorphous PE, in line with previous experimental work in methane through polystyrene109 and MD simulations for CO2 through polyethylene111 and CH4 in polybutadiene matrixes.47 At the lowest temperature (T ) 250 K), methane’s diffusivity is slower than that in the melt state (T ) 450 K), adopting a value of approximately 1.5 × 10-9 cm2 s-1 laying out of the range of MD methods. This sharp slow down in the penetrant’s kinetics unfolds the onset of the polymer’s transition to the glassy state, rendering gas diffusion a process of rare jump events (transitions) between accessible sites with the gas molecules spending most of their time trapped in the sorption clusters. For different pairs of gases (i and j), one can extract solubility (si/j ) S(i)/S(j)) and diffusivity (di/j ) D(i)/D(j)) selectivities and consequently perm-selectivity through

ai/j )

( )( )

S(i) D(i) P(i) ) si/jdi/j ) P(j) S(j) D(j)

(12)

Table 7 lists solubility, diffusivity, and permeability selectivities as a function of temperature using oxygen as the reference molecule. For the solubility selectivity of the O2/N2 pair, an increased preference to oxygen is observed as temperature decreases, whereas for the O2/CH4 pair in the melt state, O2 is more soluble in PE by a factor of around 1.5, but as temperature decreases, the trend changes, favoring the presence of CH4. At room temperature, simulation findings indicate that there is no preference in the sorption mechanisms between the two gases (O2 and CH4) as sO2/CH4 is approximately unity, in contrast to the experimental ratio of sO2/CH4 ≈ 0.42.98 This discrepancy is attributed to the TST-based under- and overestimation of the sorption coefficients for O2 and CH4, respectively. In general, present results point out that the solubility of large molecules (CH4) is favored as temperature decreases, in agreement with previous literature reports proving that a solubility decrease with temperature is more pronounced as the size of the penetrant increases.43,114 Diffusivity selectivity for both pairs (O2/N2 and O2/CH4) obviously reflect the enhanced transport behavior of gases with

Figure 11. Low-concentration solubility, S, of O2, N2, and CH4 as a function of molecular architecture (linear and short-chain branched) as obtained from TST simulations at (a) T ) 450 K and (b) T ) 300 K.

reduced size as temperature decreases. At standard conditions, oxygen diffusivity is 2.7 and 8.3 faster than the corresponding ones of N2 and CH4, respectively. Perm-selectivity, as calculated from eq 12, is shown in Table 7 as well. At room temperature, according to TST simulation data, the permeation of oxygen in amorphous linear PE is preferred to nitrogen and methane by similar ratios of 7.4 and 8.0, respectively. The simulated oxygen-nitrogen permselectivity, which steams equally from kinetic (diffusion) and thermodynamic (sorption) contributions, of aO2/N2 ) 7.4 compares reasonably well with the available experimental value98 of aO2/N2 ) 3.0. On the opposite, TST-based oxygen-methane perm-selectivity, aO2/CH4 ) 8.0, cannot capture the experimental qualitative and quantitative trend according to which oxygen and methane permeation in PE is equally favored. The cause for this discrepancy is the indirect calculation of permeability as the product of DS, which bears the accumulation of errors that are introduced in the calculation of solubility and diffusivity, especially for CH4 whose transport behavior can be marginally captured by Gusev-Suter TST. IV.5. Effect of Short-Chain Branching. Thoroughly equilibrated short-chain branched (SCB) structures are generated to study the effect of molecular architecture on the barrier properties of amorphous PE. All comparisons are made with the linear monodisperse C142C0 as the reference system. Four different monodisperse SCB analogues are chosen, with pre-


Gestoso and Karayiannis diffusion coefficient of O2 at 300 K (Figure 12b) presents a somewhat different behavior than the other penentrants at the same temperature. This is a consequence of the difficulty of the algorithm to converge for this penentrant at low temperature, causing fewer values to be obtained for each branch length and producing more scatter as a result. In contrast to the small increment recorded in the solubility coefficient of SCB structures against their linear analogues, diffusivity is seen to drop with the appearance of short branches along the main backbone, with the effect being more pronounced at standard conditions. At T ) 300 K, a minimum in the transport coefficient is related to the SCB system bearing the smallest branches (C1) (34.3% for O2, 36.1% for N2, and 46.5% for CH4), which is attributed to the maximum in the density exhibited by the same system (C132C1). However, at T ) 450 K, the minimum is observed at C5 (16.9% for O2, 19.0% for N2, and 18.5% for CH4), corresponding to the minimum in density exhibited by the same melt C112C5. Those results point out that the effect of short branching is amplified at low temperatures. V. Conclusions

Figure 12. Low-concentration diffusivity, D, of O2, N2, and CH4 as a function of molecular architecture (linear and short-chain branched) as obtained from TST simulations at (a) T ) 450 K and (b) T ) 300 K.

scribed main backbone, branch length, frequency, and number of branches so as to correspond to a total MW equal to C142C0. More details about these systems can be found in Table 2. As already mentioned in section III.2 and depicted in Figure 3, the density appears to be, in general, independent of molecular architecture, a slight increase only observed for C132C1, in agreement with past experimental studies of linear and SCB PE melts.114 Branch lengths of the SCB configurations range from C1 (single CH3 unit) up to C8, with the corresponding number of branches per chain varying from a maximum of 10 down to a minimum of 4. The dependence of solubility for O2, N2, and CH4 on molecular architecture (linear vs SCB) at T ) 450 and 300 K is depicted in Figure 11a and b, respectively. In general, shortchain branching appears to have a small effect on S for all gas penetrants, with a maximum value of solubility being recorded for C112C5. The observed increment is rather independent of temperature and accounts for an increase in the values of around 23% for oxygen and 30% for nitrogen and methane. Gas solubilities for all other SCB systems are very similar, slightly larger than that for the linear counterpart characterized by the same MW. These results are in agreement with by experimental studies on various PE systems.114 Figure 12a and b exhibits the trends of low-concentration diffusivity at T ) 450 and 300 K, respectively, for linear and short-chain branched PE systems. It can be noticed that the

We have presented simulation results for the sorption, diffusion, and selectivity of small gas molecules (O2, N2, Ar, CO2, and CH4) through PE matrixes characterized by different chain lengths and molecular architectures at infinite dilution for a wide range of temperatures. A multiscale, hierarchical modeling approach was adopted, which is based on two main steps. First, advanced chain-connectivity altering moves undertook the full-scale equilibration of purely amorphous PE systems at various conditions. Then, a number of selected representative structures served as initial configurations for the study of the barrier properties through the application of the Gusev-Suter transition state theory (TST) as implemented in the commercial software InsightII. Present simulation findings for the barrier properties of linear PE systems are in good agreement with the available experimental data. Preliminary TST-based calculations suggested that the method provides accurate predictions for molecules as large as CH4. The solubility selectivity results point to a preferential solubility of large molecules as the temperature decreases. On the other side, the diffusion of smaller gases is enhanced as the temperature decreases. Low-concentration diffusivity of all tested molecules depends on temperature through a highly non-Arrhenius function. On the other hand, solubility is found to increase with decreasing temperature; as the temperature increases, the penetrant experiences more difficulties to condense in the polymer matrix. Regarding the effect of molecular architecture, the results suggest that the existence of short branches along the main backbone has a small effect on the barrier properties of the PE system, as long as the comparison is made on the same MW and for purely amorphous samples. The results from this work show that density plays a dominant role in the barrier properties, especially on the diffusivity coefficient. It is important to point out that the changes in density appear to be correlated with more significant variations in the barrier properties, particularly at low temperatures. Finally, the strength of the Gusev-Suter TST method is proven as it manages to capture small changes in the transport behavior in two polymer matrixes characterized by the same density (C500 and C1000) but exhibiting different segmental and global dynamics.

Polyethylene Barrier Properties Future plans involve the extension of the proposed simulation methodology to more complex nonlinear PE systems (H-shaped, stars), their blends with linear chains, and the extension to chemically more complex macromolecular structure (PI, PIB, PP). In addition, the proposed scheme is extended to study polymer/penetrant systems well below the glass transition temperature. Acknowledgment. P. Gestoso acknowledges Rhodia Recherches for permission to publish this work. The authors thank Dr. James Wescott (Accelrys Ltd., Cambridge, U.K.) for his assistance throughout the entire project on the implementation of TST subroutines. Fruitful discussions with Prof. Manuel Laso (University of Madrid), Prof. Vlasis Mavrantzas (University of Patras), and Prof. Rafiqul Gani and Vipasha Soni (Technical University of Denmark) are deeply appreciated. N. Ch. Karayannis is indebted to Dow-Benelux B. V. for the financial support during the years of 2002-2005. This work was financially supported by the “Polymer Molecular Modeling at Integrated Length/time Scales” (PMILS) European Community project (Contract No. G5RD-CT-2002-00720). References and Notes (1) Zhao, J.; Wang, J. J.; Li, C. X.; Fan, Q. R. Macromolecules 2002, 35, 3097. (2) Liu, R. Y. F.; Hu, Y. S.; Hibbs, M. R.; Collard, D. M.; Schiraldi, D. A.; Hiltner, A.; Baer, E. J. Appl. Polym. Sci. 2005, 98, 1615. (3) Theodorou, D. N. In Diffusion in Polymers; Neogi, P., Ed.; Marcel Dekker: New York, 1996. (4) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. Macromolecules 2004, 37, 2978. (5) Wang, X. Y.; Veld, P. J. I.; Lu, Y.; Freeman, B. D.; Sanchez, I. C. Polymer 2005, 46, 9155. (6) Vrentas, J. S.; Duda, J. L. J. Polym. Sci., Polym. Phys. Ed. 1977, 15, 403. (7) Vrentas, J. S.; Vrentas, C. M. Macromolecules 1994, 27, 5570. (8) Vrentas, J. S.; Vrentas, C. M.; Faridi, N. Macromolecules 1996, 29, 3272. (9) Vieth, W. R.; Sladek, K. J. J. Colloid Sci. 1965, 20, 1014. (10) Petropoulos, J. H. J. Polym. Sci., Part A2 1970, 8, 1797. (11) Paul, D. R.; Koros, W. J. J. Polym. Sci., Polym. Phys. Ed. 1976, 14, 675. (12) Petropoulos, J. H. J. Membr. Sci. 1990, 53, 229. (13) Brandt, W. W. J. Phys. Chem. 1959, 63, 1080. (14) Pace, R. J.; Datyner, A. J. Polym. Sci., Polym. Phys. Ed. 1979, 17, 436. (15) Pace, R. J.; Datyner, A. J. Polym. Sci., Polym. Phys. Ed. 1979, 17, 453. (16) Pace, R. J.; Datyner, A. J. Polym. Sci., Polym. Phys. Ed. 1979, 17, 465. (17) Neogi, P. J. Polym. Sci., Polym. Phys. Ed. 1993, 31, 699. (18) Petropoulos, J. H. In Polymeric Gas Separation Membranes; Paul, D. R., Yampolski, Y. P., Eds.; CRC Press: Boca Raton, FL, 1994. (19) Pant, P. V. K.; Theodorou, D. N. Macromolecules 1995, 28, 7224. (20) Mavrantzas, V. G.; Boone, T. D.; Zervopoulou, E.; Theodorou, D. N. Macromolecules 1999, 32, 5072. (21) Karayiannis, N. Ch.; Mavrantzas, V. G.; Theodorou, D. N. Phys. ReV. Lett. 2002, 88, 105503. (22) Karayiannis, N. Ch.; Giannousaki, A. E.; Mavrantzas, V. G.; Theodorou, D. N. J. Chem. Phys. 2002, 117, 5465. (23) Karayiannis, N. Ch.; Giannousaki, A. E.; Mavrantzas, V. G. J. Chem. Phys. 2003, 118, 2451. (24) Arizzi S. Diffusion of Small Molecules in Polymeric Glasses: A Modeling Approach. Ph.D. Thesis, Massachusetts Institute of Technology, Boston, 1990. (25) Gusev, A. A.; Arizzi, S.; Suter, U. W.; Moll, D. J. J. Chem. Phys. 1993, 99, 2221. (26) Gusev, A. A.; Suter, U. W. J. Chem. Phys. 1993, 99, 2228. (27) Gusev, A. A.; Mu¨ller-Plathe, F.; van Gunsteren, W. F.; Suter, U. W. AdV. Polym. Sci. 1994, 116, 207. (28) Insight II, commercial simulation software; Accelrys Inc.: San Diego, CA; www.accelrys.com/InsightII. All simulations in this work were carried out using version 400 P+. (29) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987.

J. Phys. Chem. B, Vol. 112, No. 18, 2008 5659 (30) Mu¨ller-Plathe, F.; Rogers, S. C.; Gunsteren, W. F. Chem. Phys. Lett. 1992, 199, 237. (31) Sok, R. M.; Berendsen, H. J. C.; van Gunsteren, W. F. J. Chem. Phys. 1992, 96, 4699. (32) Pant, P. V. K.; Boyd, R. H. Macromolecules 1993, 26, 679. (33) Mu¨ller-Plathe, F.; Rogers, S. C.; Gunsteren, W. F. J. Chem. Phys. 1993, 98, 9895. (34) Tamai, Y.; Tanaka, H.; Nakanishi, K. Macromolecules 1994, 27, 4498. (35) Gee, R. H.; Boyd, R. H. Polymer 1995, 36, 1435. (36) Chassapis, C. S.; Petrou J. K.; Petropoulos J. H.; Theodorou D. N. Macromolecules 1996, 29, 3615. (37) Fukuda, M. J. Chem. Phys. 1998, 109, 6476. (38) Bharadwaj, R. K.; Boyd, R. H. Polymer 1999, 40, 4229. (39) Boshoff, J. H. D., Lobo, R. F.; Wagner, N. J. Macromolecules 2001, 34, 6107. (40) Karlsson, G. E.; Johansson, T. S.; Gedde, U. W.; Hedenqvist, M. S. Macromolecules 2002, 35, 7453. (41) Budzien, J.; McCoy, J. D.; Adolf, D. B. J. Chem. Phys. 2003, 119, 9269. (42) Alentiev, A.; Economou, I. G.; Finkelshtein, E.; Petrou, J.; Raptis, V. E.; Sanopoulou, M.; Soloviev, S.; Ushakov, N.; Yampolskii, Y. Polymer 2004, 45, 6933. (43) Raptis, V. E.; Economou, I. E.; Theodorou, D. N.; Petrou, J.; Petropoulos, J. H. Macromolecules 2004, 37, 1102. (44) Economou, I. G.; Raptis, V. E.; Melissas, V. S.; Theodorou, D. N.; Petrou, J.; Petropoulos J. H. Fluid Phase Equilibr. 2005, 15, 228. (45) Hu, N.; Fried, J. R. Polymer 2005, 46, 4330. (46) Pavel, D.; Shanks, R. Polymer 2005, 46, 6135. (47) Meunier, M. J. Chem. Phys. 2005, 123, 134906. (48) Hofmann, D.; Ulbrich, J.; Fritsch, D.; Paul, D. Polymer 1996, 37, 4773. (49) Hofmann, D.; Fritz, L.; Ulbrich, J.; Paul, D. Polymer 1997, 38, 6145. (50) Fritz, L.; Hofmann, D. Polymer 1997, 38, 1035. (51) Fritz, L.; Hofmann, D. Polymer 1998, 39, 2531. (52) Hofmann, D.; Fritz, L.; Paul, D. J. Membr. Sci. 1998, 144, 145. (53) Hofmann, D.; Fritz, L.; Ulbrich, J.; Schepers, C.; Bo¨hning, M. Macromol. Theory Simul. 2000, 9, 293. (54) Hofmann, D.; Fritz, L.; Ulbrich, J.; Paul D. Comput. Theor. Polym. Sci. 2000, 10, 419. (55) Heuchel, M.; Hofmann, D. Desalination 2002, 144, 67. (56) Hofmann, D.; Heuchel, M.; Yampolskii, Y.; Khotimskii, V.; Shantarovich, V. Macromolecules 2002, 35, 2129. (57) Hofmann, D.; Entrialgo-Castano, M.; Lerbret, A.; Heuchel, M.; Yampolskii, Y. Macromolecules 2003, 36, 8528. (58) Heuchel, M.; Hofmann, D.; Pullumbi, P. Macromolecules 2004, 37, 201. (59) Fried, J. R.; Ren, P. Comput. Theor. Polym. Sci. 2000, 10, 447. (60) Lo´pez-Gonzale´z, M.; Saiz, E.; Guzma´n, J.; Riande, E. J. Chem. Phys. 2001, 115, 6728. (61) Kucukpinar, E.; Doruker, P. Polymer 2003, 44, 3607. (62) Greenfield, M. L.; Theodorou, D. N. Macromolecules 1993, 26, 5461. (63) Greenfield, M. L.; Theodorou, D. N. Mol. Simul. 1997, 19, 329. (64) Gray-Weale, A. A.; Henchman, R. H.; Gilbert, R. G.; Greenfield, M. L.; Theodorou, D. N. Macromolecules 1997, 30, 7296. (65) Greenfield, M. L.; Theodorou, D. N. Macromolecules 1998, 31, 7068. (66) Greenfield, M. L.; Theodorou, D. N. Macromolecules 2001, 34, 8541. (67) Widom, B. J. Chem. Phys. 1963, 39, 2808. (68) Widom, B. J. Phys. Chem. 1982, 86, 869. (69) June, R.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866. (70) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. Chem. Eng. Sci. 2001, 56, 2789. (71) Tonelli, A. E. J. Polym. Sci., Polym. Phys. Ed. 2002, 40, 1254. (72) Theodorou, D. N. Comput. Phys. Commun. 2005, 169, 82. (73) Foteinopoulou, K.; Karayiannis, N. C.; Mavrantzas, V. G.; Kro¨ger, M. Macromolecules 2006, 39, 4207. (74) Samara, C. Simulation of Polypropylene of Various Tacticities with the Monte Carlo Method. Ph.D. Thesis, Department of Chemical Engineering, University of Patras, December 2000. (75) Doxastakis, M.; Mavrantzas, V. G.; Theodorou, D. N. J. Chem. Phys. 2001, 115, 11339. (76) Gestoso, P.; Nicol, E.; Doxastakis, M.; Theodorou, D. N. Macromolecules 2003, 36, 6925. (77) Wick, C. D.; Theodorou, D. N. Macromolecules 2004, 37, 7026. (78) Wick, C. D.; Siepmann, J. I.; Theodorou, D. N. J. Am. Chem. Soc. 2005, 127, 12338. (79) Zacharopoulos, N.; Vergadou, N.; Theodorou, D. N. J. Chem. Phys. 2005, 122, 244111.

5660 J. Phys. Chem. B, Vol. 112, No. 18, 2008 (80) Mavrantzas, V. G.; Theodorou, D. N. Macromolecules 1998, 31, 6310. (81) Vacatello M.; Avitabile G.; Corradini P.; Tuzi A. J. Chem. Phys. 1980, 73, 548. (82) Karayiannis, N. C.; Daoulas, K. C.; Mavrantzas, V. G. Macromolecules Submitted. (83) Dodd, L. R.; Boone, T. D.; Theodorou, D. N. Mol. Phys. 1993, 78, 961. (84) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (85) Plimpton, S. Large-Scale Atomic/Molecular MassiVely Parallel Simulator (LAMMPS); Sandia National Laboratories: Albuquerque, NM, 2001. All MD simulations in this work were carried out using version LAMMPS 2001 (Fortran 90). (86) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990. (87) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Kein, M. L. Mol. Phys. 1996, 87, 1117. (88) Nose´, S. Mol. Phys. 1984, 52, 255. (89) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (90) Andersen, H. C. J. Comput. Phys. 1083, 52, 24. (91) Karayiannis, N. C.; Mavrantzas, V. G. Macromolecules 2005, 38, 8583. (92) Dee, G. T.; Ougizawa, T.; Walsh, D. J. Polymer 1992, 33, 3462. (93) Cubero, D.; Quirke, N.; Coker, D. F. J. Chem. Phys. 2003, 119, 2669. (94) Sun, H.; Ren, P.; Fried, J. R. Comput. Theor. Polym. Sci. 1998, 8, 229. (95) Sun, H. J. Phys. Chem. B 1998, 102, 7338. (96) Theodorou, D. N.; Suter, U. W. Macromolecules 1985, 18, 1467.

Gestoso and Karayiannis (97) Materials Studio, commercial simulation software; Accelrys Inc.: San Diego, CA; www.accelrys.com/mstudio. All simulations in this work were carried out using version 3.0. (98) Michaels, A. S.; Bixler, H. J. J. Polym. Sci. 1961, 50, 413. (99) Klopffer, M. H.; Flaconneche, B. Oil Gas Sci. Technol. 2001, 56, 223. (100) Teplyakov, V.; Meares, M. Gas Sep. Purif. 1990, 4, 66. (101) Yampolskii, Y.; Wiley, D.; Maher, C. J. Appl. Polym. Sci. 2000, 76, 552. (102) Cuthbert, T. R.; Wagner, J. J.; Paulitis, M. E.; Murgia, G.; D’Aguanno, B. Macromolecules 1999, 32, 5017. (103) Cuthbert, T. R.; Wagner, N. J.; Paulaitis, M. E. Macromolecules 1997, 30, 3058. (104) Charati, S. G.; Stern, S. A. Macromolecules 1998, 31, 5529. (105) Miyake, H.; Matsuyama, M.; Ashida, K.; Watanabe, K. J. Vac. Sci. Technol., A 1983, 1, 1447. (106) Cowling, R.; Park, G. S. J. Membr. Sci. 1979, 5, 199. (107) Sato, Y.; Fujiwara, K.; Takikawa, T.; Sumarno, S.; Takishima, S.; Masuoka, H. Fluid Phase Equilb. 1999, 162, 261. (108) Lim, S. Y.; Tsotsis, T. T.; Sahimi, M. J. Chem. Phys. 2003, 119, 496. (109) Lundberg, J. L.; Wilk, M. B.; Huyett, M. J. J. Polym. Sci. 1962, 57, 275. (110) Mi, Y.; Lu, X.; Zhou J. Macromolecules 2003, 36, 6898. (111) van der Vegt, N. F. A. Macromolecules 2000, 33, 3153. (112) von Solms, N.; Nielsen, J. K.; Hassager, O.; Rubin, A.; Dandekar, A. Y.; Andersen, S. I.; Stenby, E. H. J. Appl. Polym. Sci. 2004, 91, 1476. (113) Sanchez, I. C.; Rodgers, P. A. Pure Appl. Chem. 1990, 62, 2107. (114) Lundberg, J. L. J. Polym. Sci., Part A 1964, 2, 3925.

Molecular Simulation of the Effect of Temperature ... - ACS Publications

Recommend Documents