Macromolecules 1993,26, 619-686
679
Molecular Dynamics Simulation of Diffusion of Small Penetrants in Polymers P. V. Krishna Pant and Richard H. Boyd' Department of Materials Science and Engineering and Department of Chemical Engineering, University of Utah, Salt Lake City, Utah 84112 Received September 9, 1992 ABSTRACT Molecular dynamics simulations of diffusion of a small-moleculepenetrant in an amorphous polymer matrix have been carried out for the examples of methane in polyethylene (PE) and methane in polyisobutylene (PIB). Particular attentionwaa paid to ensuringthat the nonbonded potentiale representing the polymerpolymer bead interactions in the matrix were adequately calibrated. Accurate representation of equation of state (P-V-2') behavior for the polymers was used aa the criterion for this. The simulatione were carried out over a wide range of temperatures. The results are in accord with available experimental data. The temperature variation of the diffusion coefficients for PE was non-Arrhenius in character and of the WLF type. Analysis of the diffusion trajectoriesin PE indicates that the diffusion mechanism changes significantlyover the temperaturerange studied. A t lower temperatures,the diffusantmotion is characterized by relatively long periods of quiescenceinterspersed with fairly largejumps. However, at higher temperatures, a broad spectrum of frequent jumps obtains. The diffusion rates are obviously very sensitive to the free volume and hence to the molecular packing. However, the polymer chain mobility is found to play a role aa well. The diffusion activationenergy at higher temperatureia similar to that for the decay of the torsional angle autocorrelation function for the polymer. It ie also found that increaaing the torsional barrier reduces the diffusion rate. In PIB where the free volume ie significantly lees than that in PE, the diffusion rates are much lower than those in the latter. The mechanism over the range studied resembles that of the lower temperature regime in PE.
Introduction Molecular dynamics (MD) is a useful tool for exploring, via simulation, the structure and properties of bulk amorphous po1ymers.l The length of the trajectories that can be generated in practice presently is on the order of a nanosecond. Thus the range of properties that can be studied directly is limited to those that evolve over this time scale. One of the phenomena that appears to be suitable for investigation is the diffusionof smallpenetrant molecules in an amorphous polymer matrix. That is, the diffusion coefficientsof small penetrants in many rubbery or liquid polymers are such that, at temperatures in the vicinity of room temperature and above, the average displacement of the diffusant is large enough in a nanosecond interval to be determined via MD Simulation. Obviously the ability to effect such simulations would be of practical importance in predicting diffusion coefficients and also in understanding the mechanism of diffusion. However, there have been problems in actually carrying out these simulations reliably. Initial attempts to obtain diffusion coefficients of small penetrants through MD simulation generally2+ (but not alwayslO) led to results in poor agreement with experiment. The calculated values, especially in the example of polyethylene (PE), were very high in comparison with experiment and the calculated activation energies too low. The source of the discrepancy in PE has been traced not to the diffusion simulation per se but to simulation of the polymer matrix itself." The united-atom nonbonded potentials describing the interactions of the polymer CH2 beads were insufficiently accurate for the purpose. This was particularly evident in the volume of the polymer corresponding to a given T and P being inaccurate in simulation, or, conversely if experimental volumes were invoked, the Corresponding prassureswere highly negative. Diffusion is very sensitive to the system volume. It waa found that if an improved potential was invoked and wed in diffusion simulation at the volume corresponding to P = 1 atm, rather than the experimental volume, good agreement of calculated and
experimental diffusion Coefficients could be obtained. These results have been reported in a preliminary manner.11J2 Since then, the nonbonded potential for PE has been further improved. It can be calibrated so that it leads to a very good equation of state (P-V-T) representation, Le., an accurate volume versus temperature at constant pre~sure.'~ The purpose of the present work is to present in more detail diffusion simulation results that use this potential, to discuss the diffusion mechanism, and to extend the results to polyisobutylene (PIB). The latter polymer is of interest because, experimentally,the diffusion coefficients of small penetrants in it are known to be quite low in comparison with other polymers. Thus it is a good test case for diffusion simulations.
Simulation Details Systems Simulated. The simulationswere carried out on a single polymer chain in a periodic box. In both the cases of PE and PIB the chains contained 768 unitedatom centers, Le., [-CHz-l-l7~or [-CH~C(CH&-IIQ~.The systems were initially generated by MD runs on isolated chains and then imposition of the periodic boundary conditions. For example,in the case of PIB, the molecule was first simulated as an isolated chain at 600 K for 500 pa, starting from the all-trans configuration. The chain was then enclosed in a periodic box. The box was linearly shrunk to a length of 30.A over 200 pa. Following this step, the sysbm was cooled in NPT runs to 300 K in steps of 50 K for 100 ps each. Methane was chosen as the diffusant. This was done largely on the basis of it being the closest to spherical symmetryof the commonsmall-moleculepenetrants (other than the rare gases) and thus should be represented well in the united-atom approximation. Potential Functions. The development of the nonbonded potential functions for PE has been reported on in detail elsewhere.13 The united-atom (UA)LennardJones 6-12 nonbonded potential for CH2 groups includes an offset or displacement of the force center from the
0024-9291/93/2226-oS19~~0 ~ . 1993 ~ ~ ~ American Chemical Society
Macromolecules, Vol. 26, No. 4, 1993
680 Pant and Boyd
I .25
1
1.05
Table I Potential Funotionr.
function C-C bond stretch energy = 1/2k~(R- Ro12 bond bending energy = 1/2k,(B - 00)2 - C H r (PIB) - C H r (PE)
r
1.00
)C(
torsional potential in PIB = 1/2v3(1+ COS 34) torsional potential in PE = 1/2VS(l+ COB 34) + 1/2V1(1 + COB 4) UA nonbonded potentials, Lennard-Jones 6-12c CHz4
4
" " ' " " " ' ~ " " ' ~ ~ " ~ " ' '
Figure 1. Specificvolume of PIB vs temperature at zero pressure (points) as calculated from MD simulation. The curve is the smoothed experimental data of Eichinger and F10ry.l~
constant9 kR = 663) R o = 1.54
ke = 482, eo = 122.00 ko = 482,eO = 111.6O ke = 482, eo = 109.410 V3 13.4 V3 = 13.4, Vi
3.35
0.686.R-
p 3.940, 3.510, d = 0.42 t 0.391, R- = 3.812, a = 3.460, d = 0.0 c = 0.837,It4.224, a = 3.163, d = 0.42 t = 1-18,R= 4.21, u = 3.80, d = 0.0
t
u5
1.13 h
e
1.12
3
1.11
>
1.10
)C(
c\,
-CH3'
1
CH4
a Energies are in kJ/mol, distances in A, and anglee in rad (shown above in deg). Unless otherwise notad the potentials are from ref 9. Lowered by a factor of 4 from ref 9 in order to increasethe MD time step. For the - C H r group the potential in of the 'AUA" type; the interaction center is offset form the C atom by the distance = d along the bisector of the C-C-C angle in the direction of the hydrogens. The -CH3 potential is also AUA;the offset = d is along the C C bond. From ref 13. e From the present work.
*
1.08 107
0
200
400
600
800
loo0
1200
P (am) Figure 2. Specific volume of PIB vs pressure at 383 K (points) as calculated from MD simulation. The curve is the smoothed data of Beret and Prausnitz.l6
MD Methods Employed. MD techniques have been employedin thiswork for two distinct applications: studies of the equation of state of polymer melts and the analysis of small-molecule diffusion through the melh. For the carbon atom toward the substituent hydrogen atoms. Thus, first application, it was desired that the volume be it is of the "anisotropic" united-atom type introduced by determined from the conditions of impoeed T and P.The Toxvaerd14and is referred to here as an "AUA"potential. NPT method of N o s P was utilized. In theee -P-V-T" It is to be emphasized that the potential itself is isotropic equation of state simulations, all intramolecular modes with respect to its displaced center but appears to be were considered, including the bond stretches. The bond anisotropic with respect to the carbon center. stretch vibrations were weakened by a factor of 4 from In order to carry out computations for PIB, it was realistic strengths to enable an increase in the time step necessary to ensure that the additional united-atom from 0.5 fs in our previous workg to 1 fs here. Stronger nonbonded potentialsfor this system,those involving-CH3 forms of this modification have been widely used in and >C< centers, were properly calibrated. The >C< previous simulations of molecular system.'9 The altercenter is a single atom that does not involve the unitedation is believed not to strongly perturb chain dynamics. atom approximation, and it was felt that previous values The NPT dynamic equations were solved using a fivefor its parameters would be adequate. For the -CH3group value Gear predictor-corrector, selected specifically bethis was accomplished, as in the case of PE,13by equation cause of its general applicability to a wide range of MD of state and cohesive energy calculations in which the techniques. parameters were adjusted to obtain agreement with In the second application, to simulations of penetrantexperiment. polymer diffusion, the calculationswere carried out under The volume v8 temperature at zero pressure calculated effectivelyconstant volume (i.e., NVT) conditione. These from the adjusted AUA -CH3 potential is shown in Figure computations fell into two categories with respect to 1 in comparison with the data of Eichinger and F10ry.l~ temperature. The higher temperature diffusion reeulta A calculated V vs P isotherm at 383 K comparison with reported were calculated using N d s NPT method but the data of Beret and Prausnitd6 is displayed in Figure with values of the inertial parameter W that were large 2. enough to ensure that the volume did not change at all. The heat of vaporization is a valuable check on the In this manner, NVT simulations were recovered from potential parameters. Unfortunately, for PIB this quanthe algorithm. A single penetrant molecule was used. tity is not available, as it is for linear alkanes, from At low temperatures, however, reliable predictions of extrapolation of direct measurements on oligomem. Howthe diffusivity required enhanced diffusion trajectories ever, the solubility parameter, 6, of PIB as determined with longer simulation times and a larger number of from swelling experiments is available at room temperpenetrants. At 280 and 300 K in PE, eight methane ature. T w o experimental values reported from swelling penetrants were inserted into the system. At 360and 400 experiments for 6 in PIB at 298 K are 16.0 and 16.5 K in PIB, four methane molecules were introduced. In (J/cm3)'i2 (tabulation of Grulke"). The simulation for the adopted parameters gave a value of 16.0 ( J / c ~ ~ ) ' / ~the . multiple-diffusantsystems,it was ensured that no two penetrant molecules were less than 10A apart at the time All of the potential functions used are summarized in of insertion. The slow diffusive processes under lowTable I.
Diffusion of Small Penetrants in Polymera 681
Macromolecules, Vol. 26, No. 4, 1993
50.0
40.0
9
30.0 20.0 10.0
0.0 0
200
400
600
800
1000
lux)
t @s)
Figure 3. Mean-squared displacementfor methane penetrants in PE. The data were obtained from the trajecbriw of eight methane moleculwat 280 K. Thelatter is the lowest temperature studied. temperature conditionsprecluded the coordinatedmotion of two or more penetrants on the time scale of the MD runs. Bond constraints of the "SHAKE" f o r " were incorporated into all the diffusion simulations of PIB. In PE, C-C bonds were constrained only in the low-temperature multiple-penetrant cases. This permitted an increase in the MD time step from 1to 3-4 fs. (At temperatures of 350 K or lower, a time step of 4 fs was utilized; at higher temperatures, the lower value of 3 fs was employed.) A comparison of predicted diffusivitieswas made for a onepenetrant PE system at 400 K between bond-constrained and harmonic-bond models. It was found that the diffusivity for the constrained model (2.3 X lo4 mz/s) differed only slightly from that predicted using harmonically oscillating bonds (2.4 X lo4 mz/s). The difference was well within the margin of error inherent in predictions of diffusion constants. The SHAKE procedure was implemented with a leapfrog Verlet integrator and the loose-coupling NVT method of Berendsen et aLZ1 In all, the methane-PE system was simulated at eight temperatures: 280,300,330,350,375,400,450,and 500 K. PIB was simulated at 350,400,450,500,and 600 K. Simulations were not attempted for PIB at room temperature or lower because the low diffusion r a t a at these temperatures rendered reliable predictions impossible for practical run lengths. At each condition, the penetrantpolymer systems were fmt allowed to reach equilibrium over MD runs of 100-300 pa, followed by the tracking of diffusant trajectories. The lengths of the diffusion runs from which slopes were extracted ranged from 400 to 500 ps in the high-temperature cases of PE and from 1.2 to 1.5 ns in the low-temperature simulations of PE (300and 280 K) and PIB (400and 350 K). The length of the run was determined by the time required for the clear emergence of a long-time linear region in the mean-squared displacement curve (see Figure 3). The limited times accessible in MD make it necessary to extract the most possible informationfromthe penetrant trajectories in computingthe mean-squareddisplacement ( R W ) . To this end, several points along the penetrant path are treated as origins for new subtrajectories. Optimally, the time interval between successive origins would be given by a characteristic decay time within which the penetrant motion loses correlation with itself. This time would depend on the temperature of the system and on the dynamics of the specific polymer considered. In practice, penetrant positions were collected at intervals of about 1 pe. This time is significantlyshorter than that of the positional autocorrelation decay, ensuring that no
"diffusive" information in the trajectory is lost. The time dependence of the squared displacement was computed for each penetrant and each origin in the overalltrajectory; theaverageofthesefunctionsyielded,(Rz(t)). Anexample of the resultant curve for methane in polyethylene at 280 K (the lowest temperature considered) with eight penetrants and a run length of 1.2ns is shown in Figure 3.The self-diffusioncoefficients,under the conditionsprevailing here of no correlation between the penetrant motions,were obtained as one-sixth the slope of the linear part of the (Rz(t)) curve constructed as above. The diffusion curve (Figure 3) displays several characteristic features. The offset at zero time is due to the fast ballistic motion of the penetrant between colliiions with polymer atoms. Ita value dependson the temperature (which determines the average velocity) and the average intercobion distance. At intermediate times, in the cageeffect region, the motion of the penetrant is stillcorrelated and is not yet diffusive; the slope of the curve is higher than ita long-time value. At long times, (Rz(t)) is linear in time, as expected in a diffusive process. At very long times on the order of the total runlength,the smallnumber of subtrajectoriesover which averagingis performed results in noisy deviations from linear behavior. Further details may be found in the Ph.D. Dissertation of P.V.K.P.2z
Results and Discussion Comparison with Experiment. There are experimental measurements available for the diffusion of methane in PE. Since however over part of the range of data PE is a semicrystallinepolymer, comparisonof simulation results requires further elaboration. The simulationsrefer to a totally amorphous homogeneous material. There is no doubt that the diffusion actually takes place in semicrystalline samples almost exclusively through the amorphous fraction. Assuming that the local environment for diffueion in the amorphous fraction resembles that of the homogeneous amorphous state, the diffusion coefficient of the penetrant can be estimated fromthe compoeite material measurements if a reliable mixture equation can be invoked. As an example,a simple law of mixtures for the diffusion coefficient of the composite would be D = ulD1 + u& where u1 and uz are the volume fractions of the two phases and D1and D2 are the diffusion coefficienta of the penetrant in them. If 1 is taken to be the amorphous phase and 2 the crystalline one, then DZwould be expected to be zero. Actually however this equation would be only appropriate for diffusion through a medium that consista of parallel rectilinear channels of the two types of diffusion matrix. The lamellar structure found in semicrystallinepolymers would satisfythis relationonly for diffusion locally through regions where the lamellar surfaces are parallel to the diffusion direction. The diffusion tensor is a second rank one that connects the concentration gradient with the flux. Thus it has an analogy with the dielectric tensor connecting the electric field and displacement. In that case it has been shown that simple upper and lower bounds for the dielectric constant of a macroscopically homogeneous and isotropic material that is locally of lamellar structure can be derived.z3 Transcribing these results to the diffusion case gives the upper bound estimate of the effective composite diffusion coefficient as D = (20, where
+ D,)/3
(1)
Macromolecules, Vol. 26, No. 4, 1993
682 Pant and Boyd
For DZ= 0 and solving for D1,this leads to the following lower bound estimate for the diffusion coefficient in the amorphous phase, D1, in terms of the experimental one Dexp= D and the volume fraction of amorphous material u1
D1 = 3Dexd(2v1) (2) Because of the DZ= 0 condition invoked, as opposed to a finite nonzero value, a corresponding upper bound to D1 cannot be obtained. The experimental data of Michaels and for LPE and BPE and of Kanitz and for BPE were used to compute amorphousphase DIvalues from eq 2. Data taken for the melt obviously need no ~ o r r e c t i o n . ~The * ~simulation results for PE are compared to these experiments in Figure 4. Only two points are shown for each data source; these represent the temperature extremes of the data as smoothed and tabulated by Flynnm In passing it may be noted that the diffusion coefficient for the 78% (01 = 0.22) crystalline LPE, as corrected for crystallinity, falls between the two independent studies of BPE, the latter being 55 % crystalline (01 = -0.45). This would seem to indicate some validity for the correction method used. It is also to be noted that the dependence of D on temperature, i.e., activationenergy should be independent of the method or log D( Tz)/D(TI), of crystallinity correction since the latter is strictly geometrical in nature. It may be seen in Figure 4 that the agreement between the calculated and experimentaldiffusioncoefficientsover the complete temperature range is rather good. This is especially so since the derived experimental values represent lower bounds to D1. It is also to be noted that the temperature dependenceappears to be non-Arrheniusand of Vogel-Fulcherm or WLFm form. The non-Arrhenius behavior of methane-PE diffusionis not in iteelfsurprising. The diffusion curves for nitrogen, helium, and hydrogen in natural rubber and other elastomers display a similar convex curvature on an Arrhenius plot.31 Benzene in uncross-linked natural rubber dso displays very similar behavior.32 The curvature is consistent with free-volume models of diffusion such as that of FujitaF2 The VogelFulcher parameters in
1
-9.0
5
-9.5
n
CI
(4)
The parameters become TO= 248.5 K, c1 = 5.197, and c2 = 76.4 K. While the WLF nature does not directly imply that the diffusional process follows a free-volume dependence, it is usually considered to be characteristic of dynamic phenomena in vitrifying polymers. The mechanism of diffusion is discussed in greater detail in later sections of this chapter. The simulation results for PIB are shown in Figure 5. For comparison, the simulation results for PE are also shown. There do not appear to be data for methane diffusion in PIB for direct comparison with experiment. However, data are available%for the similar penetrants 0 2 , N2, and C02. This is also shown, at room temperature,
0 0 X
m
s" -10.5
+
0 simulation WLF fit
-11.0
-11.5
2.5
2
1.5
3
3.5
4
1W/T
Figure 4. Diffusion constante for methane in PE from both experiment and simulation. BPE = branched PE and LPE = linear PE. The data sourcesare 88 follow: h4B Michaels and LWH = Lundberg et al.;= Bider;%KH = Kanitz and L = Lundberg.n Data from MB and KH were taken in eemicryshlline samples, but the values plotted here have been corrected for crystallinity as discussed in the text and refer to the amorphous fraction. Data from LWH and L were taken in the melt.
-
-log D = A / ( T - T,) + B (3) were found by fitting the simulation data to be A = 397.0, B = 6.80, and T, = 172.1 K. In WLF form, where log DO = -12.0
t
.
-7.5
*+
-9.0
P
-9.5 -8'5 -10.0
-s" 8
I
I
,
I
I
I
1
I
.
,. .
I
1
1
,
I
-10.5
1.5
2
3
2.5
3.5
4
1W/T
Figure 5. Diffusion coefficients for methane in PIB from simulation (fiied circles). Also shown for comparieon are the simulation results of F w e 4for PE (opencircleah Experimental points at room temperatureare shownfor the diffusion of 0 2 , N2, and C02 in PIB (from a tabulation by StanneP). T h e dashed line representsthe experimental temperaturedependence of all three (a common activation energy of 50 kJ/mol).
in Figure 5. The experimental activation energies for all three are found to be the eame, and the temperature dependence is indicated by a single dashed line. It may be seen that the methane-PIB simulationresults could be expected to extrapolate reasonably well to the vicinity of these experimental points. It was found above that the temperature dependence for diffusion of methane in PE is of the WLF type. The case for this is not clear for PIB. A least-aquareequadratic polynomial fit is shown as the curve in Figure 5. It does show convex curvature. However, the slower diffusion in PIB restricts the temperature range, and the scatter is sufficient that for a limited number of simulation pointa WLF behavior cannot be reliably established. Short-TimeMotion of Penetrants. Atomistic simulations permit direct access to the trajectories of penetrants through the polymer amorphous phase. These trajectoriea shouldbe able to provide a mechaniiticpicture of the diffusive process. At very short times, trajectories may be analyzed in terms of the directional correlation function &(t) = d(t)-d(t-.z),where d(t) is the unit direction vector of the penetrant motion at time t and d(t-7) is the aame direction
Diffusion of Small Penetrants in Polymers 683
Macromolecules, Vol. 26, No. 4, 1993
0.70
I , , j 6ooo
6400
7200
6800
7m
3.00
8Ooo
I(fS)
Figure 6. Collisions in the trajectory of a methane penetrant in PE at 400 K. The directional correlation function is the dot product of the direction vector at time t with the direction vector 10 fs later. Portions of the curve close to 1represent straightline pathe; spikes mark collisions with polymer beads.
vector at an earlier time. For the purposes of this analysis, a value of T = 10fs is employed. Regions of the trajectory for which &(t) remains constant at unity denote straightline paths for the penetrant; sharp downward spikes in &(t) indicate collisions with surrounding polymer atoms involving sudden changes in direction. The function &(t) is plotted over a 2-ps period in PE at 400 K in Figure 6. The penetrant motion is clearly resolved into straightline paths that last for up to several hundreds of femtoseconds,interrupted by collisionswith polymer atoms that take about 20-30 fs each. Although it is possibleto associatethisshort-timemotion with the "rattling" of penetrant molecules in voids or "cages" in the polymer structure, the view of penetrant motion as short ballistic paths interrupted by collisions reveals little about the actual nature of diffusion. The mean-squared displacement curves (Figure3) from which diffusivitiesare predicted reveal a significant intermediate region at low temperatures before diffusive behavior sets in. This indicates that the penetrant positional autocorrelation funtion decays quite slowly at low temperatures, supporting the notion that the penetrant spends fairly long times in localized voids in the polymer before making diffusivejumps to other positions. Intercollisionpaths as obtained from &(t) cannot be clearly resolved into vibrations and diffusive jumps. Consequently, other means are required for the analysis of penetrant motion. Diffusion in Terms of Random Walks. As seen, the short-time trajectories of penetrants at low temperatures are dominated by rapid vibrational motions in localized regions in the polymer. In order to observe the nature of the diffusive process, it is necessary to filter out the effect of these vibrations. This is accomplished by reducing the trajectory to that between time-averaged positions over periods of length T , where 7 is large enough that successive steps are essentially uncorrelatsd with one another. The penetrant trajectory is thus reduced to a simple random walk of time step 7. In this abstracted view of the trajectory, it should be possible to predict the total diffusivity directly from the distribution of jump sizes. It should further be possible to determine the relative contributions of large and small jumps to the diffusivity. In the random-walk process the long trajectory may be broken up into intervals. Each interval consists of a number of elementary time steps that encompass a total time = T. The p i t i o n vector of the diffusant, Ri, has some average value over the interval, (Ri)r. Let ri be
2.60 2.80
, . ,
l . l l l l l l l l l l l l l l ,
1
1
1
,
1
1
q
I
1 \
1 \ simuktion "
4
"
"
~
6
^I
1
~
"
a
7
1
"
'
1
10
"
1
12
14
0s)
Figure 7. Onset of diffusion aa seen in random-walk analyeis at two temperatures. D (in 1o-B m2/s)is the diffusivity obtained by assumingthat the trajectory is a random walk between timeaveraged positions for a given time step. A t long enough averaging times the motion is truly diffusive and the random-walk analysis gives a value for the diffusivity close to the simulated value.
defined as the progress between adjacent intervals or, ri = I(Ri), - (&-i),I. Providing that 7 is long enough to average out nondiffusive motion, the expected meansquared displacement (R2(t)) for a random walk involving n steps of duration T in three dimensions with step sizes, ri, is n
( R 2 ( t ) )= /=I
The self-diffusion coefficient is therefore given by
The time period 7 = t in this analysis is raised in small steps from the time step of integration until the "randomwalk" diffusivity D, attains an equilibrium value, which shouldcorrespondapproximatelyto the value derived from the slope of (R2(t)).The value of 7 for which equilibrium is reached is related to the decay of the velocity autocorrelation function of the diffusing molecule. D, is shown as a function of T for two different temperatures in PE (300 and 400 K)in Figure 7 to illustrate this procedure. The equilibriumvalues of diffusivity from the random-walk analysis are seen to lie quite close to the value predicted from the slope of the mean-squared displacement. A mean-squared displacement curve generated from reduced trajectories is shown in Figure 8. It amply demonstrates that no diffusional informationis lost in the averaging process; only the nondiffusive ballistic and cage part of the diffusion curve is filtered out. The distribution of jump sizes for the reduced trajectories are shown in Figure 9. It is clear that, at the lower two temperatures, the penetrant spends most of its time in localized regions of the polymer without making large "jumps". The motion in these yslow"parta of the trajectory could arise from the self-diffusion of the penetrant in a void in the polymer structure and from the self-diffusion of the void itself. The larger jumps reflect sudden
Macromolecules, Vol. 26, No. 4, 1993
684 Pant and Boyd
60.0
I-
20.0
1
20.000
10.0
1
10.000
0
0
1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 I (PS)
Figure 8. Effect of averaging the penetrant trajectory over 16ps time steps at 300 K in PE. The curve with filled circle markers denotes the original trajectory;the one with open circle markers representsthe trajectory between 16-ps averages. The markers are placed only to distinguish the curves. .
1
.
1
1
1
1
1
1
1
1
.
,
1
1
.
1
.
.
,
1
1
.
1
1R
0250
f 0.200 0.150 0.100
0.050 O.Oo0
Figure 9. Distribution of diffusive jump lengths for methane in PE at three temperatures as computed from the filtered random-walk analysis. 30.0 25.0
Q 15.0 10.0
1
1 1
5.0
0.0 0
II
II
n A II
0.0
0.m
om
L [
1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 t @SI
Figure 10. Jump map or temporal history for a single penetrant in PE at 300K. The squared displacement is computed between successive 16-ps positional averages of the penetrant.
movement to new positions in the polymer. The two kinds of motion are more clearly illustrated in Figure 10 for a penetrant in PE at 300 K. The motion of the penetrant is seen to consist of infrequent large jumps separated by long periods of quiescence. Some further illustrative calculations are instructive. The diffusivity from the motion between jumps has been computed by excisingfrom the reduced trajectory all steps of length greater than 4 2 , where Ud is the Lennard-Jones diameter of the penetrant molecule. ‘Long” jumps were thus arbitrarily defined as those that exceed Ud in length. These large jumps occurred in only 9.1% of the 16-ps steps into which the trajectory was divided. The total diffusivity is 1.09X 10-1°m2/s,and the diffusivity “between
1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 t @8)
Figure 11. Jump map or temporal history for a singlepenetrant in PE at 400K. The squared displacement is computed between successive 7.6-ps positional averages of the penetrant.
jumps” was found to be 5.4 X 10-l’ m2/s. Thus at 300 K in PE, largejumps accountfor 51% of the diffusiveprogress of the penetrant. Significantly,the self-diffusion coefficient of the polymer beads, as determined by averaging the squared displacement of all carbon atoms over a 300psperiod,was3 X 10-11m2/s.Thus,thediffusionoccurring between large jumps is of the same order as the selfdiffusion of a polymer bead. The same analysiswas performed for PE at 280 K. The step size employed was 60 ps; the overall diffusivity was 3.3 X lo-” m2/s,and the diffusivity between jumps was m2/s. Largejumps thus accounted found to be 0.81 X for 76% ’ of the totaldiffusiveprogress. The self-diffusion coefficient for polymeric carbon a t o m was found to be 1.3 X lo-” m2/s,which again was of the same order as the diffusion rate for the permeant between jumps. Thia polymer atom self-diffusioncoefficientcorresponds to the diffusion coefficient of the polymer atoms relative to the chain center of mass and is a measure of their mobility on the time scale of penetrant diffusion. It is not the longtime self-diffusioncoefficient of polymer beada in the bulk, which must ultimately approach that of the polymer chain itself. At the higher temperature of 400 K, the penetrant trajectory is considerably different. Figure 7 indicates that the diffusiveregime seta in at very short times. Figure 9 shows a broad distribution of jump sizes. The high frequency of diffusivejumps with a conspicuous lack of quiescent periods is clearly seen in Figure 11. Steps of larger than ad2 were found to occur in 60% of the 7.5-ps time steps. The net diffusivity from Figure 4 is 2.4 X 10-O m2/s. Interestingly, steps smaller than ud2 contributed a diffusivity of 1.0 X lo4 m2/s, which was much higher than the polymer bead self-diffusion coefficient of 1.9 x 10-lom2/s. The view of the penetrant motion as that of “hopping”between voids is hence not applicableto these conditions;the effectivediffusive progress of the penetrant is seen to be much faster than that of polymeric atoms even on length scalesshorter than the penetrant diameter. The analysisof penetrant trajectories thus demonstrates a transition in the diffusive process with temperature. At low temperaturea the penetrant undergoes large “hops” between long periods of localizationin voids in the polymer structure. At higher temperatures, the polymer is more mobile. The penetrant is no longer “trapped” in voids, and the diffusive regime seta in almost immediately (see also Figure 12). For comparison, a similar analysis was performed of penetrant trdectories in PIB at 400 K. It was found that the onset of diffusion took a much longer time (Figure 13) than in PE at the same temperature. Further, the
Diffusion of Small Penetrants in Polymers 685
Macromolecules, Vol. 26, No. 4, 1993 l
l
l
)
.
.
l
(
.
l
.
l
.
,
.
(
l
.
.
300K
1.6 . , , . 1.4 1.2
l
l
l
l
l
l
,
l
l
l
l
,
l
l
l
l
l
l
l
l
.
l
l
l
l
l
-
-
-
1 -
n
: 0.6
-
0.4
-
0.8
0.2
'
&
' ' '
'
' ' ' '
box. The distance from the center of each cell on the grid
to the nearest polymer atom is referred to as Rmb. The cella of the grid were deemed to be unoccupied if rmi, >
'
' ' ' '
'
' ' ' '
'
' ' ' '
'
' ' ' "
4 2 , where upis the Lennard-Jonesdiameter of the nearest polymer segment. The fraction of the total grid that was unoccupied determined the total free-volume fraction vf. Not all the free volume so defined is available to a penetrant of finite size. The volume available to a penetrant of diameter Ud was obtained from the fraction of cells on the grid that satisfied rmi, > (up + ud)/2. This smaller set of cells was then resolved into interconnected clusters of free volume, by searching the set of "empty" cells for connections. The number and size distribution of these free-volume clusters yielded information on the static polymer structure as seen by the gaseouspenetrant. This entire procedure was applied to three cases: PE at 300 K, PE at 400 K, and PIB at 400 K. At each state, five periodic box snapshots were selected for analysis. In each set of five, no two structures were less than 100 ps apart alongthe same MD trajectory. The totalfree-volume fractions for the three conditions were found to be 0.368 f 0.001 (PIB, 400 K), 0.426 f 0.001 (PE, 400 K), and 0.389 f O.OO0 (PE,300 K). These correspond to specific free volumes (in cm3/g)of 0.421,0.534, and 0.458, respectively. Both the fractional and total free volume in PIB at 400 K are even lower than those for PE at 300 K. There was an averageof 2.4 free-volumeclusters per 300 K PE periodic box, each of volume 0.23 R3. The 400 K PE systems contained 6.0 voids of volume 2.135 A3 per structure analyzed, and the PIB structures at the same temperature contained an average of 4.4 clusters of size 1.57 A3. These results confirm the strong temperature dependence of the free volume available to a penetrant of the size of methane. They also verify that the methane molecules see far less free volume in PIB than they do in PE at the same temperature. The molecular structural reasons for the more efficient packing in PIB were investigated and elaborated on ear1ier.O The conclusions of that analysis are unchanged by the present work. The computations of free volume are based entirely on a static treatment. Even in the highest free-volume case (PE, 400K), a very small fraction of the volume is available statically to the penetrant in voidsof ita size-about 0.06 9%. This does not explain the rapid diffusive motion observed at these conditions in MD simulations. Clearly, the mobility of the polymer chains must permit a rapid dynamic redistribution of the free volume in order to achieve the fast diffusion seen at this temperature. This influence of polymer mobility on diffusion is discussed in the next section. Influence of Polymer Mobility. An important consideration in mechanistic descriptions of gas-polymer diffusion is the role of chain mobility. The influence of the torsional barrier has been examined in an earlier study of gas-polymer diffusion by Takeuchi and 0kazaki.B Penetrant diffusionthrough a torsionally realistic polymer model was compared in that study with diffusion through freely rotating chains. It was found that, although diffusion did proceed faster in the absence of torsional barriers, the effect was not very large. However, the simulations on which these conclusions were based employed previous UA models for the polymethylene beads. Aa commented upon in the Introduction these result in unrealistically fast permeation of penetranta. In the work here, the effect of raking the torsionalbarrier has been investigated. The experiment was performed on PE at 400 K. The form of the torsional potential, consisting of 1-fold and 3-fold terms, was unchanged, as
Macromolecules, Vol. 26, No.4, 1993
686 Pant and Boyd 1
consistent with the lower static free volume and the lower chain mobility as evidenced by the lower frequency of diffusivejumps (Figure 14) in PIB and by the decay of the torsionalautocorrelationfunction (Figure15). The latter is defined by Takeuchi and Okazakis as
.g 0.6
1
((cos 4(0) cos 4W) - (cos 4(0))2) (7) ((cos2$40))- (cos 4(0N2) where 4 ( t )is the torsional angle at time t and the brackets indicate an ensemble average of the angles.
R4(t) =
0.4
1B
0.2
0
1
1
50
loo
1
1
I
I
150
200
250
300
t (PS)
Figure 15. Decay of the torsional angle autocorrelation function for PIB at three temperatures (circles, 400 K; diamonds,350 K; triangles, 300 K). The curve at 400 K is a KWW function, R&) = exp(-t/r)@fit, where T = 1.66 ps and B = 0.79. Similar plota
for PE may be found in ref 13.
was the energy difference between the gauche and trans minima. the gauchetrans and gauchegauche’ energy barriers were altered by raising the 3-fold potential constant V3 by 50%. All other features of the molecular model were maintained intact. The diffusivity dropped considerably from the earlier simulated value (from a realistic torsional potential) of 2.4 X m2/s to a level where it could not be reliably determined from the single penetrant molecule used. The static free volumes in the two cases were identical. The change in the diffusion coefficient was thus brought about by the restricted mobility of the polymer chains under a higher barrier. It is clear, therefore, that static analyses of structures of the form used above are insufficient for the modeling of penetrant permeation. Cooperative motions between polymer segments and penetrant molecules would thus seem to play a key role in the mechanism of diffusion. However, the connection between chain dynamics and penetrant diffusion is not asimple one. It was noted above that the temperature dependence of the diffusion coefficient in PE is non-Arrhenius (Figure 5), a characteristic usually associated with vitrifying materials. However, there is no evidence for WLF-like behavior in the chain dynamicsof PE over the temperature range studied.13The decay of the torsional angle autocorrelation function is strictly Arrhenius over the range 300-450 K (with an activation energyof 25 kJ/mol).I3 In the case of penetrant diffusion, the activation energy at the higher end of the temperature range (from Figure 4) is similar to this value, 20 kJ/mol at 450 K, but rises well above it at low temperature, 42 kJ/mol at 300 K. A slowing of diffusion from Arrheniusbehavior, derivingfrom factors other than chain mobility, is therefore indicated. The analyses of penetrant trajectories,static free volume, and polymer mobility, in conjunction with the curvature seen in Figure 4, thus suggest a change in the diffusion mechanism over the range studied. At lower temperatures, penetrant motion is characterized by periodsof localization interspersed with large diffusive jumps. Diffusion at the lower temperatures therefore seems to be limited by the number of sites available to the penetrant in the polymer matrix and the distribution of intersite energy barriers. At higher temperatures, however, the picture of penetrant diffusion as a sequence of infrequent jumps breaks down; diffusion is limited mainly by the mobility of the polymer chains. In the case of PIB, diffusion throughout the measured temperature range (Figure 5) seems to proceed along the lines of the low-temperature regime in PE. This is
Acknowledgment. The authors are grateful to the Exxon Chemical Co. for financial support of this work and especially to Mr. J. V. Fusco and Dr. G,Ver Strate for their interest and helpful comments. We are also greatly indebted to the Utah Supercomputing Institute where the computations were performed. We appreciate the assistance of Mr. Jie Han in the calculations. References and Notes (1) Computer Simulation of Polymers; Roe,R. J., Ed.; PrenticeHa& Englewood Cliffs, NJ, 1991. (2) Trohalaki, S.; Kloczkowski,A.; Mark, J. E.; Rigby, D.; Roe,R. J. Reference 1, p 220. (3) Trohalaki, S.; Rigby, D.; Kloczkowski,A.; Mark, J. E.; Roe,R. J. Polym. R e p r . (Am. Chem. SOC.,Div. Polym. Chem.) 1989, 30, (2), 23. (4) Takeuchi, H. J. Chem. Phys. 1990,93,2062. (5) Takeuchi, H. J. Chem. Phys. 1990,93,4490. (6) Takeuchi, H.; Okazaki, K. J. Chem. Phys. 1990,92,5643. (7) Takeuchi, H.; Roe, R.J.; Mark, J. E. J. Chem. Phys. 1990,93, 9042. (8) Mueller-Plathe,F. J . Chem. Phys. 1991,94,3192. (9) Boyd, R. H.; Pant, P. V. K. Macromolecules 1991,24,6325. (10)An exception,where agreementwith experimentis obtained,is for the diffusion coefficient at 300 K of He and CHI in poly(dimethylsiloxane): Sok, R. M.; Berendsen, J. C.; van Gunsteren, W. F. J. Chem. Phys. 1992,96,4699. (11) Pant, P. V. K.; Boyd, R. H. Macromolecules 1992,25, 494. (12) Boyd, R. H.; Pant, P. V. K. Polym. R e p r . (Am. Chem. SOC., Diu. Polym. Chem.) 1992, 33 (l), 635. (13) Pant, P. V. K.; Han, J.; Smith, G. D.;Boyd, R. H., submitted to J. Chem. Phys. (14) Toxvaerd, S. J. Chem. Phys. 1990,93,4290. (15) Eichmger, B. E.; Flory, P. J. Macromolecules 1968, I, 285. (16) Beret, S.; Prausnitz, J. M. Macromolecules 1975,8, 536. (17) Grulke, E. A. Solubility Parameter Values. In Polymer Handbook, 3rd ed.; Brandrup, J., Immergut, E. H., Eds.; Wiley-Interscience: New York, 1989. (18) No&, S. J . Chem. Phys. 1984,81, 511. (19) Rigby, D.; Roe,R. J. J. Chem. Phys. 1987,87,7285. (20) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J . Comput. Phys. 1977,23, 327. (21) Berendsen, H. J. C.; Postma, J. P. M.; van Guneteren, W. F.; Di Nola, A.; Haak, J. R.J . Chem. Phys. 1984,66,2821. (22) Pant, P. V. K. Ph.D. Dissertation, University of Utah,1992. (23) Boyd, R. H. J. Polym. Sci., Polym. Phys. Ed. 1988,21, 605. (24) Michaels, A. S.;Bider, H.J. J . Polym. Sci. 1961, 50, 413. (25) Kanitz, P. J. F.; Huang, R. Y. M. J . Appl. Polym. Sci. 1970,14, 2739. (26) Lundberg, J. L.; Wilk, M. B.; Huyett, M. J. J . Polym. Sci. 1962, 57, 275. (27) Lundberg, J. L. J . Polym. Sci. 1964, A2, 3925. (28) Flynn, J. H. Polymer 1982,23, 1326. (29) Vogel, H. Phys. 2.1921,22,645. Fulcher, G. S. J. Am. Chem. SOC.1925,8, 339,789. (30) Williams,M. L.; Landel, R.F.;Ferry, J. D. J. Am. Chem. SOC. 1955, 77,3701. (31) Frhh, H. L.; Stern, A. CRC Crit. Rev. Solid State Mater. Sci. 1982, 11, 123. (32) Fujita, H. In Diffusionin Polymers; Crank,J., Park, G. S.,Eds.; Academic Preee: New York, 1968, Chapter 3. (33) Stannett, V. In Diffusion in Polymers; Crank, J., Park, G. S., Eds.; Academic Press: New York, 1968; Chapter 2.
