Detailed Atomistic Simulation of the Barrier Properties of Linear and Short–Chain Branched Polyethylene Melts Through a Hierarchical Modeling Approach

Multiscale Modelling of Polymer Properties M. Laso, E.A. Perpete (Editors) © 2006 Elsevier B.V. All rights reserved.

201

Chapter 2

Detailed Atomistic Simulation of the Barrier Properties of Linear and Short-Chain Branched Polyethylene Melts Through a Hierarchical Modeling Approach Patricia Gestoso,^'^ Nikos Ch. Karayiannis"^'^ ''Accelrys Ltd., 334 Cambridge Science Park, Cambridge CB4 OWN, UK ^Rhodia Recherches, Centre de Recherches de Lyon, Saint-Fons Cedex 69192, France ^Department of Chemical Engineering, University ofPatras, Patras 26504, Greece '^Institute of Chemical Engineering and High Temperature Chemical Processes, Patras 26504, Greece

I. Introduction Sorption and diffusion of small molecules in polymers play an influential role in industrial applications. Many technologically important processes rely on the design of macromolecular materials with tailored barrier characteristics, with the permeability coefficient {P) being the key factor determining the quality of barrier end-products. Among others, gas permeation through polymers is critical to the technology of membranes for industrial gas separation and purification, the production of packaging materials, coatings, biosensors and drug implants and the removal of residual monomer or solvents from the products of polymerization. For example, a polymeric membrane designed for use in food and beverage packaging has to be practically impermeable to certain gases like oxygen (O2). In a similar fashion, storage tanks should efficiently contain gases, and thus exhibit very low permeability coefficients. Finally, in gas separation and purification processes the design goal is the development of a barrier material selective to one molecular species.

202

P. Gestoso andN.Ch. Karayiannis

To successfully confront the challenge of designing polymeric membranes with tailored barrier characteristics, a connection must be established between macromolecular morphology and chemical composition with the barrier properties of the end-product material. Towards this direction, one may resort to direct experimental measurements consisting of two main steps: first, the novel membrane is synthesized and characterized regarding the morphological features, and then its barrier properties are measured and tested under a wide variety of conditions [1,2]. In industrial practice, the material morphology is usually dictated by the manufacturing conditions and the thermodynamic properties of the material throughout its processing history. While experimental studies represent undoubtedly the most reliable and direct way to determine the ability of a novel polymer to act as an effective barrier structure, the execution of a large series of trial-and-error synthetic processes and successive permeability measurements may be proven, in practice, too expensive and/or time-consuming. Li parallel to direct calculations, large collections of experimental data over a wide range of systems and conditions may serve as the basis for the fabrication of phenomenological correlations [3]. From the theoretical point of view, sorption and diffusion in and through macromolecular systems can be explained and predicted based mainly on the concept of free volume (FV) [4-6], the dual-mode transport model [7-10] and molecular theories [11-14]. In practice, no plastic material is perfectly impermeable: Sorption and diffusion of light gases is possible, to certain extend, in any polymeric matrix as the packing of chains and their constituent atoms, even in the case of very dense polymers, leaves holes or micro voids [15], the sum of which forms the free (or unoccupied) volume (FV). The cavities of free volume can accommodate gas molecules depending on the size and shape of the penetrants relative to their own volume and shape, and the polymer-penetrant interactions in relation to the penetrant-penetrant and polymer-polymer interactions. The subset of sorption sites where the molecules can reside comprises the accessible volume (AV) of the polymer, depending on both the host matrix and the penetrant. Finally, computer simulations and modeling can play an important role in the establishment of the dependence of barrier properties on structure by elucidating, at the atomistic level, all the relevant mechanisms controlling the whole transport process. Computer-generated polymer structures built on detailed molecular models can be subjected to novel simulation techniques that allow, within modest computational time and resources, the accurate prediction of permeability under a wide range of conditions without the usually high cost of experiments or without resorting to simplifications and assumptions invoked by most theoretical studies. In recent years, the continuous increase in computational power, the affordable cost of powerful workstations and the

203

Atomistic simulation of the barrier properties

development and implementation of novel simulation techniques have drawn the attention of the polymer community, setting computer-aided design of materials in par with the well-established experimental and theoretical approaches. Li the present chapter we report our recent work concerning the detailed atomistic simulation, through a hierarchical modeling approach, of the sorption and diffusion of small gas penetrant molecules at infinite dilution (low concentration) through purely amorphous polyethylene (PE) systems characterized by different molecular architectures and molecular weights under a variety of temperature conditions. Section II presents a brief description of the molecular mechanisms and the physical aspects of gas permeation in polymers. Section III deals with the description of the methods adopted in the simulation of sorption and diffusion in rubber and glassy polymers. Results of our recent work are presented in Section IV. Finally, main conclusions and future plans are analyzed in Section V. II. Molecular mechanism of low-concentration gas transport in polymers A coarse schematic representation of the low-concentration gas permeation in a polymeric membrane is given in Fig. 1.

y |iiiiiiiiiiiiB

-<

X

Po

•

Figure 1. Coarse schematic representation of the low-concentration permeation of penetrant molecule in polymeric membrane of thickness / as a result of the pressure difference Ap.

P. Gestoso andN.Ch. Karayiannis

204

Permeation can be envisioned as a three-step process where the penetrant molecule a) is absorbed at the surface of the inlet side of the membrane, b) diffuses through the polymer matrix and c) is desorbed for the outlet surface of the membrane [16,17]. At equilibrium, penetrant concentration, c, is related to the pressure, /?, by the isothermal relation c = S{c)p

(1)

where S{c) is the solubility coefficient. At infinite dilution (low concentration of the penetrant molecule) Eq. 1 reduces to a form of Henry's law and solubility is practically independent of concentration (or pressure). Accordingly, penetrant concentration in the inlet and outlet surfaces of the polymeric membrane (Fig. 1) is given by ^in=*^An '^out=^/^out

(2)

Transport process across the membrane can be quantified by the permeation rate (flux), y, according to Pick's first law J =-Z)^ dx

(3)

where D is the diffusion coefficient. The rate of change in the concentration profile is given by

^ =± Z , ^ d/ dx V doc)

(4)

If the diffusion coefficient, D, is independent of the position in the sample, x, (i.e. the medium is homogenous [18]), then Eq. 4 corresponds to Pick's second law

At

dx^

Atomistic simulation of the barrier properties

205

At steady state, Eq. 5 is reduced to ... (6)

dc ^ d^c ^ dc ^ ^ — = 0 => —r- = 0 => — = constant

dt

dx

dx

Based on Eq. 2 and the steady-state conditions in Eq. 6, the permeation rate (Eq. 3) can be rewritten as

/

I

I

The low-concentration permeabiUty coefficient, P, is given as the product of diffusion and solubility coefficients P =D5

(8)

where the permeability coefficient, P, is related to permeation flux, J, through

P='-^ JAp

(9)

Based on Eq. 8, the low-concentration permeation of a penetrant molecule in a polymer matrix, expressed through the permeability coefficient (P), can be considered as the combined effect of kinetic (diffusion) and thermodynamic (sorption) factors, quantified by the diffusion (Z)) and solubility (5) coefficients, respectively. Solubility of a penetrant in a polymeric membrane depends mainly on the nature and magnitude of the polymer-penetrant interactions relative to the polymer-polymer and penetrant-penetrant interactions. In addition, clusters of free volume (sorption sites) should exist within the macromolecule, large enough to accommodate the penetrant molecules [19]. Penetrant diffusion is mainly controlled by the size and shape of the penetrant, the magnitude of the potential interactions with the host matrix and the size, shape, distribution and connectivity of the free-volume cavities. In the melt state, at high temperatures (well above the glass transition temperature Tg), polymer segments undergo significant thermal fluctuations leading to enhanced

P. Gestoso andN.Ch. Karayiannis

206

chain mobility. Consequently, the size, shape, distribution and connectivity of the sorption sites are continuously rearranged and the penetrant is "carried along" by density fluctuations caused by the thermal motion of the surrounding polymer chains [19]. In a glassy polymer {T < Tg) the corresponding transport mechanism is considerably different: glassy configurations are trapped in local minima of the potential energy and conformational transitions between different energy minima are prohibited by exceptionally high energy barriers. As the mobihty of the polymer segments is significantly reduced, the distribution and connectivity of the network of sorption sites remains practically unaltered within the time frame of the transport process, undergoing only minor fluctuations. The work of Takeuchi [20] revealed that the penetrant molecule spends most of its time trapped within a formed cavity of free volume and only infrequently it performs jumps from one cluster to a neighboring one through channels which are formed instantaneously via fluctuations in regions of lower density or enhanced molecular mobility. Consequently, penetrant diffusion in glassy polymers is orders of magnitude slower than in high-temperature melts. III. Simulation methods to study permeation in polymers The most straightforward way to study the transport behavior of a small molecular weight penetrants in a polymer melt is to employ Molecular Dynamics (MD) simulations [21] (see also chapter II.6 of the book), following the motion of the penetrant for long times when the hydrodynamic limit is reached and Fickian diffusion is established. Diffusion coefficient, Z), is readily calculated by invoking the Einstein equation

D = \mi\

