Gas Permeation in Semicrystalline Polyethylene as Studied by Molecular Simulation and Elastic Model

Résumé — Perméation de gaz dans le polyéthylène semi-cristallin par simulation moléculaire et modèle élastique — Dans ce travail, nous utilisons la simulation moléculaire pour étudier la perméation de deux gaz (CH 4 et CO 2 ) dans le polyéthylène. Ces simulations sont conduites à des températures inférieures à la température de fusion du polymère. Bien que dans de telles conditions, le polyéthylène soit à l’état semi-cristallin, des boîtes de simulation contenant exclusivement du polymère amorphe sont utilisées. Dans de précédents travaux [Memari P., Rousseau Polymer 51 , 4978], nous montré que les e ﬀ ets de la morphologie complexe des matériaux semi-cristallins pouvaient être pris en compte de manière implicite par une contrainte ad-hoc exercée sur la phase amorphe. Dans le présent travail, nous montrons que cette approche peut être mise en œuvre non seulement Abstract — Gas Permeation in Semicrystalline Polyethylene as Studied by Molecular Simulation and Elastic Model — We have employed molecular simulation to study the permeation of two di ﬀ erent gases (CH 4 and CO 2 ) in polyethylene. The simulations have been performed at temperatures below the polymer melting point. Although under such conditions, polyethylene is in a semicrystalline state, we have used simulation boxes containing only a purely amorphous material. We showed in previous works [Memari P., Lachet V., Rousseau B. (2010) Polymer 51 , 4978] that the e ﬀ ects of the complex morphology of semicrystalline materials on solubility can be implicitly taken into account by an ad - hoc constraint exerted on the amorphous phase. Here, it has been shown that our method can be applied not only for the calculation of equilibrium properties but also for transport properties like di ﬀ usion coe ﬃ cients. In addition, the ad-hoc constraint has been theoretically related to the fraction of elastically e ﬀ ective chains in the material by making use of Michaels and Hausslein elastic model [Michaels A.S., Hausslein R.W. (1965) J. Polymer Sci.: Part C 10 , 61]. We observe that the transport properties in amorphous regions are strongly governed by this fraction of elastically e ﬀ ective chains.


