Molecular Dynamics Study of Carbon Dioxide Sorption and

Article pubs.acs.org/Macromolecules

Molecular Dynamics Study of Carbon Dioxide Sorption and Plasticization at the Interface of a Glassy Polymer Membrane Sylvie Neyertz* and David Brown LEPMI (LMOPS), UMR 5279 CNRS, Grenoble INP, University of Savoie, University J. Fourier, Bât. IUT, Savoie Technolac, 73376 Le Bourget du Lac Cedex, France ABSTRACT: A series of large-scale molecular dynamics (MD) simulations of CO2 transport in a fully atomistic ∼50000-atom fluorinated 6FDA−6FpDA polyimide membrane were carried out under six different conditions of applied external gas pressure in order to assess the uptake mechanisms of the penetrant at the glassy polymer interface. CO2-induced volume dilation was found to occur immediately with a progressive shift of the interface location toward larger values associated with a widening of the interface and a decrease in density. However, it only leads to limited structural relaxations of the matrices and as such, the corresponding increase in mobility is small. Void-space and excess chemical potential analyses show that the holes within the matrix are getting larger as sorption proceeds, but that they are immediately occupied by penetrant molecules. The only way to accommodate more penetrant is to further swell the polymer matrix. The concave behavior of the sorption curves is well reproduced by the models and compares favorably with experimental data, except for the highest-pressure system which is in a different regime and exhibits a quasi-supercritical behavior for CO2. In all cases, the penetrants first undergo a rapid adsorption at the polymer surface. This is followed by a slower uptake mode, with the time-dependent diffusion coefficient being only accessible for short time-intervals because of the moving boundaries. All trajectories display the oscillations in voids and occasional jumps mechanism, and no transitions to the paths characteristic of liquid-like diffusion were seen, even in the most swollen systems.

1. INTRODUCTION The transport of gases in dense polymer membranes is commonly described by the sorption−diffusion mechanism,1−3 in which gas molecules in an upstream compartment enter the polymer matrix, diffuse across it and finally desorb onto a downstream gas compartment. Permeability P is the rate of transport for the penetrant through the membrane and is usually expressed as the product of an average solubility coefficient S and an average diffusion coefficient D. In the case of glassy polymers, transport of small penetrants such as O2 and N2 is not known to lead to any obvious long-term changes in the matrix.3 On the other hand, penetrant-induced plasticization phenomena including volume-swelling, higher gas permeabilities and hysteretic changes of the glassy matrix have been widely reported following exposure at medium to high CO2 activities, both in experimental studies4−12 and in molecular dynamics (MD) simulations.13,14 In glassy polymer membranes, experimental sorption isotherms for light gases are usually found to be concave to the pressure axis. This is the case for fluorinated polyimide membranes based on the 4,4′-(hexafluoroisopropylidene)diphtalic dianhydride (6FDA), which are interesting materials for gas separation applications because of their good mechanical, chemical and thermal properties combined with high permeabilities and permselectivities.15 We have earlier carried out extensive fully atomistic molecular dynamics simulations of CO2 sorption and desorption in bulk 6FDA− © XXXX American Chemical Society

6FpDA, 6FDA−6FmDA and 6FDA−DAM glassy matrices.13,14,16 The polyimide matrices were first loaded with CO2 in increments of 2% up to ∼30%, i.e., up to nominal CO2 concentrations of ∼200 cm3(STP) cm−3, in order to mimick the experimental procedure of progressive loading. Each sorption phase was then followed by a progressive desorption phase in decrements of 2%. These bulk models were able to reproduce experimental penetrant-induced plasticization and sorption−desorption hysteresis on the time scale available to MD simulations. This work also led to the development of the “Trajectory-Extending Kinetic Monte Carlo” (TEKMC) method,16 which extended the time scale of the penetrants mean-square displacements (MSD) by at least 3 orders of magnitude at negligible computing costs in order to reach the Fickian regime for diffusion. The model penetrant self-diffusion coefficients Dgas compared very well to experimental evidence.14 However, such polymer+gas bulk models can only be representative of the membrane core and thus a sophisticated iterative technique had to be used in order to estimate the vapor pressure of CO2 that would have to be applied to obtain the imposed uptake.13,17 In order to assess the uptake mechanisms at the gas/polymer interface, true membrane models with explicit surfaces will have Received: October 2, 2012 Revised: January 23, 2013

A

dx.doi.org/10.1021/ma302073u | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

used for all unlike-atom LJ interactions. Electrostatic interactions were evaluated using the Ewald summation method.41,42 The truncation radii used for both real-space electrostatic and van der Waals contributions were set to 12 Å, and (α,Kmax) parameters were chosen to obtain optimal convergence of the Ewald sums.43 Long-range corrections to the van der Waals contributions to the energy and the pressure were added using the approximation that the radial distribution functions are equal to unity beyond the cutoff.40 The time-step in the integration algorithm was set to 10−15 s. The temperature T was maintained at the required value, i.e., 308 K in the production runs, by loose-coupling to a heat bath44 with a constant equal to 0.1 ps. The pressure P was also maintained by loose-coupling with a constant equal to 5 ps.45 2.2. Preparation of the 6FDA−6FpDA Membrane. Each 6FDA−6FpDA chain had a length of 50 monomers (3302 atoms) and the model membrane contained 15 molecules, i.e. a total of 49530 atoms. Such a size is larger by a factor 3−10 than the typical MD boxes used for fully atomistic bulk systems,46,47 which makes it very expensive from a computational point-ofview. The preparation procedure of our former ODPA−ODA polyimide membrane30−33 had been directly inspired by the experimental solvent-casting process. The solvation procedure was mimicked by creating a polymer box with the density of an ideal 10% solution and a diamond-lattice impenetrable wall was subsequently added on either side of this virtual “solvated” polymer. The process of densification resulting from solvent evaporation was modeled by compressing in an affine manner the wall + polymer system in one direction. Following hightemperature relaxation and cooling down of the system, the wall was pushed away and the free-standing membrane was allowed to relax on its own. As noted before,32,33 it would have been better to dissolve the polymer in an explicit solvent and subsequently evaporate the solvent, but this would have been prohibitive in terms of computational resources. Moreover, the rigid polyimide chains could hardly have adapted to an everchanging concentration within the limited MD time scale. In the present work, we simplified the procedure by starting directly from a bulk melt rather than from a virtual polymer solution and placing a smooth flat repulsive wall (instead of a diamond-lattice one) at a distance of ∼100 Å on either side of the polymer. In practice, these changes have little effect since the glassy polyimide membranes exhibit characteristics which are similar to the former ODPA−ODA systems. A 6FDA−6FpDA bulk model in the melt was prepared using the well-documented hybrid pivot Monte Carlo−molecular dynamics (PMC−MD) single-chain sampling procedure,31,35,46,48 in which MC pivot moves of rotatable torsions are regularly attempted and the change in local energy is evaluated based on Flory’s hypothesis.49 The strength of this technique is that it has been systematically validated for a large variety of chain lengths and structures, including 6FDA− 6FpDA,35 by comparing the PMC−MD results with those of bulk melts decorrelated using MD on its own. Furthermore, the model densities of 6FDA−6FpDA models at room temperature were found to be within ∼1% of the experimental values, which validated both the force-field used and the preparation procedure.35 The melt was prepared at 650 K since the Tg of 6FDA−6FpDA has been reported as being in the range 575− 605 K.35 The ∼(85 Å)3 bulk model was run under constantvolume NVT conditions until its thermodynamic properties stabilized, i.e., for ∼300 ps.

to be created. A common approach is to extend one of the simulation box axes in order to eliminate interactions of the parent chains with their images in one direction,18−22 but this can only be used for molecules with enough mobility to relax from bulk to surface-like chains within the limited time scale of MD simulations. This is clearly not the case for glassy long polymers chains. It is also possible to use multiscale approaches with coarse-grained models,23−28 but such an approach requires intricate parametrization and reverse mapping procedures. In the present work, we apply instead an atomistic wallcompression technique which has already proven efficient for the 40-mers ODPA−ODA polyimide.29−33 A similar procedure is used here to create a fully atomistic membrane model of the 50-mers 6FDA−6FpDA (Figure 1) with ∼50000 atoms and

Figure 1. The 6FDA−6FpDA polyimide.

dimensions of ∼95 × 85 × 85 Å 3 . Six different CO 2 concentrations were then placed in the reservoirs on each side of the membrane. The basic mechanisms underlying gas adsorption, penetration and diffusion through the six models were followed with MD simulations over a total of 10000− 20000 ps. In addition, the different penetrant concentrations allowed us to study the effect of high-pressure penetrants and the associated plasticization phenomena. Results provide a picture of the early stages of CO2 sorption in such glassy matrices and have been compared to those obtained in the bulk.13,14,16 The preparation and characterization of the fluorinated membrane model is described in section 2. The insertion procedure for CO2 is given in section 3, while analyses relating to CO2 sorption, diffusion and induced-plasticization are discussed in section 4. All calculations were performed using the gmq package34 in its parallel form on the French GENCI supercomputing centers as well as on the MUST Linux server at the University of Savoie.

2. THE 6FDA−6FPDA POLYIMIDE MODEL MEMBRANE 2.1. Force-Field and Simulation Parameters. The forcefields used for the fluorinated matrices and the CO2 penetrants have already been described in detail.13,14,16,35 They will thus only be briefly summarized here. The polyimide chains were described by a combination of angle-bending, torsional, out-ofplane, excluded-volume Lennard-Jones LJ 12−6 and electrostatic potentials. The last two, so-called “nonbonded” potentials, were applied to all atom pairs separated by more than two bonds on the same chain or belonging to different chains. CO2 was modeled as a rigid three-site molecule with each atom carrying a partial charge and LJ 12−6 parameters. In all cases, the bonds were kept rigid with the SHAKE algorithm,36 while a special vector constraint37 had to be used in order to maintain the penetrant O−C−O angle at 180° and thus avoid nonequipartition of the kinetic energy. The polyimide parameters were taken from the TRIPOS 5.2 forcefield,38 with the partial charges being obtained from ab initio calculations,35 while the CO2 parameters came from the Zhang and Duan model.39 Lorentz−Berthelot combining rules40 were B


