Chapter 11. Diffusion of Gases through Polymers

11.1. Introduction

Controlling the diffusion of gas molecules through polymer membranes is a key problem in many industrial processes. It is desirable to predict the diffusion coefficient for a given gas through a polymer of given composition. If quantitative diffusion coefficients cannot be obtained, relative diffusion rates for differing gases are still very useful.

It is generally thought that crystalline regions of the polymer matrix are too dense and closely packed to allow penetration by the diffusant molecule. The amorphous or amorphous-crystalline interface regions apparently control the diffusion process. Accurately simulating diffusion at an atomistic level thus requires the simulation of amorphous polymers.

Since amorphous regions of a polymer are not well-ordered, use of a single, short chain or chain fragment and periodic boundary conditions, as is typically done in order to limit the size of the simulation, is likely to yield inaccurate results due to the imposed short-range order. Using multiple, longer chains in a much larger unit cell should give a better approximation of the true bulk amorphous behavior. This then requires large-scale molecular mechanics, in the regime from thousands to millions of atoms, the latter particularly in cases where properties are especially sensitive to molecular weight.

Key industrial gases for which experimental data are often available include CO2, O2, N2, and He. Important polymers include poly(ethylene), poly(vinyl chloride), and poly(vinylidene chloride), as well as copolymer mixtures of the latter two.

11.2. Procedure

An amorphous polymer unit cell was generated from a single chain of poly(ethylene) of 34.5 monomer units built using the BIOGRAF amorphous builder routine [1], which uses the rotational isomeric state (RIS) methodology. This single chain was equilibrated at 350 K using Nosé-Hoover constant-temperature, constant-pressure (TPN) dynamics. The final cell volume gave a density of 0.87 g/cc.

This unit cell was then replicated to build a 5x2x1 supercell consisting of 2090 atoms in 10 polymer chains. The dimensions of the supercell were approximately cubic: 27.429 x 27.046 x 27.255 Å, with angles of 90.0 deg. The supercell was relaxed using Nosé-Hoover TVN dynamics (fixed cell parameters) at 350 K for 50 ps. The resulting structure was highly disordered, with no apparent residual crystallinity. This structure is depicted in Figure 11-1.

The resulting unit cell was then scanned for voids using a program that attempts to fit a sphere of a given radius at each point on a grid within the unit cell. For CO2, the sphere size was chosen to be 3.10 Å; for He, 2.5 Å spheres were used. Random voids of adequate size were chosen as the locations of the diffusant gas molecules. Figure 11-2 shows the starting position of the CO2 molecule, along with the surrounding void spaces.

One CO2 molecule or three He atoms were placed in the unit cell. The CO2 diffusant molecule was treated as rigid using quaternions.

Dynamics was performed at 350 K for 300 ps using 1 fs timesteps; the location of the diffusant molecule was tracked at each 0.1 ps interval.

The forcefield used had been optimized for single-chain periodic poly(ethylene) [2]. The CO2 van der Waals parameters used were taken from the DREIDING II forcefield [3], with charges assigned so as to reproduce the experimental quadrupole moment. A CMM level of 3 was used; the farfield update frequency was every 10 steps. The Nosé-Hoover \tau_s constant was set to 0.01 ps.

11.3. Results

Figure 11-3 shows the path of the CO2 molecule through the unit cell during the diffusion process. The path appears to be composed of a few sections in which the molecule remains localized to an area, or trapped, interleaved with sections in which the molecule moves between areas. The lines crossing the path indicate the orientation of the CO2 molecule at 10 ps intervals through the trajectory. Unexpectedly, it appears that the axis of the molecule is generally perpendicular to the direction of motion.

Figure 11-4 shows the path of the He atoms. In this case, the more mobile He atoms appear to avoid being trapped by the polymer. Note that the three He atoms went off in different directions, and that over the course of the simulation, they managed to drift through multiple unit cells.

The mean square deviation was computed, where the angle brackets denote the ensemble average over all possible time origins t and R is the location of the molecule center of mass or atom. A graph of this function versus values of the time interval \Delta t is shown in Figure 11-5 for CO2 and Figure 11-6 for He. The slopes of these curves at large time intervals can be used to determine the diffusion constant, using

(1)

The resulting diffusion constants are given in Table 11-1 along with values for two experimental systems [4].

System              Temperature    Diffusant     Diffusion Constant (cm2/sec)
============        ===========    =========     ============================
Expt. LDPE             298 K       He            6.8 x 10^-6
(0.914 g/cc)           298 K       CO2           0.372 x 10^-6
------------        -----------    ---------     ----------------------------
Expt. HDPE             298 K       He            3.07 x 10^-6
(0.964 g/cc)           298 K       CO2           0.12 x 10^-6
------------        -----------    ---------     ----------------------------
Simulated              350 K       He            1.00 x 10^-7
(0.87 g/cc)            350 K       CO2           1.09 x 10^-8

Table 11-1. Experimental and simulated diffusion constants.

The simulated diffusion constants are approximately an order of magnitude lower than the experimental values, even for more dense systems.

Nevertheless, we did find that He diffusion is about 10 times faster than CO2 diffusion, which is reasonable given the 30-fold difference in diffusion constants for a density of 0.964 g/cc, the 20-fold difference for a density of 0.914 g/cc, and our density of 0.87 g/cc.

11.4. Conclusions

The use of the efficient, parallel molecular dynamics code on smaller systems of only a few thousand atoms, but simulated for long time periods of hundreds of ps was demonstrated with this calculation.

Multiple polymer chains produced a significantly more disordered, amorphous structure than the original single chain.

We computed a ratio of diffusion constants between CO2 and He that is of the correct order of magnitude. The low absolute values of the diffusion constants are likely related to trapping within the polymer medium. As the simulation length increases, the mean square deviation will be primarily determined by the rate of hopping between regions, rather than the rate of motion within those (presumably interstitial) regions as is currently the case.

Further investigation of the orientation of the CO2 molecule with respect to its trajectory during the diffusion process may provide greater insight into the atomic-level processes occurring in this system. Such details can only be recovered from atomistic simulations.