INTRODUCTION
In many areas of engineering, polymers are used as barriers to protect materials from gas or liquid contamination. For instance, polymers are used in food packaging, membrane separation processes, biomedical devices, offshore oil and gas production, etc. The relevant property for this process is the penetrant permeability which quantifies the flux of matter passing through the polymeric membrane. Permeability, Pe, is the product of two terms: the solubility, S , describing the solute concentration into the polymer phase and the Fickian diffusion, D, describing the solute mass transport inside the polymer: The knowledge of both solubility and Fickian diffusion coefficient (thereafter, called diffusion coefficient) is required for a complete description and understanding of permeability.
Besides experiments and theory, molecular simulation is an attractive tool to calculate such properties. Molecular simulation, by describing the system at the molecular level, may reveal the microscopic mechanisms that play an important role in the barrier properties. Thanks to the increase of computing power and to new methodological developments, several authors have investigated the modeling of gas solubility in polymers at a molecular level. For example, Gusev and Suter [19] computed the infinite dilution solubility of methane in polycarbonate. Müller-Plathe [20] computed Henry's constant of different gases in amorphous atactic polypropylene by using a test particle insertion method similar to Widom's method [21]. Later, de Pablo et al. [22] employed Gibbs ensemble Monte Carlo simulations in order to compute the solubility of short chain alkanes in polyethylene melts beyond Henry's regime. Several studies have followed, considering many different penetrant gases and more complex polymers such as polystyrene [23], poly(styrene-alt-maleic anhydride) copolymer, pols-(styrene-stat-butadiene) rubber and atactic polystyrene [24], polyimide [25], polyamide [26], polyethylene terephthalate [27,28].
In the early 1990s, Gusev et al. [29], Müller-Plathe [30,31], Müller-Plathe et al. [32], Takeuchi [33], as well as Pant and Boyd [34], used molecular simulation to investigate the diffusion of small penetrants in polymers. Pant and Boyd [34] studied the transport properties of small penetrants in amorphous polyethylene and polyisobutylene. Gusev et al. [29] and Müller-Plathe et al. [30,32] used molecular dynamics method to show that gases diffuse through polymers in a sequence of activated jumps between neighboring locations (called hopping mechanism). By using the same method, van der Vegt [35] achieved to observe a transition temperature between hopping-like and liquid-like diffusion.
Despite considerable work on solubility and diffusion coefficients, computer simulation of both quantities beyond infinite dilution limit and using the same molecular model are rare. In this paper, we will compute simultaneously solubility and diffusion coefficients in order to evaluate the permeability. Furthermore, in most of the previous works mentionned above, transport properties have been obtained in melt or purely amorphous polymer or, in the case of semicrystalline materials, within the assumptions that: the amorphous phase is the only permeable phase; the amorphous phase characteristics are not affected by the presence of the crystalline regions. The first assumption is strongly supported by experimental evidences. Experiments by Michaels and Bixler [5,6] on different polyethylene grades with different degrees of crystallinity permitted to establish the following relationship between the solubility, S , and the diffusion coefficient, D, in the semicrystalline material versus the solubility and the diffusion coefficient in the amorphous phase (denoted by subscript a): where χ is the (volume or mass) fraction of crystalline regions. In Equation (3), τ is a tortuosity factor accounting for the increase of the diffusion path caused by the presence of crystallites and β is a chain immobilization factor which takes into account the reduction of the amorphous chain segment mobility due to the proximity of crystallites. The second assumption is not always valid. In semicrystalline materials, the amorphous phase can be perturbed by the presence of crystallites. Thereby, several authors have reported large deviations between experimental data and theoretical predictions when they assume that the amorphous phase properties are not affected by the crystalline regions. Van der Vegt et al. [36] showed that gas solubilities in amorphous polyethylene, computed by molecular simulation, are 1.4 to 3.5 times greater than experimental values (S a ) corrected from amorphous content (Eq. 2). Using another molecular model, Nath and de Pablo [37] found that the solubility of small molecules in amorphous polyethylene is overestimated, especially for highly soluble gases. Hu and Fried [38] also reported solubility coefficients of small gas molecules in poly(organophosphazenes) five times larger than what would be expected by extrapolating values reported for semicrystalline samples to a 100% amorphous material.
The crystalline region effect on the amorphous phase has been accounted for by theoretical approaches [39][40][41] based on Flory and Rehner work [42,43]. This effect comes from the fact that some polymer chains leaving a crystalline region may be trapped in another crystalline region. The network formed by these tie segments may constrain the amorphous regions of the polymer. This effect, called the elastic effect, has been invoked by several authors to correct the predicted solubilities by equations of state or group contribution methods [44][45][46][47] but was never accounted for in molecular simulations.
In a previous work [48], we were able to implicitly account for the effect of crystalline regions upon gas solubilities by using Monte Carlo simulations in the osmotic ensemble in an original way. We reproduced experimental solubilities by exerting an additional uniform constraint on our purely amorphous simulation box. A single constraint value emerged, independent of the gas nature, characteristic of the semicrystalline material. We concluded that the role of this ad-hoc constraint is to reproduce the effective density of the permeable phase of the material. This method has also shown its success toward the computation of gas mixture solubilities in semicrystalline polyethylene [49].
The methodology mentioned above has been applied to polyethylene (PE), because of its wide use in industrial applications and its relative simplicity as a molecular model. In the present work, we extend the proposed approach to permeability calculations in semicrystalline PE. Two different penetrants are investigated: carbon dioxide and methane, mostly chosen because of their relevance in gas production. The obtained results are then validated against experimental data.
This article is organized as follows. In the next section, we rapidly review the methodological aspects of this work. Section 2 discusses estimations of transport properties in the amorphous region of semicrystalline polymer from experimental data. In Section 3, we present our simulation results on gas solubility, gas diffusion coefficient and gas permeation. The increase of diffusion coefficient with gas concentration in polymer is captured. From the study of transport property temperature dependency, activation energies are calculated and compared with experiments. A link between the ad-hoc constraint used in our simulations and the elastic effect is then provided in Section 4. Finally, Section 5 summarizes our findings.