([rp(O-r/0tf)] 6/

(10)

where the brackets, <>, denote averaging over all trajectories and [rp(0-rp(0)]^ is the mean square displacement (MSD) of the penetrant at time t. In the most common implementation, MD simulations track the dynamics and motion of both the polymer matrix and the penetrant molecules, which either coexist from the very beginning of the simulation (i.e. the chain segments are originally grown and relaxed in the presence of penetrants) or the penetrants are inserted a posteriori in the amorphous cell in an energetically-biased way to avoid overlaps with the existing polymer segments [22-38].

Atomistic simulation of the barrier properties

207

The main advantage of the MD method is that, by providing direct and detailed information about system dynamics, it practically constitutes a computational experiment, as the assumptions and simplifications invoked in its implementation are only related to the applied potential force field (molecular model). In parallel, MD's major limitation emanates as well from its working pattern: even in a state-of-the-art implementation executed on parallel supercomputers, thousands of CPU hours are required to study the real dynamics and transport behavior of large penetrant molecules in complex macromolecular environments. Even worse, at temperatures below Tg, penetrant mobility in glassy structures is too slow to be tracked by conventional MD simulations within reasonable computational time. An excellent alternative is the Transition State Theory (TST) as introduced by Arrizi, Gusev and Suter [39-42] for the computational study of gas transport in polymeric glasses. Li the limit of low concentration, TST is applied on detailed atomistic systems and evolves in a coarser level of representation as a succession of infrequent events. For each transition, the "reaction trajectory" leading from a local energy minimum to another through a saddle point in configuration space is tracked, and the transition rate constant is evaluated. Li the Gusev-Suter implementation [40-42], TST is initiated as a 3-D grid of fine resolution (typically 2-3 A) is laid on the amorphous cell consisting of the equilibrated polymer chains, covering every part of its volume. Next, a spherical probe with dimension identical to the one of the penetrant molecule (represented as a united-atom site) is inserted in every point (x, y, z) of the 3-D grid and the potential energy £'ins(jc, y, z) is calculated as the non-bonded interactions between the probe and the polymer segments. Thus, based on the calculated value of ^ms(^, y, z), the simulation cell is divided into clusters of accessible volume (regions of low energy, containing the local minima of the energy) and domains characterized by high energy (excluded volume regimes). If there are no potential interactions between the dissolved species (i.e. infinite dilution) the excess chemical potential, //ex, according to the Widom method, is given by

//„=/?nn

exp

^ E. ^ ms

(11)

where R is the universal gas constant, T is the temperature, k^ the Bohzmann constant and the brackets, <>, denote averaging over all the grid points and probe insertions. Additionally, the low-concentration solubility coefficient, S, is related to the excess chemical potential, /iex, through

P. Gestoso andN.Ch. Karayiannis

208

'5 = exp(-|^J

(12)

Each energy local minimum is associated with a sorption "microstate" in the configuration space of the penetrant. Adjacent microstates (i.e. sorption sites) are separated by high-energy surfaces. The penetrant molecule spends most of its time "trapped" in the microstates and only infi*equently hops from one microstate to a neighboring one. Penetrant diffusion can be envisioned as a succession of rare transitions between adjacent microstates. An elementary transition from microstate / to microstate 7 can be described by a characteristic first-order rate constant, /:,_;. For temperatures low enough, for which the polymer matrix can be considered frozen, the rate constant for the transition from microstate / to microstatej is given by [39-42]

* , . , = ^ ^


where h is Plank's constant, k^ is the Boltzmann constant and g/, Qy are the partition fimctions of the molecule in the microstate / and on the high-energy surface (separating microstates / and7), respectively. The partitionfimction,Qy on the separating surface is given by [39^2]

Qij =

7 1 ^ — \ p(x,y,z)dA

(14)

where m is the mass of the dissolved molecule, p(jc, y, z) is the Boltzmann probability density function for finding the inserted atom at point (x, y, z) and dL4 is an elementary surface area on the separating surface between / and 7, whose total area is denoted by Ay. Similarly, the partition function of the molecule in microstate / is given by [39-42]

Q,=[^^T

lp(x,y,z)jy

(15)

Atomistic simulation of the barrier properties

209

where dF is an elementary volume of the microstate i whose total volume is denoted by Vt. Finally, the probability p{x, y, z) of finding the inserted atom at point (jc, y, z) is given by

p(x,y,z)ocQxp

Einsix.y.z)^

(16)

kj

Li the original formulation of TST for glassy systems the polymer matrix was considered as completely frozen, leading to significant underestimation of the solubility and diffusivity coefficients. In the modified approach, the small 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. These displacements around the equilibrium positions follow a Gaussian probability density function [39-42]

-E^2(Ar,^) VY

fr(Ar,, Ar2...., Ar^) oc exp

\ (17)

where N is the total number of polymer atoms in the simulation box and Ar^t = fk - is the displacement of atom k from its equilibrium position, . In the simplest version, an average value independent of atom type, 3A^, is assigned to the corresponding MSD , A^ being termed as "smearing factor". In the modified version of TST that incorporates the elastic motion of the polymer segments the corresponding density function ^(x,;^, z) is given by

p(x.y.z)«Jexpl-5M|^

A

^(Ar„...,Ar^)c/V,...jV^

(18)

where W(Ari,...,ArAf) is given by Eq. 17. Knowledge of the probability density function allows the calculation of the transition rate constants through Eqs. ISIS. At equilibrium, the average residence time, r, in a microstate /, is given by [19,39-42,45]

P. Gestoso andN.Ch.

210

Karayiannis

/

where the index / runs over all sorption sites adjacent to the reference microstate /. Consequently, the corresponding probability for the transition from microstate / to the adjacent microstate7,/7i_j, is given by

"-''it.

<^"'

Based on the information about the distributions of rate constants and lengths governing the transitions, the spatial connectivity and topology of the microstates and the mobility of the dissolved molecule (as provided from Eqs. 13,19,20), a kinetic Monte Carlo (kMC) simulation is performed to generate long stochastic trajectories, each consisting of a large succession of transitions of penetrant walkers hopping in the 3-D coarse-grained lattice [19,39-42,45,46]. The low-concentration diffusion coefficient, D, can be calculated from the Einstein equation (Eq. 10) when Fickian diffiision is established. Early TST simulations [39-42] demonstrated that for small times scales the motion of the penetrant is not diffusive. Listead, an "anomalous diffiision" behavior is observed, where the mean square displacement of the penetrant, [rp(0-rp(0)]^. rp(0-rp( scales with time, t, as [Y^{tyr^{0)f ©c f with « < 1. Similar conclusions were drawn by the early MD works of Miiller-Plathe et al. for oxygen (O2) and helium (He) diffusion in amorphous polyisobutylene (PIB) [22], Pant and Boyd for methane (CH4) in polyethylene (PE) [24] and Chassapis et al. for O2 in PE [28]. Experimentally, short-time anomalous diffiision behavior has been identified through NMR measurements during diffiision of penetrants in zeolite systems by Karger et al. [47] and in molecular sieves by Ylihautala et al. [48]. An interesting observation in the original Gusev-Suter TST simulations [41] is that the crossover from anomalous to Fickian diffiision occurs when the MSD of the penetrant walkers is commensurate to the dimensions of the amorphous cell, suggesting that Fickian diffusion is artificially established in the simulation by the applied periodic boundary conditions. This finding has been confirmed by a very recent TST approach for O2 diffiision in glassy configurations of poly(ethylene terephthalate) (PET) and poly(ethylene isophthalate) (PEI) [49]. To capture correctly the duration of the anomalous behavior one should repeat TST simulation on simulation boxes with progressively larger dimensions. On

Atomistic simulation of the barrier properties

211

the other hand, as indicated in the work by Karayiannis et al. [45], the diffusion coefficient of penetrant walkers, for a given lattice with prescribed topology and connectivity, depends solely on the distribution of rate constants governing the transitions between microstates. Short-range spatial heterogeneity affects the duration of the anomalous behavior but has no effect on the estimated diffusivity. The rather coarse representation of the thermal motion of the polymer atoms as expressed by the smearing factor (same for all atomic species) limits the apphcability of TST [39-42] to relatively simple and small gas molecules (He, H2, O2, N2 and marginally CH4), whose sorption in the polymer matrix does not affect the thermal motion of the surrounding chains. Large molecules, for example CO2, may force the polymer atoms to undergo significant local relaxation in order to accommodate their presence, rendering the smearing factor of Eq. 17 obsolete to effectively capture the local dynamics of the chains. In addition, the simplified harmonic form of the smearing factor in the original TST approach is not adequate to describe more complex segmental relaxation mechanisms like phenyl ring vibrations and flips that may have a profound impact on the "dynamic flexibility" of the system [49]. Solution to this central shortcoming of the original TST can be provided by allowing the local relaxation of the polymer environment around the dissolved molecule through appropriate implementation of the multidimensional Transition State Theory, as proposed by Greenfield and Theodorou [50-54]. To this date, TST has been successfully applied on various penetrantpolymer systems providing reliable and accurate predictions for the solubility, (5) diffusivity (D) and permeability (P) coefficients (the latter being indirectly calculated from Eq. 8) of light gases (He, H2, O2, N2 and CH4) in glassy polymers. For larger and non-spherical penetrants, like CO2, TST fails to capture the correct sorption and transport behavior [40-42,49,55-68]. As the TST approach departs from the detailed atomistic level and resorts to a coarser level of description and because of the invoked simplifications related with the simplistic representation of the segmental relaxation of the polymer host, the approach should be validated against available data either from experimental studies or MD simulations. Critical preconditions of utmost importance for the successful prediction of the barrier properties, using either TST or MD, are a) the generation of realistic atomistic configurations of the polymer systems, fully equilibrated at all length scales, b) the calculation of the barrier properties as statistical averages over as many as possible totally uncorrelated computer-generated structures and c) the apphcation of a very accurate and, if possible, CPU inexpensive force field. As analyzed extensively in chapter 1.2 of the book, for chemically simple macromolecular systems (for example polyolefms, polydiens and polyethers)

P, Gestoso and N. Ch. Karayiannis

212

full-scale equilibration is provided by the application of stochastic Monte Carlo algorithnis built around advanced chain-connectivity altering moves [69-74]. In the case of more complex macromolecular systems bearing aromatic units and bulky inflexible and polar groups along the backbone (for example polyesters and polyimides), equilibration is achieved by a procedure consisting of a large succession of short MD runs in the canonical (NVT) or isothermal-isobaric (NPT) ensembles, corresponding practically to a series of annealings, coolings, compressions and decompressions at various conditions of temperature and pressure [55-65]. hi the next section we will present results about the low-concentration solubility, diffusivity and permeability of oxygen and nitrogen in linear and short-branched polyethylene melts as obtained from application of a hierarchical modeling methodology integrating simulation techniques applied on different time and length scales. IV. Multiscale Modeling of the Barrier Properties of Polyethylene Systems IV. 1. Summary of the Hierarchical Approach The adopted hierarchical methodology is divided into three phases: 1. Application of chain-connectivity altering MC moves [69-74] for the rapid and robust equilibration of purely amorphous PE samples in the atomistic level of description employing a united-atom representation. 2. Conversion of selected representative and uncorrected configurations of the simulated systems to an explicit (all-atom) representation by appropriately adding the corresponding hydrogen atoms. Application of short NVT MD runs accompanied by static energy minimizations to afford additional local relaxation due to the addition of hydrogen atoms. 3. AppHcation of the Gusev-Suter Transition State Theory (TST) [3942] to calculate the barrier properties of the computer-generated PE samples at conditions of infinite dilution (low concentration). IV.2. Systems Studied All simulated linear polyethylene samples are denoted as "C^r" where N is the total number of carbon atoms per chain. The short-chain branched (SCB) analogs are referred to as "SCB_(A^br+l)^A^freq_^br^Cb", where A^br is the number of branches distributed regularly along the main (linear backbone), Cb is the

Atomistic simulation of the barrier properties

213

branch length (i.e. the number of carbon atoms of each branch) and TVfreq is the branching frequency (i.e. the number of carbon atoms along the main backbone between two successive branch points). More details about the notation and the molecular characteristics of the SCB structures can be found in chapter II.6. Regarding the linear PE systems, at all studied temperatures a small degree of polydispersity was allowed to accelerate the equilibration phase performed through MC techniques, which is the most time-consuming stage of the whole scheme. For the short-chain branched systems all molecular characteristics were kept strictly monodisperse (branch length, Cb, branch frequency Nf^^, and number of branches, iVbr) to avoid adding extra structural parameters, such as heterogeneity in the branching frequency and in the number of branches per chain. In addition, all SCB samples bear exactly the same number of carbon atoms (corresponding to a total molecular weight of 1990g/mol, i.e. linear C142 PE as an analog) so as to attribute any detectable differences in the barrier properties solely on molecular characteristics. Table I summarizes the details for all simulated PE systems. Initial linear structures were provided by the three-stage constant-density energy minimization technique of Theodorou and Suter [75,76], which generates configurations trapped at a local energy minimum with reaUstic bonded geometry (bond length, bending angles and torsions). Initial dimensions of the cubic amorphous cell were selected so as to correspond to density values that followed the quantitative trend as a function of molecular length reported from our previous atomistic MC simulations on similar PE structures [71]. The initial configuration for each SCB system was generated through the commercial software package Materials Studio (version 3.0) [77] by Accelrys Inc at identical density values as their linear counterparts of the same molecular weight. All simulations of the equilibration phase were executed with a molecular model (potential force field) adopting the united-atom (UA) representation where each CH3, CH2 and CH unit is treated as a single spherical interacting site (pseudoatom). In previous atomistic simulations the molecular model, a hybrid combination of accurate force fields existing in the literature [78-80], has provided excellent estimates of the volumetric [71], conformational [72] and dynamic [81] properties of PE samples. More details about the employed potential force field regarding the mathematical expressions and the corresponding parameters describing all bonded and non-bonded interactions can be found in chapter 11.6 of the book.

P. Gestoso andN.Ch.

214

Karayiannis

Table I. Details of the Simulated PE Systems.

no. of chains 40 C78 22 Cl42 6 C500 3 ClOOO SCB 11 xi2; 10x1 22 22 SCB 9x 14 8x2 22 SCB 7x 16 6x5 22 SCB 5x:22 4x8 System

no. of sites 3120 3124 3000 3000 3124 3124 3124 3124

Total MW (g/mol) 1094 1990 7002 14002 1990 1990 1990 1990

Polydispersity Index 1.083 1.083/1.000 1.053 1.053 1.000 1.000 1.000 1.000

IV. 3, Step 1: Equilibration of samples in united-atom representation All linear monodisperse or polydisperse PE systems were subjected to exhaustive MC simulations based on chain-connectivity altering moves. As described in chapter 1.2 in the case of polydisperse samples, molecular weight distribution is controlled by casting the MC simulations in the semi-grand statistical ensemble through the use of an appropriately defined spectrum of chemical potentials [69,70]. Chain lengths of polydisperse systems were sampled uniformly in the closed interval [iVav(l - A), A^av(l + A)], where TVav is the number-average chain length and A denotes the half-width of the chainlength distribution reduced by A^av The polydispersity index, /, (reported in Table I for all PE samples) is related to A through / = 1 + A^ / 3. The simulation of polydisperse linear systems was carried out with the following mix of moves: internal libration (flip) (5%) [82], end-mer rotation (3%), intermolecular reptation (15%) [69,73], concerted rotation (20%) [83], end bridging (15%) [69,70], double bridging (20%) [71-73], intramolecular double rebridging (10%) [71-73], intramolecular end bridging (10%) and volume fluctuation (2%). To compare against the SCB analogs, a single strictly monodisperse C142 PE sample was simulated using a reduced set of moves (all algorithms that introduce polydispersity were excluded): internal libration (5%), end-mer rotation (3%), reptation (15%), concerted rotation (20%)), double bridging (35%), intramolecular double rebridging (10%), intramolecular endbridging (20%) and volumefluctuation(2%). The complex molecular architecture of SCB systems required the application of a more advanced MC scheme consisting of novel local or chainconnectivity altering moves [73]. Since all branching parameters were held fixed, direct application of chain-connectivity altering algorithms was not possible at this stage. For this reason all initial SCB configurations were

215

Atomistic simulation of the barrier properties

subjected to very long NPT MD simulations for a total duration of 100 and 400ns at T = 450 and 300K, respectively. All MD simulations were conducted using the LAMMPS (large-scale atomic/molecular massively parallel simulator) software (version LAMMPS 2001) [84,85]. Temperature and pressure fluctuations were controlled through the Nose-Hoover thermostat [86,87] and Andersen barostat [88], respectively. For the integration of the equations of motion, the multiple time step algorithm rRESPA [89,90] was employed with the small and large steps set equal to 1 and 5fs, respectively. More details about the dynamics of the SCB PE melt configurations can be found in chapter II.6 of the book.

5000

T"

instantaneous value -running average value selected frames

4500 4000 3500 3000 V 25002000 1500 1000

-»—r 1

10 MC Steps (10')

Figure 2. Instantaneous and running average values of the end-to-end distance, , as a function of MC steps for polydisperse, linear C142 PE system at 7 = 450K and P = latm. Filled cycles denote representative and statistically uncorrelated frames from the MC trajectory taken for successive TST calculations.

P. Gestoso andN.Ch. Karayiannis

216

instantaneous value -running average value selected frames

0.79 H

0.78 H

0.74

I

0

1

I

I

2

I

I

3

I

4

I

I

5

I

I

I

6

I

7

I

I

8

I

I

9

I

I

10

MC steps (10')

Figure 3. Same as in Fig. 2 but for the density.

At the end of the MC (or MD) simulations, a set of representative and statistically uncorrelated frames was carefully selected from the whole trajectory with dimensions (the radius of gyration for SCB, and the end-to-end distance for linear melts) and density similar to the average values obtained from sampling over the equilibrated part of the trajectory. Figs 2 and 3 present the instantaneous and running average values of the mean squared end-to-end distance, , and density, respectively, as a function of MC steps as obtained from simulations on linear polydisperse C142 system at r = 450K and P = latm. In both figures filled cycles denote representative and statistically uncorrelated frames from the MC trajectory selected for successive TST calculations. This procedure was repeated over all MC trajectories and the number of selected frames increased as the temperature decreased. Similar procedure for extracting configurations for successive studies of the barrier properties was adopted in the case of the NPT simulations of SCB analogs. Figure 4 presents the dependence of density on average chain length, A/av (or equivalently molecular weight) as obtained from MC simulations on linear polydisperse PE samples at T = 300 and 450K. Also shown are available experimental volumetric data at the same conditions of temperature and pressure as the simulations [91,92]. At room temperature {T = 300K) PE is a semi-crystalline material and the density value reported by Cubero et al. [92] is corrected so as to correspond to a purely amorphous sample. It is evident that at both temperatures the agreement between the density as calculated from MC

Atomistic simulation of the barrier properties

211

simulations and as reported experimentally is excellent, the relative error being less than 1%. As reported in a previous work [49], the precise calculation of density is a critical prerequisite for the successful and accurate prediction of the barrier properties. Even for relatively small errors of about 2-3%, one may need to impose the density to match the corresponding experimentally known values; otherwise consecutive TST simulations may fail to capture the qualitative and quantitative trends for sorption and diffusion [49]. Fig. 5 depicts the variation of density with temperature as obtained from MC simulations on a polydisperse, linear Cyg system. As shown by the linear fit in Fig. 5 it can be concluded that density increases hnearly with temperature in the whole range of appUed temperatures. Divergence from linearity at the lowest studied temperature (T= 250K) can be attributed to the onset of transition to the glassy state. Recent MC simulations suggest a glass transition temperature in the range of 190-200K for the linear C78 amorphous sample. At high enough temperatures, in the melt regime, MC is capable of equilibrating the volumetric and conformational characteristics of PE systems within two CPU hours. 0.88-1

1

1

1

1

1 —— 1 _

1

1

—1

1

1

0.86-

•

0.84-^ 0.82-

• \

£

3

A

0.80 -

g 0.78-

i ]

a> 0

0.760

0.74n 77 -

(3

• 0

simulation 300K simulation 450K 1

200

•

1

400

'

D

#

experimental 300K experimental 450K

600

1

800

]

-

1000

Molecular Length

Figure 4. Dependence of density on average molecular length as obtainedfromMC simulations on linear polydisperse PE systems at r = 300 and 450K (in both cases P = latm).

218

P. Gestoso andN.Ch.

0.62 200

250

300

350

400

450

500

550

Karayiannis

600

Temperature (K) Figure 5. Density as a function of temperature as obtainedfromMC simulations on a linear polydisperse C78 PE system at P = latm. The solid line corresponds to the best linear fit on simulation data.

The effect of molecular architecture (short-chain branching) on density is depicted in Fig. 6 as obtained from NPT MD simulations on PE systems bearing a total of 142 carbon atoms at T = 300K and P = latm. Highest density is recorded for the SCB configuration bearing the shortest (and most frequently distributed) branches (Ci), whereas for longer branches (C2-C6) density slightly decreases with increasing branch length. It can be safely concluded that shortchain branching has a very minor effect on density, at least in the range of molecular weight and temperature reported here: The relative difference between the higher and lower density values lies within the statistical error of the simulation method, about 1%. Volumetric trends of Fig. 6 are in excellent qualitative agreement with available experimental data [93]. The density values corresponding to higher temperatures {T = 450K) for molten PE either in a purely linear pattern or bearing symmetrically attached short branches along the main backbone are presented in chapter IL6 (see for example Fig. 7 and the corresponding discussion).

219

Atomistic simulation of the barrier properties

0.88

0.86 H E

c

o 0.84 H

0.82 C142

SBC 11x12 10x1

SBC 9x14 8x2

SBC_7x16_6x5

SBC 5x22 4x8

Molecular Architecture Figure 6. Dependence of density on molecular architecture as obtainedfromNPT MD simulations on short-chain branched (SCB) and linear PE systems at 7 = 300K and P = latm.

IV.4. Step 2: Conversion to explicit-atom configurations - Local relaxation Early MD simulation studies recognized the important role of the detailed representation of the constituent polymer atoms in the calculation of sorption and diffusion coefficients of penetrant molecules. According to the work of MuUer-Plathe et al. [22,94], united-atom (UA) description artificially raises transport rates by almost two orders of magnitude, even if the adopted UA model provides the correct density. An explicit-atom (EA) model packs more efficiently polymer segments of the same configuration, thus leading to a notably different distribution of free volume even if the density of the UA and EA systems is identical [17,22,94]. In the first stage of the hierarchical approach we adopted the UA description for the full-scale equilibration of the simulated samples. Besides the wellestablished asset of providing excellent predictions regarding density (Fig. 4) and intra- and intermolecular packing [71,72], the incorporation of the present UA model resulted in a tremendous saving in computational time and resources for system equilibration (a factor of three when compared against EA models), which is the time-limiting stage of the whole methodology. The critical

220

P. Gestoso andN.Ch. Karayiannis

requirement of a detailed (explicit-atom) representation of the polymer matrix in the final phase of the TST appUcation requires the conversion of the exiting UA configurations into EA structures by appropriately adding a posteriori the corresponding hydrogen atoms [95]. The resuhing EA structures were fiirther adapted to be compatible with the supported file formats of the commercial software Insightll (version 400P+) from Accelrys Inc. [96] employed for the successive TST calculations. The inclusion of hydrogen atoms and the resultant transformation from UA to EA configurations may cause energetic overlaps between added hydrogens and existing polymer segments, leading to unrealistic local packing. Therefore, a supplementary cycle of local equilibration was adopted to remove all energetically unfavorable overlaps and to impose a realistic distribution of free volume along the chains. The protocol used for the additional relaxation consisted of static structure optimizations (using an energy minimization Molecular Mechanics (MM) algorithm [75,76]) accompanied by short NVT MD simulations. The duration of each MD run was set at 50ps, the time integration step was equal to Ifs and the thermostat by Berendsen et al. [97] was employed to dump all temperature fluctuations. Casting of the short MD runs in the canonical (NVT) ensemble was decided on the basis that all simulated PE systems, as obtained from the application of the UA model during the equilibration stage, were characterized by densities in excellent agreement with available experimental data. In all simulations of EA samples all bonded and non-bonded potential interactions were described by the COMPASS [98,99] (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies) force field which is widely considered one of the most accurate, allpurpose EA potentials. Preliminary executions of the proposed hierarchical scheme on selected linear PE configurations using either COMPASS or PCFF [100,101] (Polymer-Consistent Force Field) suggested that the former provided predictions for the diffusion and solubility coefficients of oxygen and nitrogen molecules in better agreement with available experimental data. In addition, preUminary simulations indicated that a total of 3 successive combinations of MM and NVT MD runs were adequate to provide local relaxation as determined from the evohition of the total potential energy of the system and of the inter(g(r)) and intramolecular (w(r)) carbon-carbon distribution fiinctions. After the completion of the third MM/MD cycle, the potential energy was seen to reach a constant value, whereas no appreciable difference was detected in the structural features in the atomic level as evidenced from the pair distribution functions. IV, 5. Step 3: Application of Transition State Theory Gusev-Suter TST was employed, as implemented in the gsnet and gsdif subroutines of the Insightll software [96], to calculate the barrier properties of

221

Atomistic simulation of the barrier properties

PE matrices equilibrated in the first two phases of the present methodology. gsnet calculates the free energy surface of the system, thus separates the volume of the amorphous cell into sorption sites (characterized by low energy) and regimes of high-energy (excluded volume). Additionally, the connectivity of microstates is determined as well as the distribution of the first-order rate constants governing the transitions between microstates (sorption sites). Initially, a 3-D grid was laid on the amorphous cell covering the entire volume. A limitation in the maximum number of points in the 3-D grid in conjunction with the large size of the equilibrated PE configurations ranging from 40 up to 45A, notably larger than the cell dimensions encountered in typical a all-atom TST (25-30A), restricted the grid resolution to 0.5A, which is somewhat coarser than the typical value of 0.3 used in previous TST works [49]. Extensive test simulations with different grid resolutions ensured that the grid spacing of 0.5A did not affect the successive prediction of solubility, diffusivity and permeability. The penetrant molecule was represented as a united-atom, spherical site which interacts with the atoms of the polymer matrix through short-range forces characterized by a 9-6 potential (t/vdw) of the form 9

a vdW (ry) = £ 2

(

-3 G

\6"] (21)

where ry is the distance between the penetrant / and the polymer atom 7 and e and a are the well depth and the collision diameter of the COMPASS force field. For various penetrant molecules the corresponding 6: and a parameters are reported in Table II.

Table II. Well-depth, e, and collision diameter, a, for various gas penetrant molecules as implemented in Gusev-Suter TST. Penetrant

e (kcal/mol)

^(A)

O2 N2 CO2

0.2344 0.1889 0.4500

3.46 3.70 4.00

222

P. Gestoso andN.Ch. Karayiannis

The cutoff distance for the non-bonded interactions was equal to 9.5A, as test runs with larger values revealed no dependence of the calculated barrier properties on the truncation radius. The smearing factor (Eq. 17) was calculated through a self-consistent scheme with the two relevant parameters being the most probable residence time, tp, of the penetrant in the microstates and the smearing factor. Liitially, an arbitrary value was assigned for the smearing factor and the distribution of the residence times in microstates was computed using Eq. 19. Knowledge of the distribution allowed the identification of Tp and a new value was attributed to the smearing factor based on the mean square displacement (MSD) of the polymer atoms from their equilibrium positions for time equal to Tp. To extract the corresponding MSD, an additional 50ps-long NVT MD simulation was executed after the completion of the relaxation cycle. Convergence occurred when the relative error between two successive values of the smearing factor was smaller than a prescribed threshold equal to 2.5%. Based on the adopted self-consistent scheme, the calculated smearing factor depends not only on the segmental mobility of the polymer host and the applied conditions (temperature) but also on the penetrant (as revealed by the inclusion of the most probable residence time). By determining the free energy surface of the system, identifying of all sorption sites and calculating the corresponding smearing factor (a single common value for all different species of the polymer host) through the gsnet subroutine, the low-concentration solubility coefficient, iS, can be readily calculated through Eqs. 11-12. Based on the information about the connectivity and distribution of microstates, the rate constant distribution and the smearing factor, the lowconcentration diffusion coefficient, D, was extracted by employing the gsdif subroutine of Insightll. Kinetic Monte Carlo (kMC) simulations, applied on the coarse network of microstates, provided the trajectories of the penetrant walkers as a function of time. The diffusion coefficient, Z), was calculated in the longtime limit where Fickian diffusion is established through the Einstein equation (Eq. 10). The total duration of the kMC was 10"*s and the mean square displacement (MSD) was averaged over the independent trajectories of 1000 walkers. Finally, the low-concentration permeability coefficient, P, was calculated indirectly as the product DS (Eq. 8). All reported results about the barrier properties (S, D and P) of PE to O2 and N2 came as statistical averages over either four {T> 500K), five (500K >T> 350K) or eight (350K >T> 250K) uncorrelated structures for each simulated system.

Atomistic simulation of the barrier properties

223

IV. 6. TST-based results IV.6.1. Comparison with experimental data A series of TST calculations were conducted at T = 300K and P = latm to validate the accuracy of the proposed multiscale scheme and its applicability on three different gas penetrant molecules: oxygen (O2), nitrogen (N2) and carbon dioxide (CO2), wilh the latter being widely considered [62,63,65,68] as too large (and anisotropic in shape) to be effectively modeled by the present implementation of TST. As akeady mentioned, at standard conditions polyethylene is a semi-crystalline polymer, with a non negligible crystalline volume fraction, v^, especially in the case of linear (HDPE) polyethylene samples. On the other hand, computer-generated PE macrostructures are completely amorphous. To compare TST-based results with available experimental data one should "correct" the latter so as to correspond to ideal purely amorphous samples. Pant and Boyd [24] proposed the following relationship between the experimentally measured diffusion coefficient, Dexp, and the "corrected" one, Z)'exp characterizing non-crystalline PE systems

D

3 D =-—^^ "' 2 ( l - v ^ )

(22)

Similar correction can be applied on the solubility coefficient

<^^>

^•-F^

where iSexp is the experimentally measured solubility coefficient in the semicrystalline material and ^yexp is the "corrected" value corresponding to an ideal non-crystalline sample. Consequently, the permeability coefficient "corrected" for crystallinity, Fexp is given by . exp

.

e.pe.p

2D

S

3(l_vJ(l-v)

2 P 3(l-v)^

^ ^

224

P. Gestoso andN.Ch. Karayiannis

where Pexp is the experimental value of the permeability coefficient of the actual semi-crystalline PE material. All reported experimental data were corrected for crystaUinity (through Eqs. 22-24) so as to correspond to purely amorphous samples allowing the direct comparison with TST-based predictions. Table III lists the low-concentration S, D and P coefficients as obtained from application of the hierarchical scheme to the polydisperse, linear Ciooo system {T = 300K, P = latm) and the corresponding experimental data for a linear C26000 system {T = 298K, v^ = 0.43) as reported by Michaels and Bixler [100] (experimental solubility data were calculated indirectly through S = P / D based on the D and P values of Ref 100). A close inspection of the data in Table III suggests that present TST simulations provide low-concentration solubility coefficients, S, which are in very good agreement with the experimental data for the whole range of gas penetrants, even for CO2. Past simulation works concluded that solubility coefficients, especially for large penetrant molecules (CH4 and CO2) are strongly and artificially affected by system size effects and, whenever small simulation boxes were used, predicted S values departed considerably from experimental values [68,101,102]. It appears that the present TST analysis on relatively large boxes (40-45A) rectifies related issues, although for a definite statement a more systematic approach is necessary using progressively smaller simulation cells. Diffusivity values are also in good agreement with experimental findings as long as the size of the penetrant molecule is small. For example, the best agreement is observed for oxygen (relative error of around 10%), while N2 diffusion rate is underestimated by the proposed TST approach within half order of magnitude. In contrast to the small gases, TST predicts CO2 diffusivity with a larger error: the simulated value is one order of magnitude smaller than the experimental one. It is clear that in its present formulation Gusev-Suter TST is not able to capture quantitatively the CO2 transport behavior due to its large and non-spherical size and to possible alternations of the surrounding polymer segments imposed by its presence that can not be correctly represented by the smearing factor. Finally, permeability coefficients of O2 and N2 (calculated indirectly as the product DS) were overand underestimated, respectively by similar factors in the range of half order of magnitude. Given the accurate TST-based predictions of sorption and diffusion of O2 and N2 through linear PE samples at standard conditions, in the next sections we will examine the effect of molecular weight (MW), temperature (7) and molecular architecture (short-chain branching, SCB) on the barrier properties of amorphous PE to these gases.

Atomistic simulation of the barrier properties

225

Table III. Simulation and experimental results for solubility, diffusion and permeability coefficients for O2, N2 and CO2 in purely amorphous linear PE matrices (T= 300K).

Penetrant O2 N2 CO2

S lO-^cm^(STP)/ cm^Pa) Simul. Exp. 1.14 ±0.08 0.418 ±0.08 3.20 ±0.3

0.828 0.400 4.47

D 10-'cm^ / s Exp. Simul. 13.2 ±1.2 4.84 ±0.8 0.41 ±0.1

12.1 8.42 9.79

P \0'^^ cm^STP) cm/ cm^ Pa s Exp. Simul. 15.0 ±0.6 2.02 ±0.09 1.31 ±0.1

10.0 3.37 43.7

IV,6.2. Effect of Molecular Weight Figures 7 and 8 present the solubility coefficient, S, of O2 and N2 as a function of molecular length for hnear PE systems as obtained from TST calculations at T = 450 and 300K. Regardless of the temperature, TST results suggest that 8(02) > S(N2) in the whole range of simulated chain lengths, in agreement with experimental trends [100]. The results also agree with semiempirical predictions based on the correlation function proposed by Teplyakov and Meares [103] according to which log(5) scales linearly with the effective well-depth parameter of the penetrant molecule (£(02) > e(^2) as seen in Table II). A qualitative difference is observed in the dependence of S on MW as a function of the applied temperature: at standard conditions {T = BOOK), S appears to depend weakly on molecular length and a maximum appears at A^av = 500. On the other hand, at higher temperature {T = 450K), solubility decreases hyperbolically with MW in a fashion reminiscent to the increment of density as a function of molecular length (see Fig. 4). As density increases, the magnitude and distribution of free volume decreases, consequently the number and size of the sorption sites that can accommodate the penetrant molecules are also reduced. Based on the results shown in Figs. 7 and 8 it can be concluded that MW has a profound effect on solubility at higher temperatures but at lower its effect diminishes.

P. Gestoso and N. Ch. Karayiannis

226

Molecular Length Figure 7. Low-concentration solubility coefficient, S, of O2 and N2 as a function of molecular length as obtainedfi-omTST calculations on polydisperse linear PE samples at 7 = 450K. 1.6

1

1

»

1

1—

•

1

'

1 —

»

1

Exp. O2

J O 1.4-1 1.2

N2

•

Exp. N2

\

it

1.0

\

« 0.8 J

U

\

ii> 0.6

\

0.4 0.2

1

'

1

200

'

400

1

1

600

800

•

1

'

1000

Molecular Length Figure 8. Same as in Fig. 7 but at 7 = 300K. Also shown are available experimental data after Ref 100.

Atomistic simulation of the barrier properties

227

1.5 1.4 H

m 0.

1.3

o #

N. Exp. N.

1.2 ^11-1 o Q

{

i

{

0.9-1 0.8

J

0.7 0.6 200

400

600

800

1000

Molecular Length

Figure 9. Low-concentration diffusion coefficient, Z), of O2 and N2 as a function of molecular length as obtained from TST calculations on polydisperse linear PE samples at 7 = 450K. Also shown is available experimental data for N2 diffusivity after Ref [104]. 4.5-1

1

1

11

1

1

— 1

'

4.03.5-

-r• 0

n •

\

3.0-

X

Q

1.5-

J J

J J

i

f

i • 1

5

1.0-

5 j

4

0.5f\0 ' 1 i)

Oa N, Exp. 0 , Exp. N,

i

M 2.5E u "0 2.0-

—r—

1

•

1 200

1

400

1

>

600

r

BOO

••••

1000

Molecular Length Figure 10. Same as in Fig. 9 but at 7 = 300K. Also shown are available experimental data from Ref 100.

P, Gestoso andN.Ch. Karayiannis

228

Figs. 9 and 10 depict the simulated diffusion coefficient, Z), for O2 and N2 as a function of molecular length at T = 450 and 300K, respectively. As already mentioned in the previous section, simulated difiusivities at room temperature are in very good agreement with experimental data [100]. Even at the higher temperature {T= 450K), the predicted diffusivity of N2 compares very well with the corresponding experimental measurements reported by Yoshiyuki [104] with the relative error being less than 10%. Based on this resuh it is interesting to notice that Gusev-Suter TST, even if it invokes simplifications and assumptions that have to be validated (in contrast to MD) and departs from the atomistic level by resorting to a coarser representation, is able to capture quantitatively the transport behavior of small gas penetrants in polymers, even in the melt state. Qualitatively, the transport behavior of both penetrants, O2 and N2, exhibits the same dependence on the average molecular weight of the polymer matrix at different temperatures: the kinetically-driven transport process is strongly correlated to the amount and distribution of free volume (and consequently density) and the reduction of the latter has a profound effect on the penetrant diffusion rate. TST-based results further suggest that between the two polymeric samples (C500 and Ciooo) gas diffusion is faster in the matrix with the higher segmental mobility (i.e. C500) revealing the effect of the dynamical flexibility of the polymer matrix on penetrant diffusivity [81]. Obviously, as shown in Figs. 9 and 10, gas diffusivity is also strongly affected by the size of the penetrant molecule: the larger the molecule, the slower its transport in the polymer, in agreement with the correlation function proposed by Teplyakov and Meares, valid over a very wide range of polymer/penetrant systems [103]. According to with this correlation, log(Z)) drops linearly with (f, where d is the effective diameter of the penetrant. Consequently, D{0^ > D(H^ for the whole range of simulated chain lengths and applied temperatures. Table IV. Low-concentration permeability coefficients of O2 and N2 in linear polydisperse PE systems as obtainedfiromTST simulations at r = 450 and 300K. P(10-^^cm^(STP) cm / cm^ Pa s) PE System

O2(r=300K)

N2(r=300K)

O2(r=450K)

N2(r=450K)

C78

3.85 ±0.12

0.657 ± 0.20

25.0 ±2.0

9.70 ±1.0

Cl42

3.05 ±0.10

0.457 ± 0.15

19.7 ±2.0

7.69 ± 0.8

C5OO

1.82 ±0.08

0.259 ±0.10

17.1 ±1.6

6.42 ±0.7

Ciooo

1.50 ±0.06

0.202 ± 0.09

15.8±1.8

5.83 ±0.7

Atomistic simulation of the barrier properties

229

The low-concentration permeability coefficient, P, for O2 and N2 in linear PE characterized by various molecular lengths as obtained from TST simulations at T = 450 and 300K is summarized in Table IV. The decrement in permeability with increasing chain length is mainly kinetically driven (decrease in difflisivity as shown in Figs. 9 and 10) attributed primarily to the reduction in the free volume and secondarily to the reduced segmental mobility of the polymer atoms (dynamic flexibility of the polymer host). Perm-selectivity of a polymeric membrane, a,/,, quantifying the separation efficiency of the material, of critical importance in gas separation processes, can be defined as the ratio of the low-concentration permeability coefficients of the two gases i andy P(i) ^ D{i) S{i)

Pij)

DU)S(j)

where dyj and Si/j are the diffusivity and solubility selectivities of the polymer host to the pair of gases, respectively. Based on the permeability coefficients exhibited in Table IV, at room temperature ao2/N2 = 7.4, which suggests that permeation of O2 in amorphous linear PE is preferred to N2 by a ratio of 7.4, which compares reasonably well with the available experimental value of ao2/N2 = 3.0 [100]. This preferential permeation of O2 against N2 in PE steams equivalently from kinetic (diffusion selectivity) and thermodynamic (sorption selectivity) factors. IV.6.3. Effect of Temperature Fig. 11 presents the dependence of the logarithm of solubility coefficient, log(5), for O2 and N2 on reciprocal temperature as obtained from TST calculations on a polydisperse linear C78 system. Penetrant solubihty decreases exponentially with increasing temperature as the gas molecule (O2 or N2) experiences more difficulties to dissolve in the polymer matrix as the temperature increases. On the other hand, as shown in Fig. 12, diffusivity decreases with decreasing temperature. TST-based diffusion coefficients of both O2 and N2 in Unear PE follow a distinctly non-Arrhenius temperature dependence, in agreement with similar observations for methane diffusion in PE from atomistic MD simulations of Pant and Boyd [24]. Table V summarizes selectivity ratios for solubility, soim, diffusivity J02/N2 and permeability aoimi in linear Cyg as obtained from TST simulations. Based on the selectivity data listed in Table V, it is evident that temperature decrease favors the preferential sorption and diffusion of O2 against N2 in amorphous PE.

230

P. Gestoso andN.Ch.

T—

• O . W U -1

-5.25 J

1•

1

— 1 — — r ^ 1—

1

1

0.

o

^.50 J

• J

^ ^.75 J

J

Q.

\

•^

1

J

o J

•

o ^.00 J o, -6.25-^

•

^ -6.50 J

]

o

J

•

^.75 J

• •

-7.00 J

•

o o o

-7.251.5

Karayiannis

1

r-

1

•

o o

o 1

—,— —1

2.0

2.5

r-

1

3.0

1

1

3.5

4.0

1000/T(K')

Figure 11. Logarithm of solubility coefficient, log(5), of O2 and N2 as a function of the inverse of temperature as obtained from TST calculations on linear C78 system. "0.\t

-

T

1

'

-3.5-]

1

5

»— 1 — 1 — •J- - 1 -

f

r

» -r

" 1

»

••

J

•

J

0

-4.0-

| -

•

*s

0

1 -^-S-

H

•

Q

0

0 .5.0-

-5.5-

•J

• 0. 0

N. -6.0- - T — T ""—, , 1.6 1.8 2.0 2.2

0 - I — ' — r — -1

2.4

2.6

1

2.8

r—1

3.0

1

1 1 1

3.2

1000/T(K')

Figure 12. Same as in Fig. 11 but for diffusion coefficient, D.

3.4

231

Atomistic simulation of the barrier properties

Table V. Solubility 5O2/N2> diffusivity doimi and permeabilityflfo2/N2selectivities of linear C78 PE to O2 and N2 as obtained from TST calculations at various temperatures. Temperature

(K)

^011^2

^02/N2

250

2.98

300

2.73

2.71

7.44

350

2.44

1.64

3.99

400

2.28

1.37

3.11

450

2.14

1.27

2.70

500

2.00

1.20

2.41

550

1.91

1.12

2.14

600

1.82

1.10

2.01

IV.6.4. Effect of Molecular Architecture (Short-Chain Branching) The effect of molecular architecture was studied by comparing the barrier properties of short-chain branched (SCB) PE systems against the linear analogs of the same total molecular weight (C142). The molecular characteristics for all simulated SCB samples are given in Table 11. NPT MD simulations of the equilibration phase revealed that short-chain branching has a rather minimal effect on density as depicted in Fig. 6, in agreement with reported experimental findings [93] in the purely amorphous melt phase. Figs. 13 and 14 present the dependence of solubility coefficient, S, on short-chain branching at r = 450 and 300K, respectively, hi general, short-chain branching produces an increment in gas solubility to a small extend with a maximum being observed for the nonlinear structure with branches bearing 5 carbon atoms (C5). The maximum difference in all cases between TST-based values of 5 does not exceed the range of 20-30%, consequently SCB has a rather small effect on solubility of gas molecules in PE.

P. Gestoso andN.Ch.

232

T C142

Karayiannis

-I 1 1 i r SBC_S)c14_tx2 SBC_7x1«_«x5 SBC_Sx22_4x«

SBC_11ic12_1dx1

Molecular Architecture

Figure 13. Effect of short-chain branching (SCB) on low-concentration solubility coefficient, S, as obtainedfi"omTST simulations on PE systems at r = 450K. I.O -1

•

1

1.6 J 1 • 0

J

1

'

\

'

1

«

N.

1.4 J

]

1 1.2 J

1

\ \

J

^ 1.0 J * *^ o

?

0.8

"l

J

0.

j

J

s ^ 0.6 J

J

J

i

\

0.40.2- ' — 1 — C142

1 < 1 SBC_11x12_1«x1 SBC_9x14_lx2

\

1-

SBC_7x1«_«x5

Molecular Architecture

Figure 14. Same as in Fig. 13 but at 7 = 300K.

1 ' SBC_5x22_4xi

233

Atomistic simulation of the barrier properties

T 1 1 1 r SBC_11x12_10x1 SBC_Sx14_»x2 SBC_7x16_«x5

SBC_53t22^4x«

Molecular Architecture Figure 15. Effect of short-chain branching (SCB) on low-concentration diffusion coefficient, D, as obtained from TST simulations on PE systems at r = 450K. ^.O '

1

2.4-

i

1

1

'1

1

1

11

r

-T

11

o

}