Macromolecules

Article

Figure 2. (a) Symmetrized mass density distribution ⟨ρ(x)⟩ with a slab width of 1 Å for the 6FDA−6FpDA membrane as a function of the distance x from the membrane COM (dashed line with circles). The other line is a fit of ⟨ρ(x)⟩ to eq 2. (b) Symmetrized P2(cos θx) with a slab width of 5 Å for different rotatable angles as a function of the distance from the middle angle atom to the membrane COM (left axis). Also shown for comparison is the polymer ⟨ρ(x)⟩ (right axis).

The x-axis of the simulation box was then subsequently extended with respect to the origin set at the polymer center-ofmass (COM) and a smooth flat repulsive wall was placed perpendicular to the x-axis on each side of the polymer COM at a distance xwall = ± 150 Å. All atoms i interacted with the wall if their absolute xi positions were greater than xwall according to the Uwall potential given in eq 1: Uwall(xi) =

1 k wall(|xi| − xwall)2 , 2

if |xi| > xwall

conditions. Following a relaxation of 500 ps at 308 K under NVT conditions, the system was switched for the next 500 ps to conditions in which the transversal Pyy and Pzz components of the pressure tensor, P, were maintained close to 1 bar by loosecoupling to the y and z box lengths, respectively. This allows the surface area of the membrane to relax in response to changes in the transversal pressure while keeping the MD box orthorhombic. The x box length thus remained the same, with the membrane in the middle and an empty reservoir on either side. Such conditions allow for the polymer matrix to swell, or contract, in all three directions. The Ewald sum was found to converge41,42 with α = 0.20 Å−1 and Kmax = 12. Once free to move, the polymer chains relaxed slightly in the x-direction and the box dimensions changed by a maximum of 1.2 Å in both y and z directions. The final configurations were used as the starting points for the polymer +CO2 models. 2.3. Characterization of the Free-Standing 6FDA− 6FpDA Membrane. The 6FDA−6FpDA slab mass density distribution as a function of the distance along the x-axis to the membrane COM, ρ(x), is given in Figure 2a. Please note that we have arbitrarily chosen to display the results with respect to the left side of the membrane (x < 0, the COM being at x = 0) in order to be consistent with the representation of the CO2 flux moving from the left to the right of the membrane. However, results are systematically symmetrized over both the left and the right sides of the membranes. The slight density peak which appears at the interface is reminiscent of the compression procedure, but is not so different from the membrane core to resemble peaks in the vicinity of flat interfaces.29 Indeed, ρ(x) exhibits the sigmoid̈ al profile characteristic of free-standing interfaces and can be fitted to the following hyperbolic equation (eq 2):50,51

(1)

with kwall being a force constant arbitrarily set to 50 kg s−2. Such a harmonic Uwall gives a smoothly decaying soft repulsive wall, and in addition, the forces are limited by an upper bound which corresponds to a maximum displacement in one time step of 0.1 Å. The position of walls can be controlled, and in the present work, they were progressively drawn closer to the polymer melt at 650 K and at a very slow rate of 0.001 Å ps−1, until the middle of the polymer regained a density close to that of the bulk. This took almost ∼100000 ps, and during this phase, the nonbonded potentials were replaced by the much faster purely repulsive form of the excluded-volume, the Weeks−Chandler−Andersen (WCA) potential, which has been shown to be a good approximation when polyimide configurations are considered.31 In addition, the temperature coupling constant was set to 0.01 ps in order to avoid heating the system too much. Once the polymer was compressed, it was cooled down to 308 K at a rate of −1 K ps−1 and run for an additional 2000 ps with the wall still on. The full nonbonded LJ and electrostatic potentials were then reintroduced and the wall was pushed back to xwall = ± 150 Å. It should be noted that 2D (rather than 3D) periodic boundary conditions are sometimes applied in MD simulations of polymer membranes by removing periodicity in one direction in order to eliminate interactions of the parent chains with their images in that direction.18,19,21,22 However, this is usually carried out in the absence of a wall and, as pointed out in the Introduction, it can only be used for small and flexible chains with enough mobility to relax on their own from bulk to surface-like chains within the limited MD time scale. This is definitely not the case for glassy dense long chains, and in our approach, the empty reservoirs are defined by the ∼100 Å distance between each side of the polymer matrix and its impenetrable wall. As such, there is no need to switch from 3D to 2D periodic boundary conditions, and all simulations reported here were carried out under 3D periodic boundary

ρ(x) = ρmiddle

1 − tanh[2(x − h)/w] 2

(2)

where ρmiddle is the density in the middle region of the film, h is the location of the interface with respect to the membrane COM and w is generally referred to as the “interfacial width” in this structural-based definition. The parameters fitted to the 6FDA−6FpDA data (line in Figure 2) are ρmiddle = 1482 kg m−3, i.e., ∼0.3% from its average experimental density 1477 ± 3 kg m−3,35 h = ± 43.7 Å for the location of the interfaces and w = 2.7 Å for the interfacial width. It should be noted though that C


Macromolecules

Article

visualization show that the molecules are intertwined. In agreement with experimental characterizations,58,59 there are no signs of crystallinity either. As for ODPA−ODA,30−33 our 6FDA−6FpDA model is in agreement with surface-layer chains reported for much shorter, more flexible and often less-realistic polymers, in which smaller sizes and faster relaxation times allow for a gradual change from surface to bulk-like structures in the middle of the membrane models.19−22,60−72 This is also the case for coarse-grained relaxed structures28 and it confirms that our preparation technique gives reasonable configurations for glassy polyimide surface-layers.

the actual distance over which the polymer density rises from zero to values close to that at the center is 2 or 3 times the value of the fit parameter w, as can be seen clearly in Figure 2a. From this analysis, the width of this membrane model is h × 2 = 87.4 Å, although such a definition depends both on the property and on the way it is analyzed.31−33 The fluorinated interfacial width w is also quite small compared to e.g. the 4.6 Å found for ODPA−ODA.32 This could be related to the density of 6FDA−6FpDA being more than 100 kg m−3 higher than that of OPDA−ODA, and consequently fluorinated chains having more difficulties rearranging themselves at the interfaces. However, as for the location of the interface, the interfacial width is quite subjective. Indeed, the small polymer segment MSD, which are a dynamical property, display an increase by a factor of ∼2 in the vicinity of the interfaces (results not shown for the pure polymer, see section 4.5. for the polymer+CO2 systems). This would define a “dynamical” interfacial width of ∼20 Å, which is larger than the “structural” interfacial width described in eq 2 by w. The alignment of specific chain angles can be characterized by second-order Legendre polynomial functions, P2(cos θα), with θα being the angle between the α-axis and the vector between atoms i and k in a {i, j, k} triplet (eq 3): P2(cos θα) =

3 1 ⟨cos2 θα⟩ − 2 2

3. CO2 INSERTION Following the preparation of the pure 6FDA−6FpDA membrane, CO2 molecules were inserted into the empty spaces between the membrane surfaces and the soft repulsive walls situated at xwall = ± 150 Å. The size of a reservoir was ∼100 × 85 × 85 Å3. Since it is not possible to estimate a priori the amount of gas that will actually enter into the dense matrix, it was decided to carry out simulations with six different initial concentrations, i.e., 1000, 2000, 3000, 4000, 8000, and 12000 CO2 molecules per reservoir. Results obtained from both reservoirs of each membrane were symmetrized in order to improve the statistics. While the CO2 molecules entering the matrix from either side do ultimately interact with each other, the percentage of penetrants actually entering the central part remains quite low on the MD time scale and thus we are not in a regime where such interactions are likely to be significant. Although initiating the simulations with the CO2 being inserted into a single reservoir (rather than on both sides) would seem closer to the experimental permeation setup (in which molecules enter the matrix from an upstream compartment and desorb onto a downstream gas compartment),3 this leads in practice to a translation of the membrane in the x-direction, as it is not supported, and consequently to a constant change in the size of the reservoirs. Our setup is thus closer to sorption experiments, which are usually carried out from both sides of the membrane. In the remaining part of the paper, the different systems will be referred to by default as nini CO2, with nini being the number of CO2 molecules initially inserted in each reservoir of the starting configurations. It should also be noted that, in principle, instantaneous values of the pressure can be obtained at any virtual interface through the system by analyzing the momentum flux.73,74 However, in the case where Coulombic forces are present, this is not possible due to the long-range nature of the forces. Once in mechanical equilibrium, the Pxx component of the pressure tensor has to be the same everywhere. In order to estimate the time required for the initial shock waves created by the expanding gas to dissipate, tests were carried out on the same systems with the partial charges set to zero and short-range repulsive van der Waals forces. This happened typically on a time scale of about 100−200 ps. The six MD simulations were carried out under the same conditions of constant transversal pressure as described for the pure membrane, and at 308 K for up to 10000 ps for the systems with 1000, 2000, and 3000 CO2 per reservoir. Although such calculations take several months to run on an 8-processor machine, they were extended up to 20000 ps for the models that displayed a clear plastifying behavior, i.e. 4000, 8000, and 12000 CO2 per reservoir. The Ewald sums converged with α = 0.17 Å−1 and Kmax = 12, and all other parameters were

(3)