METHODOLOGY
We have employed Monte Carlo (MC) simulations in the osmotic ensemble to obtain the solubility of different gases in polyethylene. The details of these MC simulations are given elsewhere [48]. Molecular Dynamics (MD) simulations were performed in the canonical ensemble to obtain the diffusion coefficient of penetrants. Nosé-Hoover thermostat with an explicit time reversible integrator [50] has been used. A timestep of 2 fs was used in our simulations. The carbon-carbon bonds in polyethylene were constrained to their equilibrium value (1.535Å) using the Rattle algorithm [51].
Initial MD configurations are issued from the output of the osmotic ensemble MC simulations with number of penetrant molecules and volume equal to the average MC values. Initial particle velocities are taken from a Gaussian distribution. Self-diffusion coefficients, D * a , are computed from the mean square displacement msd(t) of the N penetrant molecule center of mass positions, r cm i (t): using the Einstein equation [52]: For each system, the production run lasts 80 ns. During this period, the mean displacement was larger than ten molecular diameters and we checked that penetrant molecules reached the diffusive mode (i.e. the linear regime of msd(t)). The polymer studied in this work is a model of linear polyethylene (PE) with 70 carbon atoms per chain. The force field and the mixing rule are identical to the ones used in our previous works [48,49].

ESTIMATION OF TRANSPORT PARAMETERS IN THE AMORPHOUS PHASE FROM EXPERIMENTS
There are experimental measurements available for solubility and diffusion coefficients of gas in semicrystalline PE. However the simulations refer only to the amorphous part of the material, thus comparison of simulation results with experiments requires further elaboration.
Thanks to Equation (2), the conversion of measured solubility coefficients S to solubility coefficients in amorphous phase S a is straightforward. The only required parameter is the crystallinity ratio χ which can be measured independently by several techniques (Differential Scanning Calorimetry (DSC), density measurement, Raman spectroscopy, Transmission Electron Microscopy (TEM) and Nuclear Magnetic Resonance (NMR) [53]).
The determination of the diffusion coefficient in the amorphous region, D a , from experimental value D, is not straightforward. Parameters τ and β involved in Equation (3) are not easily measured experimentally. Therefore, we could only obtain estimates for D a from three different methods.
In the first method, D a is extrapolated from the diffusion coefficient in polymer melts using the following equation: where T is a temperature above the polymer melting point, D is the diffusion coefficient in polymer melt at T , ΔE d is an activation energy for diffusion and T is an arbitrary temperature below the melting point. We assumed that ΔE d is constant over the whole temperature range. The second method uses diffusion coefficient values in a 100% amorphous polymer grade. For polyethylene, there isn't a 100% amorphous grade at low temperature. Other polymers, such as natural rubber, are usually considered as its amorphous analogue [6]: The third method, proposed by Pant and Boyd [34], gives a lower bound for the diffusion coefficient using the following equation: In the present work, we compare our simulation results with the values estimated from these three methods (Eq. 6-8).