S* 1.6E o

A

2

IA

2.0-

N

N,

I

*

"o 1.2O

0.8-

i 5

0.40.0- 1 1 0142

,

•

§]

I

^

'•

1

SB C_11x12.10x1 SBC_Sx14_8x2

'

\

SBC_7x1i_6x5

Molecular Architecture Figure 16. Same as in Fig. 15 but at 7= 300K.

1

1

'

SBC_5x22_4x«

P. Gestoso andN.Ch.

234

Karayiannis

Figs. 15 and 16 present the dependence of diffusion coefficient on molecular architecture as obtained from TST simulations at T = 450 and 300K, respectively. In contrast to solubility, transport rates of penetrant molecules in short-chain branched systems appear to be somewhat smaller than those in linear counterparts characterized by the same molecular weight. Given that density is not significantly affected by the addition of branches (as shown in Fig. 6), an explanation can be given in terms of slower mobility of the polymer segments in SCB structures, especially in the vicinity of branch points as analyzed in chapter 11.6 of the book. The decrement in diffusion rates of the gas molecules, as a result of the short-chain branching, is more pronounced at lower temperatures. Low-concentration permeability coefficients, P, of O2 and N2 in PE matrices characterized by various molecular architectures are reported in Table VI. Based on the relative differences of the TST-predicted permeabilities for different structures it is evident that short-chain branching has a minimal effect on the barrier characteristics of purely amorphous PE systems. Accordingly, differences encountered between end-product LLDPE and HDPE materials could be attributed to the different degrees of crystallinity.