The limiting values of P2(cos θα) with respect to the normal to the interface are −1/2 for a perfectly perpendicular alignment, 1 for a perfectly parallel alignment and 0 for a random alignment of the vectors defining θα . The P2(cos θx) were calculated for “rotatable” angles, i.e., the Car−C−Car angles between both aromatic rings in the 6FDA and 6FpDA moieties, the CF3−C− CF3 angles in their middle and the C−N−C 6FDA−6FpDA linkages (Figure 1). The same analyses were carried out in the bulk models,35 and as expected from isotropic systems, all these P2(cos θx) were close to zero. On the other hand, Figure 2b shows that the membrane ring−ring Car−C−Car and C−N−C pivot angles tend toward −1/2 in the vicinity of the surface, which indicates a parallel alignment with respect to the surface (i.e., a perpendicular alignment with respect to its normal). On the other hand, the CF3−C−CF3 P2(cos θx) are positive, and consequently will have a tendency to align more perpendicular to the surface than in the bulk. Since the CF3−C−CF3 angles are already perpendicular with respect to the direction of the main backbone, the Legendre analyses confirm that polyimide chains at the vicinity of the interface align in a parallel fashion with respect to the surface, a feature which is well-known in free-standing membranes.29−31,52−57 In addition, the average mean-square 6FDA−6FpDA chain radii of gyration ⟨S2⟩ as a function of x were calculated and, despite the limited statistics on the number of chains, their perpendicular components, ⟨S⊥2⟩, were less than 20% at the surfaces. This indicates a flattening of the chains, which are sometimes referred to as “pancake”-like structures in the literature.29−31,54,55 As found before, both the alignment and flattening persist to a certain point into the membrane cores, albeit attenuated (Figure 2b). This is clearly related to the compression step, the low mobility, and the relatively small length scale if compared to experimental sizes. Consequently, the chains in such models can be considered as closer to surface-like than to bulk-like structures, which is the reason why bulk models should be studied separately.13,14,35 It is also important to note that the chains do not really lead to separate layers, and that their D


Macromolecules

Article

while the location of the membrane interface h is shifted toward larger values with time. Along with the deformations in the yand z-directions, this leads to consequent volume dilations (Figure 4d). Interestingly, the evolution of the interface location h in the x-direction seems to be almost independent of the initial number of CO2 in the reservoir (Figure 4b). On the other hand, the swelling in the y- and z-directions, i.e., the swelling of the surface area, is directly correlated to the number of CO2, as shown by the subsequent %volume swelling (Figure 4d). There is no real preference for either y-swelling or zswelling since differences between both amount to less than 1%. In terms of volume swelling, 10000 ps corresponds to 3.3% in 1000 CO2, 6.9% in 2000 CO2, 7.1% in 3000 CO2, 7.8% in 4000 CO2 and 8.9% in 8000 CO2. 20000 ps correspond to 12.1% in 4000 CO2 and 16.9% in 8000 CO2. Such values are fully consistent with the range of swelling studied in the bulk models, i.e. up to 16%,13 and with the experimentally estimated 20% volume dilation of 6FDA−6pDA following exposure to CO2 at 60 atm.8 The swelling of the membrane consequently leads to lower ρmiddle values as displayed by Figure 4a and a progressive “opening” of the glassy structures. A closer comparison of the swelling effects in either the longitudinal (along x) or transversal (along y or z) directions can be obtained by calculating the respective strains defined as:

the same as described in section 2. Configurations were stored at 20 ps intervals, and thermodynamic and conformational data every 1 ps. In addition, a series of short 100 ps-simulations were carried out with the configurations being stored at shorter 1 ps intervals in order to better follow the initial adsorption phase. Schematic representations of some systems under study are displayed in Figure 3.

Figure 3. The polymer and the CO2 penetrants displayed along the xdirection in space-filling models with the VMD 1.8.2 software.75 The color code is as follows: for the polyimide, cyan = C, red = O, blue = N, white = H, green = F; for the gas, yellow = C, red = O. Full boxes of length ∼300 Å in the x-direction are shown for (a) the 1000 CO2 system at 10000 ps and (b) the 12000 CO2 system at 20000 ps. (c and d) ∼25 Å2 close-ups of part a in the vicinity of an interface and in the center of the membrane, respectively.

% longitudinal strain =

% transversal strain =

(V (t ) − V0) × 100 V0

( A (t ) − A0

A0 )

(5a)

× 100 (5b)

where h(t) and A(t) are the location of the interface h and the surface area A at time t, respectively, and h0 and A0 correspond to those of the pure polymer matrix. The time-evolution of these strains is shown in Figure 4e for the longitudinal (eq 5a) and in Figure 4f for the transversal (eq 5b) components. The longitudinal strain amounts to 8% at 20000 ps, but as seen in Figure 4e and as noted previously in Figure 4b, it is very similar for all five lower CO2 concentrations, in spite of the increased uptake of CO2 into the membrane that accompanies the increase in CO2 content (see later). It suggests that the tendency to swell is partly counteracted by the increase in the gas pressure in the reservoir, which compresses the membrane in the x-direction. In the case of the transversal strain, the membrane does not only expand because of the swelling response of the matrix, but also because of the increased tension in the film. At constant applied transversal pressure, this increase in the film tension counteracts the increased pressure in the reservoir due to the gas. Under the imposed model conditions, the transversal strain is always smaller than the longitudinal one; i.e., it only amounts to 4% at 20000 ps for the 8000 CO2 system, which means that the polymer finds it easier to swell along x in spite of the reservoir pressure. Biaxial mechanical tests were carried out on the pure membrane model, i.e., in the absence of CO2, at very high applied tensions (∼600 bar) and they only resulted in deformations of ∼1% (not displayed). This indicates that the major part of transversal strain is rather due to the swelling effect. In spite of similar h, the width of the interface in the x-direction w is also correlated to the reservoirs with the polymer mass density profiles becoming smoother and w larger when CO2 pressure increases. This local interface relaxation is mostly seen at the early stages of the simulations when the gas penetrates very fast into the

4. CO2 PERMEATION 4.1. Volume Swelling. CO2-induced volume dilation is well-known to occur in many glassy polymers, including 6FDAbased polyimides.10 As noted for bulk models,13 the gasinduced volume swelling at time t can be directly measured in molecular simulations by using eq 4: % volume swelling =

(h(t ) − h0) × 100 h0

(4)

where V(t) is the volume of the membrane at time t and V0 is the volume of the pure polymer matrix, i.e., before any gas has been absorbed. In isotropic bulk models with periodic boundary conditions, a polymer expansion/contraction along an axis simply leads to an increase/decrease of the corresponding box length. This is also the case here for the membrane model in the y- and z-directions. However, while it can similarly expand or contract along the x-direction, this will slightly affect the size of the reservoirs but the x-axis box length remains constant (Figure 3). As such, the dimensions of the membranes along the x-direction at time t were estimated from the polymer symmetrized slab mass density profiles ⟨ρ(x)⟩ at time t by using the same hyperbolic form (eq 2) as for pure 6FDA−6FpDA (Figure 2a). The time evolution of the density ρmiddle in the middle region of the film, the location h of the interface with respect to the membrane COM and the interfacial width w are shown in Figure 4, parts a−c. The resulting % volume swelling of the membranes (eq 4) calculated from the product of 2h with the box lengths along the y- and z-directions is displayed in Figure 4d. It is clear from Figure 4 that for the five lower CO2 systems, the density of the middle region of the film ρmiddle decreases E


Macromolecules

Article

Figure 4. Time evolution of the (a) ρmiddle, (b) h, and (c) w parameters extracted using eq 2 from the time-dependent polymer along the xdirection. (d) Corresponding membrane total % volume swelling. (e) Time-evolution of the longitudinal strain (eq 5a) and (f) that of the transversal strain (eq 5b) for all models under study except for the 12000 CO2 for which the results are on a different scale (see text). The legends in parts b−f are the same as that in part a. It refers to the numbers of CO2 molecules initially placed in each reservoir.

supercritical CO2 (scCO2) on polyimide membranes.76,77 An interesting result is that, at pressures above the supercritical point, CO2 permeability decreases in un-cross-linked polymer membranes, thus pointing toward a major structural reorganization toward a new pseudoequilibrium state with reduced free volume.76 On the other hand, scCO2-conditioned samples show higher sorption capacities than unconditioned films. It has been suggested by Koros et al. that these different structural responses might be related to the concentration gradient of the penetrant in the polymer being different in permeation and sorption experiments.76 Although the MD time scale is much too short for any large-scale polymer configurational changes to occur at the temperature under study, the setup of our model is here closer to that of a sorption experiment, and as such, we expect to see an increase in sorption capacity rather than a decrease in free volume. According to Wessling et al.,4 CO2-induced dilation on related 6FDA-based polyimides starts almost immediately on the experimental time scale. This is supported here by the fact that it occurs in all six models under study, even on the very limited MD time scale. However, such simulations are

matrix (see later) and w tends to remain constant afterward (Figure 4c). The highest-pressure 12000 CO2 system exhibits a very different behavior from the others. The reservoir concentrations are so high that the CO2 molecules first force the polymer membrane to contract in the x-direction before allowing it to come back to its initial dimension when the gas penetrates into the membrane (Figure 4b and the associated longitudinal strain obtained from eq 5a, which is ∼0.2% at 20000 ps). While they still prevent the swelling of the membrane in the x-direction, there is an important increase in the surface area (over 20% swelling in both y and z directions) and in the width of the interface. The amount of relaxation is so large that the cohesion of the glassy membrane is much more reduced than in the other systems (Figure 4a), the volume expansion amounts to 38.5% at 10000 ps and to 51.9% at 20000 ps, while the transversal strain (eq 5b) reaches 23% at 20000 ps. Such values are clearly typical of a different regime for this membrane + CO2 system, and the huge volume increase suggests a quasisupercritical behavior for CO2. To our knowledge, there have only been very few experimental studies on the effect of F