SIMULATION RESULTS
In Figure 1, we present CO 2 and CH 4 concentrations in the amorphous phase of polyethylene at 293 K and 298 K respectively, obtained from Monte Carlo simulations in the osmotic ensemble. During these simulations, an additional isotropic constraint σ PE of 80 MPa is imposed on the simulation box to account for the effect of crystallites by reproducing the average density of the permeable phase [48]. As can be seen, a good agreement between simulation and experiments is observed. It is also shown that the gas concentration in polyethylene is linearly proportional to the pressure in this pressure range. In such cases, the Fickian diffusion coefficient (D a ) is identical to the Maxwell-Stefan diffusion coefficient [54]. In addition, if the polymer chain mobility is neglected comparing to solute velocities, the Maxwell-Stefan diffusion coefficient can itself be approximated by the penetrant self-diffusion coefficient (D * a ). (For Fig. 1, see also [5,[55][56][57][58]). Therefore, diffusion coefficient D a can be calculated from Equation (5). A msd(t) curve and its time derivative for CO 2 diffusion in PE is given in Figure 2. The average diffusion coefficient is computed in the region where the slope of msd(t) is constant.
In Figures 3 and 4, we present the simulation results for CO 2 diffusion at 293 K and CH 4 diffusion at 298 K. (For Fig. 3, see also [6,7,59,60] and for Fig. 4, see [6,61]). It is shown that the use of an additional constraint of 80 MPa on the amorphous box leads to a decrease of the diffusion coefficient. This observation is related to the increased density of the amorphous region under such an additional constraint, reducing gas mobility.
Experimental data are also reported in these graphs. Laguna et al. [7] reported two different experimental values  for CO 2 diffusion coefficient: a diffusion coefficient measured by permeation cell technique and a self-diffusion coefficient measured by NMR. The value obtained by NMR is an order of magnitude higher than the one obtained by permeation experiments. Except this NMR value and the estimate from the extrapolation method, all experimental values corrected either by Equation (7) or by Equation (8)  Diffusion coefficient of methane versus concentration in the amorphous phase of PE at 298 K. Two series of simulations are presented: with (black) and without (red) an additional constraint of 80 MPa. Dashed lines are exponential fits of the simulation results (Eq. 9). Experimental data are from Michaels and Bixler [6] and Lundberg [61]. a Estimated from Equation (6), b estimated from Equation (7), c estimated from Equation (8). need of an additional constraint to allow for a more realistic description of the amorphous region. It is also concluded that the extrapolation method is not reliable to estimate the diffusion coefficients at low temperatures. Main reasons are: the diffusion activation energy is not constant with temperature; the effect of crystalline regions on the gas diffusion in the amorphous phase is completely ignored by the extrapolation approach. We also observe that the diffusion coefficient increases with gas concentration in the polymer. This is related to the plasticization of polymer chains due to the presence of penetrant molecules. Vittoria [62] describes this behaviour using an exponential function: In this relationship, D ∞ a refers to the diffusion coefficient at infinite dilution which is related to the fraction of free volume and γ is the concentration coefficient which is a measure of the plasticization effect. D ∞ a and γ values obtained from our simulation results are given in Table 1.
A diffusion slightly faster is obtained for CO 2 compared to CH 4 . Such behaviour has been already reported in the literature with other polymers: polypropylene [63] and polystyrene [64], and has been interpreted on the basis of molecular structure considerations, attributing the faster diffusion to faster axial motions in case of linear molecules.
Finally, permeabilities have been computed from our molecular simulation diffusion and solubility coefficients. Calculated values are presented in Table 2 and shown in Figure 5 for the two studied gases. Data are in good agreement with experimental results of Michaels and Bixler at infinite dilution [6].
As the morphology of semicrystalline polyethylene evolves with temperature, we showed in our previous works [48,65] that the additional constraint must decrease with increasing temperature: from 80 ± 10 MPa at room temperature (293-298 K) to 60 ± 10 MPa at 333 K and 40 ± 10 MPa at 353 K. We have thus computed CO 2 and CH 4 solubility and diffusion coefficients at these three temperatures (see Tab. 2, 3) and we have then estimated activation energies for solubility, diffusion and permeability processes (see Tab. 4). Calculated activation energies are compared with experimental values of Flaconnèche et al. [66] and Michaels and Bixler [6]. The orders of magnitude are well reproduced. However, the activation energy for the diffusion process are slightly underestimated by our simulations. This difference can be related to the chain immobilization factor (β) whose contribution to the activation energy is not taken into account in our simulations.

INTERPRETATION OF THE ADDITIONAL CONSTRAINT BASED ON THE ELASTIC MODEL
The non crystalline regions of semicrystalline polymer are composed of elastically effective and elastically ineffective chain segments. The elastically effective chain segments consist exclusively of intercrystalline tie segments. The network formed by these tie segments has been considered as crosslinks constraining the amorphous phase of the polymer. This network results in an additional elastic constraint [39][40][41]. Several authors have invoked this elastic effect to correct their calculated gas solubility coefficients in the amorphous phase of semicrystalline polymers. The correction term is usually introduced as a contribution to the activity coefficients of the different species. Some modified equations of state are proposed that allow the calculation of thermodynamics properties including the elastic effect contribution [40,[44][45][46][47]67]. Different elastic models exist in the literature. The Flory and Rehner model, which is dedicated to the study of crosslinks rubbers [40,42,43], assumes a Gaussian distribution of the tie segments. A few years later, this model has been extended to semicrystalline polymers [39][40][41] considering a Hookean behaviour of the tie segments. Some authors have tried to compare these two elastic models [44,47]. Doong and Ho [44] showed that the model Permeability in PE amorphous phase for different pressures of carbon dioxide (top) and methane (bottom) at room temperature. Simulation results are obtained with an additional constraint of 80 MPa. Experimental data are from Michaels and Bixler [6]. b Estimated from Equation (7), c estimated from Equation (8).
proposed by Michaels and Hausslein [39], especially due to its temperature dependency, is more relevant in semicrystalline PE. In the elastic model with Hookean segments, the energy cost associated with the stretching of an elastically effective tie segment is described by [39,46]: where r is the average distance between adjoining crystallites, u the average number of methylene units in the tie segment and K a spring constant. The force due to an elastically effective tie segment, f elastic a , is given by −∇(ΔU elastic a ). By summing over all ν elastically effective segments, the total force on the polymer amorphous phase can be obtained: Kr ur (11) wherer is the unit vector in the r direction. The constraint on the amorphous phase due to the elastically effective tie segments is then: where V a is the total volume of the amorphous regions. Michaels and Hausslein [39] and more recently Banasazak et al. [46], based on the thermodynamic equilibrium between crystalline and amorphous phases, have shown the following relationship (see Eq. 26 in [39] where ΔH f is the polymer molar heat of fusion, v P a the partial molar volume of polymer in the amorphous polymer phase, φ the volume fraction of polymer in the amorphous polymer phase, T m the melting point of crystallites in the presence of penetrant molecules and f the fraction of elastically effective chains in the amorphous regions. Substituting Equation (13) into Equation (12), we obtain the following expression for the additional constraint: It is interesting to note that Equation (14) predicts a nonzero constraint value even at zero penetrant concentration, i.e. for pure polymer. In this case, φ = 1, v P a = 1/ρ a where ρ a is the number density of the pure amorphous region and T m = T m , the melting point of pure semicrystalline polymer. Equation (14) can thus be rewritten: In the following, we compare the values of σ PE computed using Equation (15) with the values used in our molecular simulations. We recall that our constraint values were calibrated from comparison between experimental CO 2 solubiliy coefficients on PE samples with crystallinity χ = 50% and simulated values computed in the osmotic ensemble [48]. All quantities required in Equation (15), except f , have been taken from experimental results of Flaconnèche et al. obtained on a medium density polyethylene sample [66]. Few authors have reported values for the empirical factor f in order to reproduce solubility data of penetrants in different PE samples [39,44,45,47]. f values are typically in the range 0.27-0.50. In Figure 6, we present the additional constraint value, σ PE , calculated from Equation (15) as a function of f at 298 K, 333 K and 353 K. At these temperatures, the additional constraint values used in our simulations are respectively 80 ± 10 MPa, 60 ± 10 MPa and 40 ± 10 MPa [48]. Figure 6 shows that a value of f , around 0.65-0.75, can reproduce the additional constraints at all temperatures. This fraction of elastically effective chains is slightly larger than the values reported in the litterature. Such a large value might be related to the use of rather small chain lengths (70 carbon atoms) to model polyethylene. The use of longer chains would indeed lead to lower additional constraint values, and thus to lower fractions of elastically effective chains.

CONCLUSIONS
We have shown in this work that our recently proposed technique based on osmotic ensemble Monte Carlo simulations [48] provides an effective method for evaluation of gas solubility in the amorphous phase of semicrystalline polymer. The real material density is reproduced when an ad-hoc constraint is applied to the polymer amorphous phase. This density value can then be imposed in molecular dynamics simulations without the use of a priori experimental results. The resulting diffusion coefficients are in good agreement with those obtained from experiments (Michaels and Bixler [6] among others). Our methodology permitted the computation of the transport parameters (diffusion coefficient and permeability) of polyethylene amorphous region toward carbon dioxide and methane. The study of diffusion coefficient at different gas concentrations allowed us to capture the plasticization effect in polyethylene. In addition, we reproduced reasonable activation energies.
Elastic model of Michaels and Hausslein [39] is seen to be able to justify the effective constraint values employed in the simulations. It is shown that intercrystalline tie segments exert an additional constraint on the polymer amorphous phase, what should be accounted for in any molecular simulation of the amorphous phase. Our findings suggest that the fraction of elastically effective chains is independent of the gas pressure (in the range of 0-8 MPa) or nature (for CO 2 and CH 4 ) or temperature (in the range of 293-353 K). Rather, it is a characteristic of the polymer with a given crystallinity and history. For this approach to be fully predictive, an a priori knowledge is required for this fraction of elastically effective chains.

ACKNOWLEDGMENTS
PM would like to thank IFP Energies nouvelles for financial support. Jean-Marie Teuler (Laboratoire de Chimie Physique) is gratefully acknowledged for support in code development.