Table VL Low-concentration permeability coefficients of O2 and N2 in monodisperse short-chain branched (SCB) and linear PE systems as obtainedfromTST simulations at r = 450 and 300K. All systems are characterized by MW = 1990g/mol. P(10-^^cm^(STP) cm / cm^ Pa s) PE System

O2(r=300K)

N2(r=300K)

O2(r=450K)

N2(r=450K)

Cl42

2.59 ±0.10

0.382 ± 0.06

20.3 ±2.0

7.90 ± 0.8

SCB_llxl2_10xl

1.59±0.15

0.232 ± 0.08

19.4 ±2.5

7.64 ±0.9

SCB_9xl4_8x2

2.17 ±0.12

0.315 ±0.02

21.6±1.6

8.39 ±0.7

SCB_7xl6_6x5

2.80 ±0.15

0.359 ±0.06

20.9 ±1.8

8.13 ±0.7

SCB_5x22_4x8

2.39 ±0.10

0.381 ±0.05

21.7 ±3.0

8.66 ±0.6

Atomistic simulation of the barrier properties

235

V. Conclusions - Future Plans We have presented a study of the barrier properties of purely amorphous polyethylene samples to oxygen and nitrogen gases at infinite dilution as obtained from a hierarchical, multiscale modeling approach integrating state-ofthe-art simulation techniques. In the first stage, atomistic MC simulations built around chain-connectivity altering moves [69-74, an analysis of this technique is given in chapter 1.2 of the book] provide, within modest computational time, vast trajectories of representative and uncorrected confiigurations of the studied systems. Properly selected structures from the MC-generated trajectories with representative long-range characteristics, atomic packing and volumetric properties serve as excellent configurations for successive simulation studies based on Transition State Theory (TST) [39-42], as implemented in Insightll commercial software package [96]. Regarding computational time, the proposed hierarchical scheme, by integrating atomistic MC and Gusev-Suter TST (as implemented in Insightll), outperforms any conventional approaches by many orders of magnitude, especially for the modeling of the barrier properties of truly long, entangled systems (i.e. Ciooo PE) at low temperatures. Because of its simplifying assumptions, TST at its present formulation can be applied to small spherical gases at low concentration of the penetrant molecule. TST-based results revealed very good agreement between the calculated solubility, diffusivity and permeability coefficients of O2 and N2 in PE and available experimental data for a wide range of conditions. Sorption and diffusion are strongly correlated with the amount of free volume and the segmental mobility of the polymer atoms as identified by the effect of molecular weight on the barrier properties. Solubility decreases with increasing temperature as the penetrant molecule exhibits difficulties in condensing in the polymer matrix, while diffusivity shows a distinctly non-Arrhenius behavior as a function of temperature. Short-chain branching has a rather small effect on the barrier properties if the comparison is made between linear and branched structures bearing exactly the same molecular weight. Solubility in short-chain branched samples is to some extent favored when compared to linear analogs, whereas the opposite trend is observed regarding diffusivity. Finally, based on TST results, perm-selectivity of PE favors the permeation of O2 against N2 by a factor of around 7, in reasonable agreement with experimental trends which suggest a ratio of 3. As revealed from present findings, preferential permeation of O2 steams equivalently from kinetic (diffusivity selectivity) and thermodynamic (solubility selectivity) factors. Further work concerns the modeling of the barrier properties of glassy PE samples (at temperature in the vicinity and below Tg) and the expansion of the