Macromolecules

Article

Figure 5. (a) CO2 concentrations C(t) in the membrane models as a function of the square root of time. The 12000 CO2 C(t) is displayed on a different scale (right axis) in order to better visualize the other curves (left axis). (b) Comparison of the sorption (left axis) and dilation (right axis) kinetics for 8000 CO2 (c) C(t) in the membrane as a function of reservoir pressure P(t) at t = 10000 ps (symbols). The line is a fit of the data to the dual-mode sorption isotherm model (eq 8). (d) S(t) as a function of reservoir pressure P(t) at t = 10000 ps (symbols). The line shows the experimental data for 6FDA−6FpDA taken from Coleman and Koros.6

Figure 5a displays the time-dependent CO2 concentration C(t) in the membrane (from both reservoirs) calculated using eq 7 as a function of the square root of time, t1/2, for the six 6FDA− 6FpDA + CO2 models. In the associated 6FDA−6FpDA bulk models,6,13,14 the nominal CO2 concentration range was 0−200 cm3(STP) cm−3, which corresponded roughly to the 0−60 atm CO2 pressure range studied experimentally by Coleman and Koros.6,8 The order of magnitude of the C(t) is similar here with maxima at 20000 ps of ∼160 and ∼250 cm3(STP) cm−3 for the 8000 CO2 and the 12000 CO2 systems, respectively. As found for O2 sorption in polyimides,33 the gas C(t) uptake curves and the underlying nCO2(t) curves (not displayed) for the 1000, 2000, 3000, 4000, and 8000 CO2 systems (left axis of Figure 5a) first exhibit an initial jump due to the penetrants undergoing a rapid adsorption and saturating the polymer interfacial regions, while establishing a quasi dynamic equilibrium with the gas in the reservoirs. This is followed by a slower uptake rate due to the gas molecules at the interfaces accessing progressively the denser parts of the matrix, which is characterized by a more linear dependence with respect to t1/2. It is often assumed that a linear relationship between the weight gain of the polymer and t1/2 is the signature of Fickian diffusion.82,83 However, at such short times, the nCO2(t) and C(t) curves also appear to be linear with respect to t and the molecular time scale is clearly too limited to attain Fickian diffusion. Nevertheless, we still favor here the representation of the uptake curves vs t1/2, as this is usually used in experimental mass transport kinetics.4,78 In the case of the 12000 CO2 system, the initial jump is missing in the nCO2(t) curve, which is linear with respect to t1/2 from the start.

obviously not able to access possible slower penetrant-induced volume relaxations.4,12,78−80 4.2. Gas Uptake Curves. The amount of CO2 that has entered a polymer matrix at time t is also clearly correlated to this time-dependent evolution of the membrane width. As noted before,31−33 the definition of the membrane thickness for characterizing sorption is not only a source of difficulty in molecular models but also in experimental measurements.4,81 In the present case, the width of each membrane model at time t was once again obtained from the time-dependent location of the interface h parameter (eq 2 and Figure 4b). The number of CO2 molecules located inside the membrane at time t, nCO2(t), could then be recorded using 2h(t) as the thickness of the polymer. Experimentally nCO2(t) is often expressed as an equivalent volume of gas VSTP g (t), i.e. it is the volume that the nCO2(t) molecules would occupy in the pure state if held at standard temperature and pressure (STP: 273.15 K; 1.013 × 105 Pa) and if they behaved as an ideal gas:33 V gSTP(t ) =

nCO2(t )kBT STP P STP

(6)

The penetrant concentration at time t, C(t), expressed in cm3(STP)/(cm3 polymer) which is usually referred to for simplicity as cm3(STP) cm−3, can then be obtained from the equivalent volume of the gas, VSTP g (t), divided by the volume of the polymer at time t, Vpol(t) (Figure 4d): C(t ) =

V gSTP(t ) Vpol(t )

=

nCO2(t )kBT STP P STPVpol(t )

(7) G


Macromolecules

Article

This suggests that there is a facilitated immediate flow of gas molecules from the interfaces to the denser parts of the matrix. Along with the nonlinear Vpol(t), it leads to a nonlinear C(t) curve (right axis of Figure 5a), thus confirming that CO2 is here in a different regime than in the other systems. When comparing CO2 sorption and dilation kinetics in 6FDA-based polyimides, Wessling et al. reported a small delay at low pressures between the gas uptake and the polymer strip length increase.4 The C(t) uptake curve for the 8000 CO2 system is displayed as a function of t1/2 along with the associated volume swelling in Figure 5b. Similar plots were obtained for the other systems, and none of them show any delay between the sorption and the dilation kinetics on the MD time scale. Wessling et al.4 explained the apparent delay by the fact that dilation occurs simultaneously with the CO2 increase in the middle of the film in their experimental setup, but that at earlier times, it is a complex combination of swelling in the length direction and in the thickness of the film. Indeed, as also pointed out by Wessling et al.,4 any CO2 entering does appear to swell the matrix (Figure 5b). The effective volume dilation was a bit more difficult to see in some of the 6FDA−6FpDA bulk models, since it was found to be less than 1% for concentrations under ∼40−60 cm3(STP) cm−3.6,13,14 Interestingly, this CO2 concentration range can be considered experimentally6,9 as the start of plasticization, i.e., it is the point where the permeability passes through a minimum before the increasing diffusivity starts compensating for the decreasing solubility of the penetrant.9,10,78 However, volume dilation starts before that, and in our membrane models, all C(t) are above 40−60 cm3(STP) cm−3 at 10000 ps (Figure 5a), which results in swelling well over 1%. All six 6FDA−6FpDA + CO2 models can thus be considered as being plasticized by the penetrant. CO2-sorption into our model membrane has another effect, i.e., the gas density in the reservoirs changes constantly with time. Of course, it is not possible to ever reach steady-state within the MD time scale (experimentally, steady-state permeabilities at conditioning pressures were obtained over 2−3 weeks),6 but it is fairly reasonable to assume that the reservoir pressure at the later stages will be closer to the equilibrium pressure than at the very start. In order to estimate the time-dependent reservoir pressures, the average reservoir gas densities at time t were measured in the |x| = 60−80 Å regions of each model, i.e. where the gas density profiles are flat. A series of constant-volume NVT simulations of pure 512molecule CO2 systems were then carried out at 308 K in the range of the various reservoir gas densities. The corresponding average pressures were used to obtain a pressure vs density calibration curve, from which the pressures in the gas reservoirs at time t, P(t), could be estimated using the average reservoir gas densities in the |x| = 60−80 Å regions. Initial pressures are difficult to extract directly due to the inhomogeneity of the gas density at the start. By the time the gas density has started to stabilize (∼10−20 ps, see also section 4.6), some of the penetrant molecules have already adsorbed onto the surface and have started to migrate into the interior of the membrane. Nevertheless average densities in the gas reservoir region from the 60 < |x| < 80 Å over the period from 20 to 100 ps have been extracted and converted into equivalent pressures based on the calibration curve. The values obtained were ∼36, ∼63, ∼80, ∼87, ∼130, and ∼760 bar for the 1000, 2000, 3000, 4000, 8000, and 12000 CO2 systems, respectively. It should be noted that at 308 K, i.e., just above the critical

temperature of CO2, the pressure vs density curve of CO2 is highly nonlinear.39 The pressure is relatively insensitive to the density in the range from 200 to 600 kg m−3 but rises rapidly above densities of ∼800 kg m−3. For this reason, the initial pressure of the 12000 CO2 per reservoir system is particularly high. This high sensitivity to the density also means that the pressure in this system falls rapidly as the CO2 penetrates the membrane and the system expands transversally. By 20 ns, its pressure is even marginally lower than that of the 8000 CO2 per reservoir system. In agreement with the experimental behavior for such glassy polymers,3 the C(t) are concave to the pressure axis P(t) (Figure 5c). Experimentally obtained uptake curves are very often fitted to the “dual-mode” (DMS) sorption model:84 C(t ) = CD + CH = kDP(t ) + CH′

bP(t ) (1 + bP(t ))

(8)

where the subscript D refers to Henry’s law mode sorption and the subscript H to Langmuir mode sorption. Although experimental data are obtained over much longer time scales and also in a slightly different way by fixing the external gas pressure, rather than the initial amount of gas, it is nevertheless interesting to make comparisons. In the idealized DMS model, kD is Henry’s solubility coefficient and characterizes sorption into the densified equilibrium matrix of the glassy polymer. Langmuir sorption is related to the nonequilibrium excess volume associated with the glassy state and is described by two parameters, the sorption capacity C′H and the affinity parameter b. As seen in Figure 5c, the simulated C(t) as a function of P(t) estimated at t = 10000 ps do fit very well to the form of eq 8, with the exception once again of the 12000 CO2 system. However, although the idealized DMS model does correlate with various factors related to gas transport,3,6,8 it is also known to suffer from some fundamental shortcomings. Experimentally reported parameters for the same polymer show often scatter,85 which could be due in part to the fact that they depend on sample processing conditions,86,87 and the “constant” kD, C′H and b parameters are not constant, but vary systematically with the pressure range used.88 In addition, molecular simulations have already shown that the penetrant uptake curves can be fitted to the form of eq 8 without any evidence of two different environments for the penetrants at the molecular level.13,33 As summarized in Chapter 1 of ref 3, the DMS model is an empirical, phenomenological model that fits the experimentally measured isotherms but it still lacks a molecular foundation and a predictive ability. It is clear from Figure 5c that it does fit very well model isotherms too. From the C(t) and the P(t) values at time t, it is possible to define a time-dependent solubility coefficient at time t, S(t) = C(t)/P(t). The S(t) at t = 10000 ps corresponding to the data of Figure 5c are displayed in Figure 5d. They are shown along with the experimental data obtained at 308 K in the 0−60 atm range by Coleman and Koros,6 which are in good agreement with those of other groups.89,90 Although the experimental data have been obtained with steady-state pressures and our simulations are not in equilibrium, it is reassuring to note that the solubility range exhibited by our models is rather close to that of the real systems. 4.3. Void Space. The polyimide void space as a function of CO2 concentration was analyzed using a simple geometric technique similar to the phantom sphere approaches.91,92 The percentage of probe accessible volume (PAV) is obtained by H


Macromolecules

Article

Figure 6. Symmetrized percentages of CO2-excluded probe accessible volume (%PAV) obtained with a probe of radius 1.8 Å as a function of x and at different times t for the 0, 2000, 4000, 8000 (a) and 12000 (b) CO2 systems. The black triangles in part b are the only example shown here of CO2included % PAV.

Figure 7. Symmetrized ⟨μex(x)⟩ for the gas in the polymer at t = 10000 ps, along with those for the pure membrane. The dashed lines indicate roughly the interface regions. Insertion energies were calculated by either taking into account the polymer and all other CO2 penetrants (a) or the polymer on its own (b). The legend is the same on both graphs and the standard errors are displayed. Most of them in part b are smaller than the size of the symbols.

are all smaller than in the pure membrane. This indicates that all empty holes are already occupied by pre-existing gas molecules at t = 10000 or 20000 ps. The only way to accommodate more CO2 is to further swell the polymer matrix. 4.4. Excess Chemical Potentials. In molecular simulations, a test-particle insertion (TPI) method95 is often used to calculate the excess chemical potentials of gas molecules in a polymer matrix. A probe CO2 molecule is repeatedly inserted at random positions in the system and the changes in potential energy associated with the virtual insertion, ΔΦ, are used to obtain the probability of insertion, which is known to be related to the excess chemical potential of the gas in the polymer, μex by:13,31−33,47

repeated trial insertions of virtual spherical probes of preset radius (here 1.8 Å, i.e., half the Chung diameter of CO2,3 as in the bulk models13) into a membrane without any preassumption on the actual form of the holes. An insertion is “accepted” when the center of the probe does not overlap with the other atoms which are represented by hard spheres with standard van der Waals radii.93,94 The % PAV is then simply the fraction of “accepted” insertions with respect to the total number of trials. In a polymer + gas system, the %PAV can be calculated in two ways, i.e., by including or excluding the other gas molecules already present. For our membrane models, the symmetrized % PAV were further resolved as a function of x for slabs of width 5 Å and averaged over the last 1000 ps of either the 10000 or 20000 ps runs. They are shown in Figure 6 for some of the 6FDA−6FpDA + CO2 systems. In Figure 6, the reference is the low % PAV in the dense part of the pure membrane (squares with diagonal lines). Over the simulations, the CO2-excluded PAV increases both with reservoir pressure and time (Figure 6a), thus showing that the holes within the polymer matrix are getting larger. There is a progressive void-space swelling as more CO2 penetrates into the membrane, starting from the low-density interfacial regions and slowly getting toward the COMs. In the case of the 12000 CO2 system (Figure 6b), the matrix opening is much larger and more uniform over the whole system, in agreement with a supercritical behavior for the penetrant. It should be noted that the CO2-included %PAV (one example is shown in Figure 6b)

⎛ ΔΦ ⎞ μex = − kBT ln exp⎜ − ⎟ ⎝ kBT ⎠

(9)

with kB being the Boltzmann constant and T the temperature. μex can be averaged over several configurations of the production run, thus leading to ⟨μex⟩, which can itself be resolved as a function of x, ⟨μex(x)⟩. The ⟨μex(x)⟩ of CO2 in 6FDA−6FpDA were analyzed as a function of time t for slabs of width 5 Å using 10 million trials per system. The other CO2 present in the systems were either taken into account (Figure 7a) or excluded (Figure 7b) . For the pure polymer (Figure 7), the ⟨μex(x)⟩ range from the expected value of zero in the empty reservoirs to more negative I


Macromolecules

Article

Figure 8. (a) Spatially resolved Boltzmann-weighted probability density distributions in the 4000 CO2 system at t = 10000 ps for changes in potential energy upon repeated insertions of a CO2 probe. The slabs shown are situated in the reservoir far from the interface (|x| = 60−65 Å), in the reservoir close to the interface (|x| = 50−55 Å), in the interfacial region (|x| = 45−50 Å) and in the core (|x| = 0−5 Å). (b) raw ρ(ΔΦ/kBT) for the interface slab (|x| = 50−55 Å) and all systems at t = 10000 ps.

The slabs which are displayed in Figure 8a correspond to different parts of the membrane model. The energies associated with the peak in the reservoir far from the interface are close to zero since the CO2 probe interacts there only with other gas molecules. The weighted probability density distributions then get more displaced toward the left as |x| decreases, which is due to the progressive increase in matrix density toward the center of the membrane. It leads to more negative insertion energies and consequently more favorable sites for probe insertions. However, as seen in Figure 8a, a drawback of the Boltzmannweighted representation is that the exponential term becomes very large as ΔΦ gets more negative, and leads to an amplification of the statistical noise in the data. In order to compare all systems, it is thus often easier to display the raw ρ(ΔΦ/kBT), as shown in Figure 8b for the interface slab (|x| = 50−55 Å) at t = 10000 ps. As the reservoir pressure increases, a larger concentration of CO2 molecules at the interface decreases ρ(ΔΦ/kBT), since more of the free volume is occupied. However, although accessible sites are harder to find, they are also energetically more favorable as the penetrant can now interact with both the matrix and more of the other penetrants. This explains the shift of the peaks toward the left as the pressure increases. The swelling contribution at the interface can be studied by carrying out the same analyses while excluding the other CO2 penetrants (not shown). Adsorption sites with lower energies are indeed created with respect to the pure matrix, but Figure 8b suggests that these additional sites are immediately occupied by CO2 and the pressure-effect remains dominant. Moreover, when the polymer swells too much (12000 CO2), the interactions start getting less favorable in spite of the other CO2 and the curve is slighly shifted to the right with respect to 8000 CO2. Such differences tend to level off in the dense cores of the membrane (not shown) as the various distributions are found roughly in the same energy range than that displayed for the 4000 CO2 system (Figure 8a). If the TPI analysis has been carried out while excluding the other CO2 molecules, the distributions show an increase in the number of available sites for adsorption which is proportional to the swelling of the matrix. On the other hand, if the other CO2 molecules are taken into account, the void-spaces are being filled up and very few sites are left for adsorption, thus confirming the results of the void-space analyses. In addition to the dual-mode sorption model (eq 8), several other models have been proposed to describe sorption in glassy

values within the membrane. Insertions in the denser parts are most favorable as CO2 maximizes the number of attractive interactions with the matrix. On the other hand, those in the interfacial region only interact with the polymer on one side. This contributes to the gradient in the total chemical potential, which is known to be the underlying driving force for diffusion.33,96−98 In the polymer + gas systems, the presence of other CO2 affects ⟨μex(x)⟩. It increases gradually as the increasing number of penetrants in the membrane lowers the solubility, an effect which is clearly related to both reservoir pressure and time. While Figure 7 is at t = 10000 ps, the ⟨μex(x)⟩ for the 4000 and 8000 CO2 systems at t = 20000 ps (not shown) resemble the 12000 CO2 curve in Figure 7a. As such, the difference in ⟨μex(x)⟩ between the reservoirs and the membrane diminishes. However, it does not mean that there is no net flux anymore since the latter is proportional to the total and not just to the excess chemical potential of the gas. Indeed, equal solubilities imply that the density gradient is the main driving force for diffusion as the ideal part of the chemical potential is related to the density of the gas. In the case of analyses carried out by excluding interactions of the probe with the other CO2 molecules, solubility is higher than in the pure matrix (Figure 7b). As expected, the spaces in the swollen matrices are even more favorable insertion sites for the penetrants. The exception is 12000 CO2, in which the voids get so large (Figure 6b) that the decreasing density leads to a diminution of the attractive interactions. It is well-known that there is a delicate interplay between densities and solubilities:99 a lower polymer density should lead to more space for the penetrants, i.e. higher solubility, but it also leads to a lower cohesive energy density, i.e. less solubility, which seems to be the dominant effect in the highest-pressure system. As for ⟨μex⟩, the changes in the potential energy associated with the virtual insertion of the CO2 probe can be resolved spatially as a function of x. Normalized probability density distributions for ΔΦ/kBT, ρ(ΔΦ/kBT), were calculated and symmetrized for all slabs of width 5 Å at different times t. In order to obtain distributions more representative of the likely sites of adsorption, the ρ(ΔΦ/kBT) were multiplied by exp(ΔΦ/kBT) to obtain Boltzmann-weighted probability density distributions.33,35 They are shown in Figure 8a for the 4000 CO2 system at t = 10000 ps. J


Macromolecules

Article

Figure 9. Trajectories along x for CO2 molecules in the (a) 4000 and (b) 12000 CO2 systems. The dashed lines indicate the time-dependent boundaries of the respective membranes. (c) CO2 trajectory displayed with VMD 1.8.2.75 Every yellow line segment spans the 20 ps storage timeinterval and the polymer at 20000 ps is shown in the background. (d) CO2 MSDs as a function of x for a time-interval of 20 ps averaged around 10000 ps. The dashed line indicates the initial boundary of 6FDA−6FpDA. (e) Polymer segment MSDs as a function of x for a time-interval of 100 ps averaged around 10000 ps. (f) Same as part e for the 4000 and 8000 CO2 systems but compared with results averaged around 20000 ps.