236

P. Gestoso andN.Ch,

Karayiannis

proposed multiscale methodology to other chemically simple polymeric systems. Application of the present technique to chemically complex macromolecules requires the adaptation of a coarse-graining strategy where the atomistic structure of the polymer is mapped into a coarse-grained representation involving fewer degrees of freedom [105-109]. Li parallel, current efforts focus on the generalization of the TST approach to handle large, non-spherical gas molecules (for example CO2) by allowing the polymer matrix to undergo structural relaxation to accommodate the presence of penetrants. Acknowledgments P. Gestoso acknowledges Rhodia Recherches for permission to publish this work. N. Ch. Karayiannis expresses his gratitude to Prof. Doros Theodorou (University of Athens) for his guidance in the modeling of diffusion in polymers. Stimulating discussions with Prof. Manuel Laso (University of Madrid), Prof V. Mavrantzas (University of Patras), Prof. R. Gani and V. Soni (Technical University of Denmark), Dr. James Wescott (Accehys Ltd.), Dr. Nikolas Zacharopoulos and Niki Vergadou (National Research Center "Demokritos") and all partners of the PMILS project are deeply appreciated. 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 [1] R. Y. F. Liu, Y. S. Hu, M. R. Hibbs, D. M. Collard, D. A. Schiraldi, A. Hiltner and E. Baer, J. Appl. Polym. Sci. 98 (2005) 1615. [2] R. Y. F. Liu, Y. S. Hu, M. R. Hibbs, D. M. Collard, D. A. Schiraldi, A. Hiltner and E. Baer, J. Appl. Polym. Sci. 98 (2005) 1629. [3] J. H. Petropoulos, Polymeric Gas Separation Membranes, Eds. D. R. Paul and Y. P. Yampolski, CRC Press, Boca Ration, 1994. [4] J. S. and J. L. Duda, J. Polym. Sci. Polym. Phys. Ed. 15 (1977) 403. [5] J. S. Vrentas and C. M. Vrentas, Macromolecules 27 (1994) 5570. [6] J. S. Vrentas, C. M. Vrentas andN. Faridi, Macromolecules 29 (1996) 3272. [7] W. R. Vieth and K. J. Sladek, J. Coll. Sci. 20 (1965) 1014. [8] J. H. Petropoulos J. Polym. Sci. A-2 8 (1970) 1797. [9] D. R. Paul and W. J. Koros J. Polym. Sci. Polym. Phys. Ed. 14 (1976) 675. [10] J. H. Petropoulos, J. Membr Sci. 53 (1990) 229. [11] W. W. Brandt, J. Phys. Chem. 63 (1959) 1080. [12] R. J. Pace and A. Datyner, J. Polym. Sci. Polym. Phys. Ed. 17 (1979) 436. [13] R. J. Pace and A. Datyner, J. Polym. Sci. Polym. Phys. Ed. 17 (1979) 453. [14] R. J. Pace and A. Datyner, J. Polym. Sci. Polym. Phys. Ed. 17 (1979) 465. [15] P. Neogi, J. Polym. Sci. Polym. Phys. Ed. 31 (1993) 699.

Atomistic simulation of the barrier properties

237

[16] D. Homann, L. Fritz, J. Ulbrich, C. Schepers and M. Bohning, Macromol. Theory Simul. 9 (2000) 293. [17] F. Muller-Plathe, Acta Polymerica 45 (1994) 259. [18] E. L. V. Lewis, R. A. Duckett, I. M. Ward, J. P. A. Fairclough and A. J. Ryan, Polymer 44 (2003) 1631. [19] D. N. Theodorou, Diffusion in Polymers, Ed. P. Neogi, Marcel Dekker, New York, 1996. [20] H. Takeuchi, J. Chem. Phys. 93 (1990) 2062. [21] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [22] F. Muller-Plathe, S. C. Rogers and W. F. Gunsteren, Chem. Phys. Lett. 199 (1992) 237. [23] R. M. Sok, H. J. C. Berendsen and W. F. van Gunsteren, J. Chem. Phys. 96 (1992) 4699. [24] P. V. K. Pant and R. H. Boyd, Macromolecules 26 (1993) 679. [25] F. Muller-Plathe, S. C. Rogers and W. F. Gunsteren, J. Chem. Phys. 98 (1993) 9895. [26] Y. Tamai, H. Tanaka and K. Nakanishi, Macromolecules 27 (1994) 4498. [27] R. H. Gee and R. H. Boyd, Polymer 36 (1995) 1435. [28] C. S. Chassapis, J. K. Petrou, J. H. Petropoulos and D. N. Theodorou, Macromolecules 29 (1996)3615. [29] M. Fukuda, J. Chem. Phys. 109 (1998) 6476. [30] R. K. Bharadwaj and R. H. Boyd, Polymer 40 (1999) 4229. [31] J. H. D. Boshoff, R. F. Lobo and N. J. Wagner 34 (2001) 6107. [32] G. E. Karlsson, T. S. Johansson, U. W. Gedde and M. S. Hedenqvist, Macromolecules 35 (2002) 7453. [33] J. Budzien, J. D. McCoy and D. B. Adolf, J. Chem. Phys. 119 (2003) 9269. [34] A. Alentiev, L G. Economou, E. Finkelshtein, J. Petrou, V. E. Raptis, M. Sanopoulou, S. Soloviev, N. Ushakov and Y. Yampolskii, Polymer 45 (2004) 6933. [35] V. E. Raptis, I. E. Economou, D. N. Theodorou, J. Petrou and J. H. Petropoulos, Macromolecules 37 (2004) 1102. [36] L G. Economou, V. E. Raptis, V. S. Mehssas, D. N. Theodorou, J. Petrou and J. H. Petropoulos, Fluid Phase Equilibr. 228-229 (2005) 15. [37] N. Hu and J. R. Fried, Polymer 46 (2005) 4330. [38] D. Pavel and R. Shanks, Polymer 46 (2005) 6135. [39] S. Arizzi, Diffusion of small molecules in polymeric glasses: a modeling approach. Ph.D. Thesis, Massachusetts Institute of Technology, Boston, 1990. [40] A. A. Gusev, S. Arizzi, U. W. Suter and D. J. Moll, J. Chem. Phys. 99 (1993) 2221. [41] A. A. Gusev and U. W. Suter, J. Chem. Phys. 99 (1993) 2228. [42] A. A. Gusev, F. Muller-Plathe, W. F. van Gunsteren and U. W. Suter, Adv. Polym. Sci. 116(1994)207. [43] B. Widom, J. Chem. Phys. 39 (1963) 2808. [44] B. Widom, J. Phys. Chem. 86 (1982) 869. [45] N. C. Karayiannis, V. G. Mavrantzas and D. N. Theodorou, Chem. Eng. Sci. 56 (2001) 2789. [46] R. L. June, A. T. Bell and D. N. Theodorou, J. Chem. Phys. 95 (1991) 8866. [47] J. Karger, G. Fleischer and U. Roland, Diffusion in Condensed Matter, Eds. J. Karger, P. Heitjans and R. Haberlandt, Vieweg, Wiesbaden, 1998. [48] M. Ylihautala, J. Jokisaari, E. Fischer and R. Kimmich, Phys. Rev. E. 57 (1998) 6844. [49] N. Ch. Karayiannis, V. G. Mavrantzas and D. N. Theodorou, Macromolecules 37 (2004) 2978. [50] M. L. Greenfield and D. N. Theodorou, Macromolecules 26 (1993) 5461. [51] M. L. Greenfield and D. N. Theodorou, Mol. Simul. 19 (1997) 329. [52] A. A. Gray-Weale, R. H. Henchman, R. G. Gilbert, M. L. Greenfield and D. N. Theodorou, Macromolecules 30 (1997) 7296.