polymers,3,100−102 among which the site-distribution (SD) model of Kirchheim and co-workers.103,104 This model describes sorption in terms of a Gaussian distribution of site energies, which are correlated to a volume distribution of the intermolecular space yielding different elastic distortion energies when occupied by solute molecules. Although our statistics are limited, our weighted energy probability distributions (Figure 8a) are indeed reminiscent of Gaussian forms and in the membrane core distributions, there are no evidence of the two-energy-level of the DMS model.103 Similar conclusions were obtained for O2 in ODPA−ODA33 and for CO2 in the fluorinated bulk models.35 A critical discussion of the DMS and the SD models has also been presented earlier within the framework of experimental and simulated CO2 sorption data in polysulfones. The SD model was found to describe the diffusive/elastic fraction of experimental sorption

and dilation data rather well, in agreement with modeling results.79 4.5. Gas Trajectories. As found for other small gas penetrants in glassy dense matrices,105−107 trajectories of individual CO2 molecules within the 6FDA−6FpDA membrane can be described as a combination of oscillations within available polymer free volumes and occasional jumping events upon opening of temporary channels, with some jumps being productive and others unproductive in terms of diffusion. In our models, no CO2 manages to cross either of the 1000, 2000, and 3000 CO2 membranes within the 10000 ps MD simulations, although a few of them do reach ∼60−70 Å along the x-direction. On the other hand, the 20000 ps time scale allows for some gas molecules in the 4000, 8000, and 12000 CO2 systems to cross the entire membranes. Their number logically increases with the swelling response of the matrix from ∼5 for the 4000 CO2, ∼20 for the 8000 CO2 to K


Macromolecules

Article

Figure 10. Left axis: time evolution of the average mass density distribution ⟨ρ(x)⟩ for CO2 at the start of the simulation for the (a) 4000 and the (b) 8000 CO2 system. Right axis: the corresponding polyimide ⟨ρ(x)⟩ at 100 ps. On that very short time scale, the swelling of the matrix is limited in both systems. The legend is the same on both graphs.

∼300 for the 12000 CO2 systems. Parts a and b of Figure 9 display the x-positions of some molecules as a function of time, while Figure 9c gives a direct visualization using the VMD 1.8.2 software. In spite of the fact that MD configurations have been accumulated every 20 ps and that some individual jump events are likely to occur on a shorter time scale,14 Figure 9 does illustrate the overall mechanism of CO2 mobility in glassy polymer matrices. While there are many possible paths for the molecules depending on the microvoids heterogeneity, the natural fluctuations and the swelling extent of the polymer matrices, there are no real signs of diffusional anisotropy in either the x, y or z directions. This is coherent with the mobility being restricted by the local relaxations of the polymer. As such, the motion mechanisms are clearly very close to those found in the bulk,14,16 with some regions occupied by the polymer from which the penetrants are excluded and some temporary “channels” through the system that form an interconnecting and percolating network. The CO2 mean-square displacements, MSD = ⟨(ri(t + t0) − ri(t0))2⟩, with t0 being a time-origin and t a time-interval were averaged over all gas molecules in slabs of width 5 Å around 10000 ps. They are shown as a function of x in Figure 9d for a time-interval of t = 20 ps. As expected, the overall mobility is restricted in the membrane cores. However, most of the MSD curves are found to overlap and compare very well with those in the bulk with average jump lengths of ∼5−6 Å.14 On the other hand, if larger time-intervals are considered (not shown), molecules diffuse faster in the concentrated systems. This is once again consistent with the bulk;14 i.e., molecules diffuse with shorter time-of-residences and there are more frequent jumps in the dilated systems, but the amplitudes of the individual jumps remain fairly similar. The only exception is the 12000 CO2 system in which the swelling is so large that jumps can occur on longer distances both on a short and on a large time scale. However, even in that very highly plasticized system, there are no transitions to the very smooth paths characteristic of liquid-like diffusion since the jump-mechanism remains clear to see (Figure 9b). This suggests that plasticization in the experimental pressure range6,13 is not likely to result in a change of mechanism for CO2 motion. Parts e and f of Figure 9 show the polymer segment MSD for a time-interval of t = 100 ps, averaged either around 10000 ps

or around 20000 ps. Interestingly, the enhancement of chain mobility in the interfacial region (with respect to that at the COM) seems to extend further into the membrane than the mere perturbation of slab mass density. The increase of the local diffusivity at a free surface has been reported in several studies of both rubbery and glassy polymers.30−33,61,67,108,109 It was attributed to the lower density prevailing in the interfacial region associated with a transmission of enhanced mobility through molecular connectivity.61 As mentioned in section 2.3, this “dynamical interfacial thickness” is usually at least twice as large as the thickness obtained from the mass density profile.30−33 This is clear to see in Figure 9, parts e and f, and it is further enhanced by the penetrant entering the matrix and the boundaries moving as plasticization occurs. However, as found for the bulk systems,14 the consequent dilations of the simulation boxes do not lead to a very significant increase in the mobility of the matrices, at least not in the MD time scale. As such, in the absence of CO2, the systems will probably not be able to recover their initial structure and thus the excess freevolume that has been created by the swelling will remain to a certain extent. This explains the well-known hysteresis effect seen upon desorption following conditioning of the samples by CO2.6,8 However, while Figure 9e shows the expected tendency for chain motions to increase at the higher CO2 concentrations and associated volume dilations, Figure 9f suggests that these local motions actually get slower with time. This is particularly relevant with respect to the work of Visser, Wessling, and coworkers4,12,78 who describe CO2-sorption-induced relaxations by the sum of a Fickian and two relaxational sorption regimes, albeit on much larger time scales. The first relaxational regime is a fast one and the second one is much slower. Indeed, it is known experimentally that the permeability in plasticized polymers increases as a function of time, but that the extent of this increase actually decreases with proceeding times.78 While our molecular models cannot be directly compared to experiments lasting several hours,4,12,78 it is interesting to note that, after adjusting to the entrance of the gas at the start, the relaxation modes of the polymer do tend to slow down with time afterward (Figure 9f) . 4.6. The Initial Adsorption Phase. The adsorption phase was characterized by the time evolution of the x-dependent symmetrized average CO2 mass density, ⟨ρ(x)⟩, in the first stages of the simulations for all the systems under study. In all L


Macromolecules

Article

Figure 11. Time evolution of the average mass density distributions ⟨ρ(x)⟩ for CO2 (left axis) and for the 6FDA−6FpDA polymer (right axis) in the (a) 4000 and the (b) 8000 CO2 systems. The legend is the same for both graphs. (c) Symmetrized averaged % trans as a function of x for the 6FpDA Car−C−Car−Car angles at t = 10000 ps. (d) Symmetrized averaged P2(cos θx) for the 6FpDA Car−C−Car angles at t = 20000 ps as a function of x. In parts c and d, the slab width is 5 Å and the dashed lines indicate roughly the interface regions at the relevant times.

cases, the slab width resolution for ⟨ρ(x)⟩ was 1 Å. Figure 10 compares the gas ⟨ρ(x)⟩ at 0, 2, 4, 6, 10, 50, and 100 ps for the 4000 and the 8000 CO2 systems. As found for O2 sorption,33 the penetrants reach the vicinity and enter the membranes within very few ps, while the evolution of the density profiles is slower after that. In the lowpressure systems, a gas adsorption peak is clear to see at the interface, where the penetrants gather within less than 100 ps (Figure 10a). This tends to saturate the interfacial region with gas, but visual examinations show that there are continuous exchanges between nonadsorbed gas molecules in the reservoir and those adsorbed at the polymer surface, in agreement with a Langmuir mode based on gas molecules which constantly condense and evaporate at the surface.110 In the higherpressure systems, the adsorption peak is smoother (Figure 10b) as the interface is dynamically saturated while many molecules are still left in the reservoirs. Indeed, in the 12000 CO2 system (not shown), it almost looks like a continuous plateau. In all cases, the rate of advancement in the x-direction region at the vicinity of the membranes is on the order of ∼1−2 Å ps−1, which is close to a gas phase behavior. 4.7. The Slower Diffusion Phase. The time evolution of the symmetrized average mass densities ⟨ρ(x)⟩ for the gas and the polymer analyzed over a longer time scale than in Figure 10 is given in Figures 11a-b for the 4000 and the 8000 CO2 systems. The concomitant changes in the plasticized polymer structures are displayed in Figures 11c-d.

The gas density profiles at longer times have the same characteristics than those at shorter times (Figures 10 and 11), i.e. the penetrant adsorption peaks are well visible at the interface for the lower pressure systems while they tend to be smoother for the higher pressure systems. However their evolution is clearly slower and the polymer ⟨ρ(x)⟩ react to the entry of the penetrant by the increase of the interface location and a decrease in density. As noted before, the local interface relaxation is seen at the early stages of the simulations (w parameter in eq 2 and Figure 4c), and indeed, the differences between 10000 and 20000 ps are much smaller than those between 1000 and 10000 ps. Experimentally, phenomena at the basis of plasticization are more pronounced in the beginning,9 which seems to be also the case in the models. In addition, the change in interface location results in shifts in the gas adsorption peaks. In order to assess the effects on the matrix structure, Figure 11c displays the symmetrized average percent of Car−C−Car− Car angles between both aromatic rings in the 6FpDA moieties in the trans conformation (defined as being trans if −60° < torsional angle 0) = C0. In the present work, the boundaries are moving, so we cannot possibly apply these analyses to large time intervals. This implies that the infinite condition on the right side; i.e., the penetrants should not exit on the other side of the membrane is de facto obeyed. However, this is also means that the statistics will be poor and that it is unlikely that the diffusion coefficient Dpenetrant(t), which is known to vary with t for short time-intervals, will be able to reach its limiting value at large t.32,33 Analyses were thus carried out on various small time-intervals ranging from 100 to 2000 ps and defining the width of the membrane as being the 2h value (Figure 4b) in the middle of the time-interval. Gas concentrations vs time curves were obtained by artificially labeling gas probes as being on the “left”, xi(t0) < −h, or on the “right”, xi(t0) > h at a time origin t0 and then following and symmetrizing the x-evolution of these t0-labeled atoms. Such curves obtained on the 8000−10000 ps time-interval with h being that at 9000 ps are displayed in Figure 12a. The lines are N