238 [53] [54] [55] [56] [57] [58] [59] [60]

P. Gestoso andN.Ch.

Karayiannis

M. L. Greenfield and D. N. Theodorou, Macromolecules 31 (1998) 7068. M. L. Greenfield and D. N. Theodorou, Macromolecules 34 (2001) 8541. D. Hofmann, J. Ulbrich, D. Fritsch and D. Paul, Polymer 37 (1996) 4773. D. Hofmann, L.Fritz, J. Ulbrich and D.Paul, Polymer 38 (1997)6145. L. Fritz and D. Hofinann, Polymer 38 (1997) 1035. L. Fritz and D. Hofinann, Polymer 39 (1998) 2531. D. Hofmann, L. Fritz and D. Paul, J. Membrane Sci. 144 (1998) 145. D. Hofmann, L. Fritz, J. Ulbrich, C. Schepers and M. Bohning, Macromol. Theor. Simul. 9 (2000) 293. [61] D. Hofmann, L. Fritz, J. Ulbrich and D. Paul, Comput. Theor. Polym. Sci. 10 (2000) 419. [62] M. Heuchel and D. Hofinann, Desalination 144 (2002) 67. [63] D. Hofmann, M. Heuchel, Y. Yampolskii, V. Khotimskii and V. Shantarovich, Macromolecules 35 (2002) 2129. [64] D. Hofinann, M. Entrialgo-Castano, A. Lerbret, M. Heuchel and Y. Yampolskii, Macromolecules 36 (2003) 8528. [65] M. Heuchel, D. Hofinann and P. Pullumbi, Macromolecules 37 (2004) 201. [66] J. R. Fried and P. Ren, Comput. Theor. Polym. Sci. 10 (2000) 447. [67] M. Lopez-Gonzalez, E. Saiz, J. Guzman and E. Riande, J. Chem. Phys. 115 (2001) 6728. [68] E. Kucukpinar, P. Doruker, Polymer 44 (2003) 3607. [69] P. V. K. Pant and D. N. Theodorou, Macromolecules 28 (1995) 7224. [70] V. G. Mavrantzas, T. D. Boone, E. Zervopoulou and D. N. Theodorou, Macromolecules 32 (1999)5072. [71] N. Ch. Karayiannis, V. G. Mavrantzas and D. N. Theodorou, Phys. Rev. Lett. 88 (2002) 105503. [72] N.. Ch. Karayiannis, A. E. Giannousaki, V. G. Mavrantzas and D. N. Theodorou, J. Chem. Phys. 117(2002)5465. [73] N. Ch. Karayiannis, A. E. Giannousaki and V. G. Mavrantzas, J. Chem. Phys. 118 (2003) 2451. [74] D. N. Theodorou, Bridging Time Scales: Molecular Simulations for the Next Decade, Eds. P. Nielaba, M. Mareschal and G. Ciccotti, Springer-Verlag, Berlin, 2002. [75] D. N. Theodorou and U. W. Suter, Macromolecules 18 (1985) 1467. [76] D. N. Theodorou and U. W. Suter, Macromolecules 19 (1986) 139. [77] Commercial simulation software Materials Studio (version 3.0) by Accelrys Inc. http://www.accelrys.com/products/mstudio/. [78] M. G. Martin and J. L Siepmann, J. Chem. Phys. B 102 (1998) 2569. [79] S. K. Nath and R. J. Khare, J. Chem. Phys. 115 (2001) 10837. [80] S. Toxvaerd, J. Chem. Phys. 107 (1997) 5197. [81] N. C. Karayiannis and V. G. Mavrantzas, Macromolecules 38 (2005) 8583. [82] V. G. Mavrantzas and D. N. Theodorou, J. Chem. Phys. 31 (1998) 6310. [83] L. R. Dodd, T. D. Boone and D. N. Theodorou, Mol. Phys. 78 (1993) 961. [84] S. J. Plimpton, J. Comput. Phys. 117 (1995) 1. [85] LAMMPS software distributed by Dr. S. Plimpton at Sandia National Laboratories, US. All equilibration NPT MD simulations were carried out using version LAMMPS 2001 (Fortran 90). [86] S Nose, Mol. Phys. 52 (1984) 255. [87] W. G. Hoover, Phys. Rev. A 31 (1985) 1695. [88] H. C. Andersen, J. Comput. Phys. 52 (1983) 24. [89] M. E. Tuckerman, B. J. Berne and G. J. Martyna, J. Chem. Phys. 97 (1992) 1990. [90] G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Kein, Mol. Phys. 87 (1996) 1117. [91] G. T. Dee, T. Ougizawa and D. J. Walsh, Polymer 33 (1992) 3462. [92] D. Cubero, N. Quirke and D. F. Coker, J. Chem. Phys. 119 (2003) 2669.