Macromolecules

Article

upon desorption following conditioning of the samples by CO2, which is linked to the systems not being able to recover their initial structure and thus some of the excess free-volume that has been created by the swelling remaining. Because of a strong chemical potential gradient, the penetrants underwent a very fast sorption at the polymer surface in the very first stages of the simulations. This saturated the interfacial region with gas, while a quasi-equilibrium was quickly established with the reservoir. The penetrant adsorption peaks were well visible at the interface for the lower pressure systems while they tended to be smoother for the higher pressure systems. Once the initial adsorption phase was completed, a second slower diffusion mode then took over. It could be described by fitting an erfc solution of the onedimensional diffusion equation in a semi-infinite system to the density of those gas molecules having entered the membrane over specific time-intervals, providing that these time-intervals were small enough because of the moving boundaries. The parameters obtained from the fits were the time-dependent concentration at the interface, which follows a Langmuir model as a function of pressure, and the time-dependent diffusion coefficient of the penetrant in the membrane. The latter varies with t for short time-intervals since the diffusion is governed by the anomalous regime and we could not reach its limiting value at large t. However, the results were very reasonable if compared to experimental and bulk values. Such membrane models can only realistically represent the very early stages of CO2 sorption, but compared to experiment, all significant features are being reproduced in our models. Other un-cross-linked polyimide matrices are currently under study in order to assess the effect of the polymer structure on the sorption and dilation. However, since CO2-induced swelling and plasticization are detrimental to membrane performance as the polymer loses its cohesion and shape-discriminating ability, much experimental work has been devoted to reducing it by, e.g., thermal annealing or cross-linking.11,113 These more complicated structures should also be envisaged at a later stage.

case, diffusion is still governed by the anomalous behavior and such membrane models can only realistically represent the initial immediate stages.

5. CONCLUSIONS A series of large-scale MD simulations of CO2 transport in a fully atomistic glassy 6FDA−6FpDA polyimide membrane model were carried out in order to assess the uptake mechanisms of the penetrant at the polymer interface. A ∼50000-atom pure 6FDA−6FpDA membrane model with explicit interfaces was first created by using a wall-compression technique directly inspired from the experimental solventcasting process. Six different CO2 concentrations were then placed in the reservoirs on each side of the membrane and the basic mechanisms underlying gas adsorption, penetration and diffusion through the membrane models were followed over a time scale of 10000 to 20000 ps. CO2-induced volume dilation was found to occur immediately in all six models under study, with a progressive shift of the interface location toward larger values associated with a widening of the interface and a decrease in density. The volume swelling ranged from about 3% to 17% in the various systems, with the exception of the highest-pressure 12000 CO2 system, for which a much larger volume increase (up to 50%) suggested a quasi-supercritical behavior for CO2. Sorption was characterized as a function of time within the framework of the moving interface locations. There was no apparent delay between sorption and dilation kinetics. The gas uptake curves C(t) for the five lower-pressure systems were typical of small-gas penetrants into glassy matrices, and the CO2 concentration range was 0−160 cm3(STP) cm−3, i.e. very close to that of the associated 6FDA−6FpDA bulk models. The C(t) were concave to the time-dependent pressure axis P(t) and fitted well to the functional form of the DMS sorption model, while the time-dependent solubility coefficient S(t) was in good agreement with experimenta data. However, there was no evidence of two different environments for the penetrants at the molecular level. Indeed, the Boltzmann weighted insertion energy probability density distributions, which are Gaussian, support the site-distribution SD model. Void-space and excess chemical potential analyses confirmed that the holes within the matrix are getting larger as sorption proceeds. However, since they are very favorable insertion sites, they are immediately occupied by penetrant molecules. The only way to accommodate more penetrant is to further swell the polymer matrix. In all cases, the 12000 CO2 system was found to behave differently from the other systems, thus confirming that it is in a regime more characteristic of a solvent-like than a gas-like behavior for the penetrant. In terms of local mobility, the CO2 trajectories displayed the oscillations in voids and occasional jumps mechanism typical for small penetrants in glassy matrices. In the swollen systems, the molecules diffused with shorter time-of-residences and there were more frequent jumps, but the amplitudes of the individual jumps remained fairly similar. Even in the highly plasticized system, there were no transitions to the paths characteristic of liquid-like diffusion. Dilations led to very limited structural relaxations of the matrices rather than to larges changes in the chain structures. As such, the increase in the mobility of the matrices was also quite small. In addition, after adjusting to the entrance of the gas at the start, the relaxation modes of the polymer tended to slow down with time. This is probably at the basis of the hysteresis effect seen

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] Telephone: +33 4 79758697. Fax: +33 4 79758614. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was granted access to the HPC resources of CCRT/ CINES/IDRIS under the allocations 2011- and 2012-095053 made by GENCI (Grand Equipement National de Calcul Intensif), France. The MUST cluster at the University of Savoie (France) is also acknowledged for the provision of computer time

■

REFERENCES

(1) Paul, D. R. Sep. Purif. Rev. 1976, 5, 33−50. (2) Koros, W. J.; Mahajan, R. J. Membr. Sci. 2000, 175, 181−196. (3) Yampolskii, Y.; Pinnau, I.; Freeman, B. D. Materials Science of Membranes; John Wiley & Sons Ltd.: Chichester, U.K., 2006. (4) Wessling, M.; Huisman, I.; Van der Boomgaard, T.; Smolders, C. A. J. Polym. Sci., Part B: Polym. Phys. 1995, 33, 1371−1384. (5) Bos, A., Ph.D. Thesis. University of Twente: Twente, The Netherlands, 1996. O


Macromolecules

Article

(6) Coleman, M. R.; Koros, W. J. Macromolecules 1997, 30, 6899− 6905. (7) Coleman, M. R., Ph.D. Thesis. University of Texas: Austin, TX, 1988. (8) Coleman, M. R.; Koros, W. J. Macromolecules 1999, 32, 3106− 3113. (9) Bos, A.; Pünt, I. G. M.; Wessling, M.; Strathmann, H. J. Membr. Sci. 1999, 155, 67−78. (10) Ismail, A. F.; Lorna, W. Sep. Purif. Technol. 2002, 27, 173−194. (11) Wind, J. D.; Sirard, S. M.; Paul, D. R.; Green, P. F.; Johnston, K. P.; Koros, W. J. Macromolecules 2003, 36, 6433−6441. (12) Visser, T.; Wessling, M. Macromolecules 2007, 40, 4992−5000. (13) Pandiyan, S.; Brown, D.; Neyertz, S.; Van der Vegt, N. F. A. Macromolecules 2010, 43, 2605−2621. (14) Neyertz, S.; Brown, D.; Pandiyan, S.; Van der Vegt, N. F. A. Macromolecules 2010, 43, 7813−7827. (15) Ghosh, M. K.; Mittal, K. L. Polyimides: fundamentals and applications; Marcel Dekker, Inc.: New York, 1996. (16) Neyertz, S.; Brown, D. Macromolecules 2010, 43, 9210−9214. (17) Van der Vegt, N. F. A.; Briels, W. J.; Wessling, M.; Strathmann, H. J. Chem. Phys. 1999, 110, 11061−11069. (18) Natarajan, U.; Mattice, W. L. J. Membr. Sci. 1998, 146, 135−142. (19) Chang, J.; Han, J.; Yang, L.; Jaffe, R. L.; Yoon, D. Y. J. Chem. Phys. 2001, 115, 2831−2840. (20) Ayyagari, C.; Bedrov, D. Polymer 2004, 45, 4549−4558. (21) Ijantkar, A. S.; Natarajan, U. Polymer 2004, 45, 1373−1381. (22) Prathab, B.; Aminabhavi, T. M.; Parthasarathi, R.; Manikandan, P.; Subramanian, V. Polymer 2006, 47, 6914−6924. (23) Doruker, P.; Mattice, W. L. Macromolecules 1998, 31, 1418− 1426. (24) Queyroy, S.; Neyertz, S.; Brown, D.; Müller-Plathe, F. Macromolecules 2004, 37, 7338−7350. (25) Harmandaris, V. A.; Adhikari, N. P.; Van der Vegt, N. F. A.; Kremer, K. Macromolecules 2006, 39, 6708−6719. (26) Hess, B.; León, S.; Van der Vegt, N. F. A.; Kremer, K. Soft Matter 2006, 2, 409−414. (27) Harmandaris, V. A.; Reith, D.; Van der Vegt, N. F. A.; Kremer, K. Macromol. Chem. Phys. 2007, 208, 2109−2120. (28) Marcon, V.; Fritz, D.; van der Vegt, N. F. A. Soft Matter 2012, 8, 5585−5594. (29) Neyertz, S.; Douanne, A.; Brown, D. Macromolecules 2005, 38, 10286−10298. (30) Neyertz, S.; Douanne, A.; Brown, D. J. Membr. Sci. 2006, 280, 517−529. (31) Neyertz, S. Soft Mater. 2007, 4, 15−83. (32) Neyertz, S.; Brown, D. Macromolecules 2008, 41, 2711−2721. (33) Neyertz, S.; Brown, D. Macromolecules 2009, 42, 8521−8533. (34) Brown, D. The gmq User Manual Version 4: available at http:// www.lmops.univ-savoie.fr/brown/gmq.html, 2008. (35) Pandiyan, S.; Brown, D.; Van der Vegt, N. F. A.; Neyertz, S. J. Polym. Sci., Part B: Polym. Phys. 2009, 47, 1166−1180. (36) Hammonds, K. D.; Ryckaert, J.-P. Comput. Phys. Commun. 1991, 62, 336−351. (37) Ciccotti, G.; Ferrario, M.; Ryckaert, J. P. Mol. Phys. 1982, 47, 1253−1264. (38) Clark, M.; Cramer, R. D., III; Van Opdenbosch, N. J. Comput. Chem. 1989, 10, 982−1012. (39) Zhang, Z.; Duan, Z. J. Chem. Phys. 2005, 122, 214507. (40) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, England, 1987. (41) Ewald, P. P. Ann. Phys. 1921, 369, 253−287. (42) Smith, W. Comput. Phys. Commun. 1992, 67, 392−406. (43) Fincham, D. Mol. Simul. 1994, 13, 1−19. (44) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (45) Brown, D.; Clarke, J. H. R. Comput. Phys. Commun. 1991, 62, 360−369. (46) Neyertz, S.; Brown, D. Macromolecules 2004, 37, 10109−10122. (47) Neyertz, S. Macromol. Theory Simul. 2007, 16, 513−524.