Atomistic simulation of the barrier properties

239

[93] J. L. Lundberg, J. Polym. Sci. Part A 2 (1964) 3925. [94] F. Muller-Plathe, S. C. Rogers and W. F. van Gunsteren, Macromolecules 25 (1992) 6722. [95] O. Ahumada, D. N. Theodorou, A. Triolo, V. Arrighi, C. Karatasos and J. P. Ryckaert, Macromolecules 35 (2002) 7110. [96] Commercial simulation software Insightll (version 400P+) by Accelrys Inc. http://www.accelrys.com/products/insight/. [97] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Di Nola and J. R. Haak, J. Chem.Phys. 81(1984)3684. [98] H. Sun, P. Ren and J. R. Fried, Comput. Theor. Polym. Sci. 8 (1998) 229. [99] H. Sun, J. Phys. Chem. B 102 (1998) 7338. [100] A. S. Michaels and H. J. Bixler, J. Polym. Sci. 50 (1961) 413. [101] T. R. Cuthbert, J. J. Wagner, M. E. Paulitis, G. Murgia and B. D. Aguanno, Macromolecules 32 (1999) 5017. [102] T. R. Cuthbert, N. J. Wagner and M. E. Paulaitis, Macromolecules 30 (1997) 3058. [103] V. Teplyakov and M. Meares, Gas Sep. Purif. 4 (1990) 66. [104] Y. Sato, K. Fujiwara, T. Takikawa, Sumamo, S. Takishima and H. Masuoka, Fluid Phase Equilb. 162(1999)261. [105] W. Tschop, K. Kremer, J. Batoulis, T. Burger and O. Hahn, Acta Polym. 49 (1998) 61. [106] W. Tschop, K. Kremer, O. Hahn, J. Batoulis and T. Burger 49 (1998) 75. [107] T. Aoyagi, F. Sawa, T. Shoji, H. Fukunaga, J. Takimoto and M. Doi, Comput. Phys. Commun. 145 (2002) 267. [108] F. Muller-Plathe, Chem. Phys. Chem. 3 (2002) 754. [109] N. Zacharopoulos, N. Vergadou and D. N. Theodorou, J. Chem. Phys. 122 (2005) 244111.