(48) Neyertz, S.; Brown, D. J. Chem. Phys. 2001, 115, 708−717. (49) Flory, P. J. The Statistical Mechanics of Chain Molecules; Hanser Publishers: New York, 1988. (50) Helfand, E.; Tagami, Y. J. Chem. Phys. 1972, 56, 3592−3601. (51) Helfand, E.; Tagami, Y. J. Chem. Phys. 1972, 57, 1812−1813. (52) Ten Brinke, G.; Ausserre, D.; Hadziioannou, G. J. Chem. Phys. 1988, 89, 4374−4380. (53) Bitsanis, I.; Hadziioannou, G. J. Chem. Phys. 1990, 92, 3827− 3847. (54) Vacatello, M. Macromol. Theory Simul. 2001, 10, 187−195. (55) Vacatello, M. Macromol. Theory Simul. 2002, 11, 53−57. (56) Mischler, C.; Baschnagel, J.; Binder, K. Adv. Colloid Interface Sci. 2001, 94, 197−227. (57) Izumisawa, S.; Jhon, M. S. J. Chem. Phys. 2002, 117, 3972−3977. (58) Coleman, M. R.; Koros, W. J. J. Polym. Sci., Part B: Polym. Phys. 1994, 32, 1915−1926. (59) Costello, L. M.; Koros, W. J. J. Polym. Sci., Part B: Polym. Phys. 1995, 33, 135−146. (60) Mansfield, K. F.; Theodorou, D. N. Macromolecules 1990, 23, 4430−4445. (61) Mansfield, K. F.; Theodorou, D. N. Macromolecules 1991, 24, 6283−6294. (62) Jang, J. H.; Ozisik, R.; Mattice, W. L. Macromolecules 2000, 33, 7663−7671. (63) Clancy, T. C.; Jang, J. H.; Dhinojwala, A.; Mattice, W. L. J. Phys. Chem. B 2001, 105, 11493−11497. (64) Kikuchi, H.; Kuwajima, S.; Fukuda, M. J. Chem. Phys. 2001, 115, 6258−6265. (65) Doruker, P.; Mattice, W. L. Macromol. Theory Simul. 2001, 10, 363−367. (66) Jain, T. S.; De Pablo, J. J. Macromolecules 2002, 35, 2167−2176. (67) Xu, G.; Mattice, W. L. J. Chem. Phys. 2003, 118, 5241−5247. (68) Kikuchi, H.; Fukura, M. KGK Kautsch. Gummi Kunstst. 2004, 57, 416−422. (69) Baljon, A. R. C.; Van Weert, M. H. M.; Barber DeGraaff, R.; Khare, R. Macromolecules 2005, 38, 2391−2399. (70) Ochoa, J. G. D.; Binder, K.; Paul, W. J. Phys.: Condens. Matter 2006, 18, 2777−2787. (71) Peter, S.; Meyer, H.; Baschnagel, J. J. Polym. Sci., Part B: Polym. Phys. 2006, 44, 2951−2967. (72) Morita, H.; Tanaka, K.; Kajiyama, T.; Nishi, T.; Doi, M. Macromolecules 2006, 39, 6233−6237. (73) Haile, J. M. Molecular Dynamics Simulation: Elementary Methods; Wiley: New York, 1992. (74) Brown, D.; Neyertz, S. Mol. Phys. 1995, 84, 577−595. (75) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38. (76) Kratochvil, A. M.; Damle-Mogri, S.; Koros, W. J. Macromolecules 2009, 42, 5670−5675. (77) Kratochvil, A. M.; Koros, W. J. Macromolecules 2010, 43, 4679− 4687. (78) Wessling, M.; Huisman, I.; Van der Boomgaard, T.; Smolders, C. A. J. Appl. Polym. Sci. 1995, 58, 1959−1966. (79) Hö lck, O.; Siegert, M. R.; Heuchel, M.; Bö hning, M. Macromolecules 2006, 39, 9590−9604. (80) Böhning, M.; Springer, J. Polymer 1998, 39, 5183−5195. (81) Hiltner, A.; Liu, R. Y. F.; Hu, Y. S.; Baer, E. J. Polym. Sci., Part B: Polym. Phys. 2005, 43, 1047−1063. (82) Crank, J. C. The Mathematics of Diffusion, 2nd ed.; Oxford University Press: Oxford, U.K., 1979. (83) Tsige, M.; Grest, G. S. J. Chem. Phys. 2004, 120, 2989−2995. (84) Paul, D. R. Ber. Bunsen-Ges. 1979, 83, 294−302. (85) Kanehashi, S.; Nagai, K. J. Membr. Sci. 2005, 253, 117−138. (86) Jordan, S. M.; Koros, W. J. J. Membr. Sci. 1990, 51, 233−247. (87) Jordan, S. M.; Henson, M. A.; Koros, W. J. J. Membr. Sci. 1990, 54, 103−118. (88) Bondar, V. I.; Kamiya, Y.; Yampolskii, Y. P. J. Polym. Sci., Part B: Polym. Phys. 1996, 34, 369−378. P


Macromolecules

Article

(89) Wang, R.; Cao, C.; Chung, T.-S. J. Membr. Sci. 2002, 198, 259− 271. (90) Hibshman, C.; Mager, M.; Marand, E. J. Membr. Sci. 2004, 229, 73−80. (91) Boyd, R. H.; Pant, P. V. K. Macromolecules 1991, 24, 4078− 4083. (92) Pant, P. V. K.; Boyd, R. H. Macromolecules 1993, 26, 679−686. (93) Bondi, A. J. Phys. Chem. 1964, 68, 441−451. (94) Van Krevelen, D. W. Properties of polymers: their correlation with chemical structure; their numerical estimation and prediction from additive group contributions, 3rd completely revised ed.; Elsevier: Amsterdam, 1990. (95) Widom, B. J. Chem. Phys. 1963, 39, 2808−2812. (96) Krishna, R.; Wesselingh, J. A. Chem. Eng. Sci. 1997, 52, 861− 911. (97) Thompson, A. P.; Ford, D. M.; Heffelfinger, G. S. J. Chem. Phys. 1998, 109, 6406−6414. (98) Thompson, A. P.; Heffelfinger, G. S. J. Chem. Phys. 1999, 110, 10693−10705. (99) Van der Vegt, N. F. A.; Briels, W. J.; Wessling, M.; Strathmann, H. J. Chem. Phys. 1996, 105, 8849−8857. (100) Sanchez, I. C.; Lacombe, R. H. Macromolecules 1978, 11, 1145−1156. (101) De Angelis, M. G.; Merkel, T. C.; Bondar, V. I.; Freeman, B. D.; Doghieri, F.; Sarti, G. C. J. Polym. Sci., Polym. Phys. Ed. 1999, 37, 3011−3026. (102) De Angelis, M. G.; Merkel, T. C.; Bondar, V. I.; Freeman, B. D.; Doghieri, F.; Sarti, G. C. Macromolecules 2002, 35, 1276−1288. (103) Kirchheim, R. Macromolecules 1992, 25, 6952−6960. (104) Gotthardt, P.; Grüger, A.; Brion, H. G.; Plaetschke, R.; Kirchheim, R. Macromolecules 1997, 30, 8058−8065. (105) Müller-Plathe, F. Acta Polym. 1994, 45, 259−293. (106) Gusev, A. A.; Müller-Plathe, F.; van Gunsteren, W. F.; Suter, U. W. Adv. Polym. Sci. 1994, 116, 207−247. (107) Gusev, A. A.; Suter, U. W.; Moll, D. J. Macromolecules 1995, 28, 2582−2584. (108) Jang, J. H.; Mattice, W. L. Polymer 1999, 40, 4685−4694. (109) Torres, J. A.; Nealey, P. F.; De Pablo, J. J. Phys. Rev. Lett. 2000, 85, 3221−3224. (110) Moore, W. J. Physical Chemistry, 5th ed.; Prentice-Hall Inc.: Englewood Cliffs, NJ, 1972. (111) Küntz, M.; Lavallée, P. J. Phys. D: Appl. Phys. 2004, 37, L5−L8. (112) Coleman, M. R.; Koros, W. J. J. Membr. Sci. 1990, 50, 285− 297. (113) Omole, I. C.; Miller, S. J.; Koros, W. J. Macromolecules 2008, 41, 6367−6375.

Q


Molecular Dynamics Study of Carbon Dioxide Sorption and

Recommend Documents