A molecular dynamics simulation study of surface effects on gas permeation in free-standing polyimide membranes

Journal of Membrane Science 280 (2006) 517–529

A molecular dynamics simulation study of surface effects on gas permeation in free-standing polyimide membranes Sylvie Neyertz ∗ , Anthony Douanne, David Brown Laboratoire Mat´eriaux Organiques a` Propri´et´es Sp´ecifiques (LMOPS), UMR CNRS 5041, Universit´e de Savoie, Bˆat IUT, 73376 Le Bourget du Lac Cedex, France Received 16 December 2005; received in revised form 7 February 2006; accepted 9 February 2006 Available online 6 March 2006

Abstract A 141100-atom model of a glassy ODPA–ODA polyimide free-standing membrane, corresponding to a thickness of two average radii of gyration for the 40-mers chains, has been studied using molecular dynamics (MD) simulations. Due to the large-scale of the fully atomistic model, a parallelized particle-mesh technique using an iterative solution of the Poisson equation had to be implemented for the efficient evaluation of the electrostatic interactions. With flattened-chain configurations, the density was found to adjust itself naturally in the middle of the membrane to ∼95% of the ODPA–ODA experimental value. At the free-standing surfaces, the density profile became sigmo¨ıdal, indicating surface roughness. For comparison, two isotropic bulk models, one at the “normal” density as obtained for ODPA–ODA under ambient conditions and the other one at 95% of the normal-density, were built. Small gas probes were inserted into all three models in order to investigate whether the interfacial structure of the glassy free-standing membrane can influence penetrant transport. Gas diffusion in the low-density part of the interface was found to be very fast. The limiting value for the gas diffusion coefficient Dmembrane is only attained when the probes enter more dense regions in the membrane. Indeed, Dmembrane compares well with Dbulk obtained for the 95%-density bulk system, i.e. about twice that in the normal-density bulk. Solubility is larger in the membrane than in both bulk models, thus suggesting an effect of chain flattening in addition to the density. Adsorption is particularly high at the free-standing interfaces. © 2006 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics simulations; Large-scale fully atomistic polyimide; Glassy; Free-standing surfaces; Gas permeation

1. Introduction Polyimides are high-performance macromolecules which are often used in thin dense membrane applications [1–3], such as gas-separation of oxygen or nitrogen from air or the purification of natural gas [1]. However, their permeation properties are strongly affected by the considerable number of parameters that can influence diffusion, selectivity or other physical properties over the entire processing stage [4–8]. Polyimide membranes are often prepared by the so-called “solvent casting method”, in which they are cast on a substrate from solutions prepared by dissolution of the polyimide into a solvent such as m-cresol or N-methylpyrrolidone [9]. Different casting surfaces will potentially lead to different properties. It is known for example that membranes which are cast on glass result in stiffer, stronger and

∗

Corresponding author. Tel.: +33 4 79 75 86 97; fax: +33 4 79 75 86 14. E-mail address: [email protected] (S. Neyertz).

0376-7388/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.memsci.2006.02.011

more oriented films than when they are cast on liquids, such as mercury [8]. We have recently carried out a comparative study of gas permeation in two fully atomistic models of a glassy polyimide using molecular dynamics (MD) simulations [10] in order to test for possible “skin-effects” [3,11–13] related to the presence of interfaces with specific features. The 40-mer ODPA–ODA homopolyimide (Fig. 1) was chosen as a test case since experimental data for both gas permeability and diffusion in this specific polyimide are available [14,15]. In addition, the preparation procedure for atomistic models of ODPA–ODA in the pure bulk had been earlier optimized and validated with respect to other experimental properties [16–18]. The first system studied [10] was an isotropic 27-chain 56025-atom bulk model of the ODPA–ODA amorphous phase [10,18], which was periodic in three dimensions. The second system was a 68-chain 141100atom model of an actual membrane, which was created using an original procedure based on the experimental solvent casting process [10]. The membrane had a width of 13 nm, i.e. two

518

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

Fig. 1. The (ODPA–ODA) polyimide chain.

average radii of gyration for the ODPA–ODA molecules. It was periodic in the x- and y-directions, and confined between two solid walls in the z-direction in order to be consistent with the presence of the glass surface. Despite changes in configurations and high-density interfaces in the vicinity of the walls, the permeability of small gas probes was found to be similar in both bulk and membrane systems. In the former case [10], calculations were speeded-up by a factor of 70 using a soft-repulsive Weeks–Chandler–Andersen (WCA) potential for the non-bonded interactions. While this allowed for diffusive equilibrium and good statistics to be attained in both bulk and membrane models, the polymer membrane had to remain constrained between both walls in order to compensate for the absence of attractive forces, and thus only the presence of a high-density interface was tested [10]. In the present work, we include more realistic Lennard–Jones (LJ) and Coulombic interactions for non-bonded interactions. This allows us to remove both walls around the membrane, so that the density at the interface and the interior will adjust itself naturally. Such a model is likely to be more representative of reality where the film is peeled off the glass plate and local rearrangement, particularly at the surface, can occur. Several examples of atomistic simulations of free-standing polymer membranes can be found in the literature. Here, we restrict ourselves to polymer/vacuum and polymer/small penetrants interfaces, and do not consider, for example, polymer/polymer interfaces where the key features are adhesion and miscibility [19–23]. Simulations of polymer/vacuum in the melt [24–26] or in the glassy state [27,28] have shown that, although densities decrease monotonically across the polymer–vacuum interface, the chain structures are qualitatively similar to the flattened configurations which develop near a solid wall [29], albeit of more moderate amplitude. Usually, the interfacial width is increased in a free-standing film and flexible polymer chain ends have a tendency to segregate [25,30], unless they are sufficiently attractive [31]. Mobility is also affected and the apparent glass transition, which is dependent on the thickness of the film, is usually lower than in the bulk [32–34]. Crystallization can even be initiated in the surface region and propagate into the interior of the thin film [35,36]. MD simulations have been used to characterize the surface interactions of thin films of random copolymers based on styrene, butadiene and acrylonitrile with solvents [37], as well as self-associating polymers [38]. However, in actually very few cases up to now, and to our knowledge only in rubbery matrices [39,40], have fully atomistic free-standing membrane models been used to study gas transport. In the present paper, we make a comparative study of gas permeation in bulk models of the glassy ODPA–ODA polyimide with that in the corresponding free-standing membrane model. The latter originates from the earlier 141100-atom membrane

model of the same polymer confined between two walls [10] and the relevant computational details are given in Section 2. Removing the walls results in the relaxation of the polyimide matrix and the high-density interfaces; these effects are analysed in Section 3.1. In Section 3.2, diffusion and solubility results for highly mobile gas probes in the bulk and free-standing membrane models are discussed. It should be noted that the free-standing polyimide model is over 20–50 times larger than the typical sizes used for these glassy systems [17,41], and that consequently these calculations with full excluded-volume and electrostatic interactions have proven extremely expensive in terms of computational time (∼26,000 h monoprocessor on an IBM Power 4 for 1.1 ns of simulation). Within this context, relatively slow penetrants such as oxygen [15] could not be considered, but the use of small gas probes allowed for significant diffusion to be attained through the highly rigid membrane on the limited timescale accessible to the fully atomistic MD simulations. 2. Computational details 2.1. Force-field The chemical structure of an (ODPA–ODA) chain is shown in Fig. 1 [9,17,18,42]. Each chain has a total of 40 monomers and 2075 atoms. The gas probe was modeled as a single atom. All calculations were performed using the MD code of the gmq package [43] in its parallel form, ddgmq [44], on the French IDRIS (Orsay) and CINES (Montpellier) supercomputing centres as well as on local resources at the University of Savoie. The force-field and the parameters for the polyimide are the same as described before [10,17,18], with “bonded” anglebending, torsional and out-of-plane interactions arising from near-neighbour connections in the structure. The so-called “non-bonded” excluded-volume and electrostatic potentials are applied to all atom pairs separated by more than two bonds on the same molecule or belonging to different molecules. The model gas probe is also identical to that used previously [10]. Its small ˚ allows it to be highly mobile in the membrane. size, σ = 1.88 A, Standard combination rules are used for all cross terms [18], except for gas–gas interaction parameters which are set to zero to represent an ideal gas. As explained in Section 1, “non-bonded” interactions in the gas permeation study of the confined membrane were described solely by the purely repulsive WCA potential [10], which is a reasonable approximation to the more realistic excluded-volume Lennard–Jones and electrostatic forms when weakly interacting gases and rigid matrices are being investigated [18]. The low computational cost of the WCA form allowed a steady state uptake of gas to be achieved in these MD simulations [10], but

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

the polymer membrane had to be constrained between rigid walls to maintain its density. Removing the walls and allowing for the polymer to relax in a realistic way means that we have to revert to the much slower LJ (Eq. (1)) and electrostatic (Eq. (2)) combination of non-bonded potentials in order to maintain cohesion in the membrane:
   6  σ 12 σ ULJ (rij ) = 4ε (1) − rij rij Ucoul (rij ) =

qi qj 4πε0 rij

(2)

where rij is the distance between atoms i and j, ε the well-depth and σ the distance at which the LJ potential is zero, qi and qj the partial charges on atoms i and j and ε0 is the vacuum permittivity. ULJ (rij ) and Ucoul (rij ) are the corresponding energies. Considering the size of the polymer membrane, this amounts to a single integration time-step of 10−15 s taking about 85 s monoprocessor CPU time on an IBM Power 4. The Ewald summation method, which is routinely used to evaluate the electrostatic interactions in reasonably small models [45,46], was applied to the bulk models under study. On the other hand, a parallelized particle-mesh technique using an iterative solution of the Poisson equation was implemented for the membrane since it is much more efficient than the Ewald sum for large systems and it is also naturally applicable to models which are periodic in two dimensions [47]. Using this technique, both the short-range and the long-range contribution of the electrostatic potential are calculated in real space. The former is a sum of pair interactions, while the latter is handled with a smooth projection of discrete point charges to a grid [48] followed by a multigrid-approach to resolving Poisson’s equation in order-N operations [49]. 2.2. The free-standing membrane and the bulk models The free-standing membrane originates from the procedure described previously [10]. The average radius of gyration S2 1/2 ˚ for an ODPA–ODA chain was shown to be equal to ∼58.4 A at 700 K by using a well-documented “hybrid Pivot Monte Carlo-Molecular Dynamics (PMC/MD) single-chain sampling” approach in the melt [10,16–18,50–55]. Since S2 1/2 is usually quoted as being the minimum distance for the influence of the interface [28,29,56–61], the size of the model was based on a length of ∼2S2 1/2 for the cell, i.e. a total of 141100 atoms in 68 chains. The membrane was then created using an original approach designed to mimick loosely the experimental solvent casting process [10]. The polymer chains were generated with the hybrid PMC/MD single-chain sampling method in the melt and placed at a density corresponding to a 10% (w/w) solution. Following the progressive introduction of excludedvolume (using the WCA soft-repulsive form for the non-bonded potentials), an impenetrable wall was added on either side of the polymer and the membrane was compressed in the z-direction until the density in its middle reached the experimental value of 1368 kg m−3 [14], while the basis vectors in the x- and ydirection remained constant. Internal stresses were then released

519

with a constant-volume (NVT) run of 3000 ps at 700 K. At that stage, the Lennard–Jones and Coulombic non-bonded potentials were introduced at 700 K for 250 ps. Finally, the system was cooled down to 300 K and allowed to settle for 100 ps, the drift in coulombic energy being only 0.03% over the last 50 ps. In all cases, the temperature was maintained through loosecoupling to a heat bath [62] with a coupling constant equal to 0.1 ps. The final size of the confined membrane model was ˚ × 113 A ˚ × 158 A. ˚ 112 A In the free-standing membrane, the part of the wall which was interacting directly to constrain the polymer, but was invisible to the gas probe, was removed. Only the middle four layers of the initial 10-layer tetrahedral diamond lattice arrangement [10], i.e. a total of 1680 atoms linked through flexible bonds and bends, were retained in order to prevent gas probes from escaping into the vacuum. However, the wall atoms were moved away from the ˚ effectively prepolymer at such a distance that the cut-off of 10 A vented them from interacting. The total width of the remaining ˚ At that stage, the parameters associated wall amounted to ∼14 A. with the lattice diffusion multigrid procedure for the calculation ˚ −1 for the β of electrostatic interactions [47] were set to 0.24 A ˚ ˚ ˚ for coefficient and to hx = 0.700 A, hy = 0.707 A and hz = 0.702 A the mesh spacing in the x-, y- and z-directions, respectively. Particle charges were interpolated onto the 64 closest grid points of the mesh using a Gaussian charge spreading function. The optimized value for the charge-spreading diffusion coefficient ˚ 2 per iteration and the number of steps for the was D = 0.0578 A diffusion scheme was Nt = 72. Upon removal of the wall constraint, the polymer was found ˚ to relax in the z-direction from a total length of ∼130 to ∼140 A. The density of the free-standing membrane relaxed within less than 100 ps and the final size of the membrane + wall model was ˚ × 113 A ˚ × 180 A. ˚ A schematic representation adjusted to 112 A of its left interface is shown in Fig. 2a. For comparison, the left interface of the confined membrane [10] is presented in Fig. 2b. It is clear that in the free-standing membrane, the interface has expanded significantly and has become a lot rougher than its confined counterpart. The ODPA–ODA bulk simulations are based on the 6225atom model also used previously [18], which had a density of 1377 kg m−3 , i.e. within 0.7% of the experimental value of 1368 kg m−3 [14]. The Ewald summation method was used to evaluate electrostatic interactions with the same parameters as before [18], the only difference with our former calculations being the gas probe interaction parameters. This model will be referred to afterwards as the “normal-density” bulk. In order to make comparisons with the lower density of the free-standing polymer membrane, the “normal-density” bulk ˚ was driven smoothly from a side length of ∼41.4 to ∼42.1 A, and was allowed to relax for 100 ps in the pure state. This second model will be referred to as the “95%-density” bulk with a final value of 1301 kg m−3 . 2.3. Insertion of gas probes The number of gas probes np to insert [18,63] into the bulk models was set to 25. This is sufficient to give good statistics

520

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

Fig. 2. Compared close-ups of the free-standing ODPA–ODA membrane surface (a) on the left side and that of the corresponding confined membrane surface (b) ˚ from the polyimide center-of-masses. [10]. Both snapshots have been taken at the end of their respective production runs and show the left interfacial zone from ∼50 A ˚ along the y-axis and the color code is the following: cyan, C; red, O; blue, N and white, H. These schematic representations are They have dimensions of ∼30 A displayed using the VMD 1.8.2 software [90]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

whilst implying reasonable external pressures [18]. The number of gas probes absorbed into the polymer, at a temperature T, is related to the external pressure P applied and is directly linked to the probability of insertion for the gas into the polymer, pip . For small penetrants, this probability can be obtained from Widom’s test-particle insertion method [18,64,65],    U pip = exp − (3a) kB T where U is the change in potential energy associated with a trial insertion of the gas probe into the polymer matrix. pip is directly related to the excess chemical potential of the gas in the polymer, µex by:    U µex = −RT ln exp − (3b) = −RT ln(pip ) kB T Considering the volume of the polymer Vp , the external pressure is thus: P≈

np kB T pip Vp

(4)

As explained before [18], in a simulation using 3D periodic boundary conditions, the external pressure and the gas concentration can be (and often are) set independently. In our calculations, we aim to have a consistency between these variables by using Eq. (4). In a membrane model, the external pressure and the gas concentration are linked as in reality. One thus has to also consider the number of gas probes in the gas phase ng as well as

the volume of the gas phase, Vg . Assuming a probability of insertion of 1 into the gas phase, i.e. ideality, it can be shown that: P≈

(np + ng )kB T Vg + pip Vp

(5)

Widom’s test particle insertion method was used on all the stored configurations of each production run for the pure polymer systems with 500,000 trials per configuration using the parameters associated with the specific gas probe of this work. pip  was found to be 0.14 for the “normal-density” bulk model, which corresponds to an external pressure of ∼110 × 105 Pa (∼110 bars) when 25 gas atoms are included. The 95%-density bulk was run under constant-volume conditions, so no external pressure was necessary to maintain the density. For the membrane model, pip  was obtained as a function of z for bin ˚ and revealed a much more complicated profile widths of 1 A ˚ × 113 A ˚ ×1A ˚ (see Fig. 11 later). Each slab of volume 112 A was thus multiplied by its respective pip  and the total contribution summed. From the reverse use of Eq. (5) and considering the same P of ∼110 × 105 Pa (∼110 bars), this led to a total of np + ng = 3280 gas probes to be inserted into the free-standing membrane. Since it has been shown that introducing penetrants (which do not interact with each other) on one side or on both sides of a membrane model leads to analogous permeation curves [40], 1640 probes were randomly inserted into the free space on one side and the remaining 1640 into the free space on the other side of the polymer membrane.

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

It is important to note that the traditional experimental setup uses a pressure gradient in order for the permeation rate to be measured [66]. However, this is not necessary in the case of an atomistic simulation, as the gas atoms are continually entering, leaving and diffusing through the membrane even in the absence of an actual pressure gradient; indeed this is the fundamental spontaneous process that is probed by applying the relatively small pressure gradients used in real experiments. As we record the positions of all atoms as a function of time, it is straightforward to see how the density distribution of a particular subset of the gas atoms changes with time. Consequently, no additional external pressure gradient was applied. Following equilibration, the production run for each bulk model was carried out. The “normal-density” + gas model was run for 3000 ps under NPT conditions with an applied pressure P of 110 × 105 Pa (110 bars), while the “95%-density” + gas model was run for 1000 ps under NVT conditions. The much more expensive membrane + gas model was simulated for a total period of 1100 ps. In all cases, configurations were stored at 1 or 5 ps intervals, and thermodynamic and conformational data every 1 ps for post-analysis.

3. Results and discussion The free-standing membrane model is obviously very similar in terms of chain structure and configurations to the confined membrane it originates from [10]. It is clear that the influence of the preparation procedure for such a glassy fully atomistic film is a key factor, but, to date, relatively few techniques have been described in the literature. One common approach is that of extending the z-axis in order to eliminate interactions of the parent chains with their images in one direction [22,25,26,37,38]. However, it is very unlikely that rigid chains, such as polyimides, would relax much at ambient temperatures in the timescale of a few nanoseconds. This also applies to the so-called healing method where a snapshot from an equilibrated bulk film is duplicated in the z-direction and the two films are merged, thus increasing the thickness of the film [25]. A procedure related to ours has been used for example by Kikuchi et al. [67,68] who compressed a bulk polyisoprene chain with a repulsive wall, albeit starting from a box big enough to contain a single-chain under the periodic boundary conditions. The future in the field certainly lies in multiscale approaches, which involve reverse mapping from equilibrated coarse-grained models using dynamic MC [21,25,69] or MD [70,71] simulations. However, they require well-parameterized coarse-grained models of the polymers in question, which are well beyond the scope of this work. In the present paper, we investigate if the interfacial structure of a free-standing membrane can influence gas permeation. Section 3.1 will point out the actual differences between the freestanding membrane and its confined counterpart [10]. The gas probe solubility and diffusion through the free-standing model will then be compared to those of the corresponding bulk models in Section 3.2.

521

Fig. 3. Mass density distribution for the polyimide as a function of the distance z from the center-of-mass of the membrane. The curve has been symmetrized about the membrane COM and normalized with respect to the density of the ˚ The black reference bulk system, ρexp = 1377 kg m−3 . The slab width is 1 A. circles are fits of the sigmoidal part of ρ(z) to the form of Eq. (6) (see text for details).

3.1. Characterization of the free-standing polyimide membrane 3.1.1. Densities and interfacial thickness The previously prepared [18] 6625-atom “normal-density” bulk model has a density of ρbulk = 1377 kg m−3 . The confined membrane was prepared with a similar density at the center of the membrane but had a higher-density (up to 1840 kg m−3 ) layer of ˚ thickness in the vicinity of each planar surface [10]. Upon ∼6 A removal of the walls, the polymer matrix loses these high-density surface layers and the middle part of the model relaxes also towards a lower average value of 1301 kg m−3 . This decrease is a direct result of the flattening of the chains during the membrane construction phase [10], and results in densities approximately 5% lower than the experimental bulk value of 1368 kg m−3 [14]. For this reason, the “95%-density” bulk model with a density of 1301 kg m−3 was also prepared and examined. Fig. 3 shows the mass density distribution, ρ(z), of the polymer chains in the membrane model, relative to the bulk density ρbulk , as a function of the distance, z, from the middle of a given ˚ to the center-of-mass (COM) of the membrane. slab of width 1 A The flattening of the chain configurations leads to more apparent fluctuations, of ∼30 kg m−3 , as a function of z than for the con˚ fined membrane, but the ρ(z) remain fairly linear up to ∼60 A from the COM. No smooth oscillations of the density profile are evident and, consequently, there is no ordered chain layering in the core of the membrane. As seen before [22,25–28,30,31,33,34,38,40,67,69,72], the density profile at the free-standing surface is sigmo¨ıdal. Fig. 3 exemplifies the difficulties associated with precisely determining the boundaries of the membrane, and there are several possible definitions for the interfacial thickness. For example, it has been defined as being the distance over which mass density falls from its bulk value to zero [27,28]. Despite its very large magnitude, our membrane does not attain true bulk behaviour at the centre due to the method of preparation, the large S2 1/2 of the

522

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

ODPA–ODA molecules, and the rigidity of the chains. If we use instead the average relaxed density of 1301 kg m−3 for the middle of the membrane as our reference point, then ρ(z) smoothly ˚ thus leading to an interfacial thickdecrease from 60.5 to 72 A, ˚ ness of 12.5 A. As reported elsewhere [27], it is significantly smaller than the overall chain dimensions. Unlike rubbery polymers [10,25–27,33,34,69], no tendency for certain species, such as chain ends or special groups, to migrate to the surface was identified in the case of our glassy model. Another way to identify the interfacial thickness [30,31,38,72] is to fit the sigmo¨ıdal part of the mass density curve ρ(z) to the following hyperbolic equation [73,74]: ρ(z) = ρmiddle

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

(6)

where ρmiddle is the density of the middle region of the film, h the location of the interface and w is the interfacial width. ρ(z) does fit very well to the form of Eq. (6) and the resulting curve is displayed as filled circles in Fig. 3. The optimized parameters are ˚ for the location of the interface ρmiddle = 1306 kg m−3 , h = 66.9 A ˚ for the interfacial width. The value of ρmiddle is in and w = 4.5 A good agreement with the actual average of 1301 kg m−3 .A third way to locate the interface [30] is to determine the point where dρ(z)/dz is at its minimum value near the free surface. Here, the ˚ In all cases, the width of the minimum is situated at h = 65.5 A. free-standing interface is indeed larger than that of the confined interface [10,30]. 3.1.2. Structure To characterize chain alignment, the second order Legendre polynomial functions, P2 (cos θ α ), were calculated where the angle θ α for a triplet of atoms {i, j, k} is that between the αaxis and the vector between atoms i and k [75]. This analysis was carried out for all so-called “pivot angles” which includes (see Fig. 1) the flexible C–O–C ether bridges of the ODA and ODPA moieties and the rigid C–N–C ODA–ODPA linkages [17]. The P2 function is especially interesting since its limiting values, 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 θ α . Fig. 4 shows that the free-standing membrane retains the features of the confined system, i.e. a tendency towards a parallel alignment of the chains with respect to the surface in the z-direction originating from the flattening induced in the chains during the compression step of the model preparation. However, in the free-standing membrane, there is not the same precipitous fall in P2 (cos θ α ) towards −1/2 at the interface as found in the confined model (see Fig. 4 of our previous paper [10]) where the presence of the wall forces a near perfect parallel alignment of the chains to it. While a pronounced alignment of the chains in the vicinity of the interface has been extensively described in the literature for confined systems [10,27,29,56,60,61], it is also known to be present but attenuated in a free-standing membrane [22,25–27,31,33], as shown in Fig. 4. Indeed, the tendency for parallel orientation has been shown to decrease with increased roughness [72], which can clearly be related to Fig. 2. No specific features such as end-beads which orient perpendicular to

˚ They Fig. 4. Average P2 (cos θ z ) for pivot angles situated in slabs of width 2 A. are displayed as a function of the z distance from the position of the pivot-angle middle atom to the membrane COM. Ether bridges and ODPA–ODA links are COC and CNC angles, respectively. The curves have been symmetrized about the membrane COM. The standard errors are smaller than the size of the symbols.

the surface [25] or become almost random [26] in flexible polymers were identified in the glassy polyimide. Unlike P2 (cos θ z ), the P2 (cos θ x ) and P2 (cos θ y ) functions were superimposable, thus showing that there was no preference with respect to either x or to y in the membrane model. As far as polymer chain configurations are concerned, the average mean-square radius of gyration in the free-standing ˚ 2 , is practically the same as that found membrane, S2  = 2750 A ˚ 2 [10], indicating that in the confined membrane, S2  = 2760 A the relaxation that takes place is not due to changes in the conformations of these stiff chains. It should be noted that, in both cases, the component of the average mean-square radius of gyra2 , represents only ∼1% tion perpendicular to the interface, S⊥ of the total of S2 , which confirms the flattened nature of the polyimide chains. 3.1.3. Dynamics In their investigation on free surfaces of glassy atactic polypropylene [28], Mansfield and Theodorou had reported enhancements in the small mean-square displacements (MSD) of atoms in the vicinity of the interface. While it is out of the question to follow long-range atomic polyimide diffusion given the timescale of our MD simulations, there is indeed a consistent trend in the local mobility depending on the distance of the atoms to the interface. This is displayed in Fig. 5, where the MSDs, (ri (t + t0 ) − ri (t0 ))2 , are recorded for t = 10, 50 and 100 ps for all ˚ averaged polymer atoms situated, at t0 , in slabs of width 10 A over all possible time origins. For comparison, the polyimide density profile is also given in Fig. 5. While the absolute values for the polyimide MSDs are very small, there is undoubtedly an enhancement of atomic mobility close to the interface, which Mansfield and Theodorou attributed to the lower density prevailing there [28]. However, similar to their and other work [75], this dynamical feature extends further into the polymer than the mere perturbation of density, and what has been called the “dynamical interfacial thickness” is indeed about twice as large as the thickness obtained from the mass density profile. Such a

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

523

in polyimides [18], albeit with shorter times of residence in the holes [10]. As in the confined membrane, the z-component of the gas probe MSDs for penetrants situated in the middle of the membrane is found to be lower by ∼20% than the x- and y-components. This was attributed to the anisotropic configurations of the membrane and to a slight channeling effect along the x–y plane [10]. The x-, y- and z-components of the corresponding MSDs in the bulk systems are lower for the normal-density bulk whereas they are slightly higher in the 95%-density bulk with respect to the membrane. However, both bulk systems do display an isotropic diffusion since the three components are the same for a given bulk model. Fig. 5. Left axis: polyimide mean-square displacements (MSD) as a function of the z distance to the membrane COM for all atoms situated in a slab of width ˚ at time t0 . Three time-intervals are shown, i.e. t = 10, 50 and 100 ps. The 10 A MSDs have been averaged over all polyimide atoms and all possible time origins of the production run. Standard errors are smaller than the size of the symbols. Right axis: polyimide mass density ρ(z) as a function of z with a slab width of ˚ All curves have been symmetrized about the membrane COM. 1 A.

transmission of enhanced mobility has been attributed to molecular connectivity [28] and the increase of local diffusivity at the surface has been reported in several other studies of both rubbery and glassy polymers [30,32,34]. The concept of a fluidlike interfacial region has even been used to explain free-standing film behaviour at the vicinity of the glass transition temperature [33], but in the present work, the temperature is much too low to see such phenomena. 3.2. Permeation of small gas probes 3.2.1. Trajectories Typical trajectories for individual gas probes traversing the free-standing membrane are displayed in Fig. 6. They are quite characteristic of small mobile penetrants in glassy matrices [16,18] with the gas probes oscillating within available free volumes and jumping between different voids when the limited fluctuations of the glassy matrix allow for the temporary opening of channels [65,76,77]. Although the gas probe dimensions are smaller, Fig. 6 is qualitatively similar to the motion of helium

Fig. 6. Typical trajectories of gas probes along the z-direction of the freestanding membrane model.

3.2.2. Diffusion coefficients in the bulk Gas diffusion coefficients in the bulk models were obtained by averaging the gas MSDs over all penetrants and over all possible time origins of the production runs, t0 . The self-diffusion coefficient D is obtained from Einstein’s equation [18], and is a good approximation for the true diffusivity [65] as long as there are no interactions between the gas probes [78,79]. In the present case, the crossover from the anomalous to the Fickian diffusion regime occurs in less than 100 ps, as shown by the fact that the slope of log (MSD) versus log (t) goes to one [18]. The resulting diffusion coefficients are Dbulk = 1.0 × 10−4 cm2 s−1 for the normal-density bulk and Dbulk = 2.2 × 10−4 cm2 s−1 for the 95%-density bulk. These values were further confirmed by using the probability density distribution of displacement vector components, which can be fitted by a single Gaussian curve in the Fickian regime [18,80], as well as by calculating the self part of the van Hove correlation function, Gs (r, t) [18,81,82]. As expected, both approaches give the same Dbulk than those evaluated using the MSD versus time curves. 3.2.3. Diffusion coefficient in the free-standing membrane Considering the set-up of the membrane model, the gas concentration, c, is expected to vary as a function of time and of the z-position; it will thus be referred to as c(z, t). However, the simplest and most common way of solving the one-dimensional Fick’s law diffusion equation is through the identification of a regime in which the diffusion coefficient D is assumed to be independent of both z and concentration [79,83]. In terms of time t, this is indicated by a linear relationship between the weight gain of the polymer and t1/2 . The number uptake of gas probes in the membrane over the whole simulation run is presented as a function of t1/2 in Fig. 7. As before, such an analysis requires the definition of the width of the membrane L [10]. Since several possible definitions for the membrane thickness have been put forward (Section 3.1.1) and the surface of the membrane is highly heterogeneous (Fig. 2), we ˚ corresponding to the start of the dense considered L/2 = 60.5 A ˚ obtained from Eq. (6) and part of the membrane, L/2 = 66.9 A ˚ referring to the point where dρ(z)/dz passes through L/2 = 65.5 A a minimum at the interface. In addition, the number uptake was ˚ which studied under the hypothetical condition that L/2 = 50 A, would be situated well inside the dense part of the model. The uptake curves differ at short times but the limiting slopes are very similar, thus showing that the definition of the location of

524

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

Fig. 7. The number uptake of gas probes by the membrane as a function of the square root of time, t. Due to the uncertainty in the definition of the membrane thickness, L, several values related to the interface have been considered (see text ˚ which is situated and legend for details), in addition to an hypothetical L/2 = 50 A well into the dense part of the model. The actual uptake curves are displayed with symbols and lines whereas the straight lines are the corresponding linear fit to t1/2 .

the interface does not fundamentally affect the analyses. When ˚ the curves do not L/2 is defined as being 60.5, 65.5 or 66.9 A, 1/2 scale linearly as t at short times (t < 225 ps) with the critical differences being at very short times (t < 20 ps). This is related to the initial filling-up of the interface by the gas where the diffusion coefficient varies with z. As the definition of L/2 gets closer to the dense part of the interface, this non-t1/2 dependence of diffusion disappears, as can be seen from the uptake curve ˚ This indicates that the diffusion of the hypothetical L/2 = 50 A. coefficient settles to a constant value in the dense part of the ˚ uptake membrane. The progressive tendency of the L/2 > 60.5 A 1/2 curves to become linear to t means that the transition from a large value for the diffusion coefficient in the gas phase to a low constant value inside the membrane occurs over a certain range. It is clear that, despite the choice of a fast-diffusing model gas probe, the available simulation time did not allow for the steady state regime to be reached. Using results from a numerical procedure (see later), this can be estimated at 3.5 ns. However, the equilibrium and non-equilibrium uptake curves have been shown to be superimposable in the case of the confined membrane [10], thus suggesting that the non-equilibrium experiment is sufficient for our purpose. If D is independent of the concentration and of the z-position within the membrane, and diffusion in the gas phase is considerably higher, the diffusion coefficient can be evaluated from fits of the time-dependent density distributions to the following solution of the one-dimensional diffusion equation in a semi-infinite system [79,83]:    z c(z , t) = c0 erfc √ (7) 4Dt In Eq. (7), z is the coordinate in a reference system where is the left edge of a medium which extends to z = +∞. In the case of a finite-length membrane of length L, this solution can only be used up until gas probes start to exit from the far

z = 0

Fig. 8. Mass density distributions of labeled gas probes as a function of the z-position in the membrane at time-intervals of 50, 100, 200 and 300 ps. The actual profiles (displayed with symbols) have been averaged over all possible time origins by labeling gas probes which are inside or outside of the membrane, ˚ at each time origin. They are shown as a function of defined by L/2 = 66.9 A, ˚ and the left and right profiles have been z = z + (L/2). The slab width is 1 A symmetrized. The lines are fits to the form of Eq. (7) with a constant diffusion coefficient of D = 1.9 × 10−4 cm2 s−1 .

side of the membrane. With the parameters for the gas probe used in this work and the actual thickness of the membrane, this restricts the analyses to time-intervals shorter than 500 ps. In addition, the following boundary conditions are required: c(z ≥ 0, t = 0) = 0; c(z = 0, t > 0) = c0 . Gas concentration versus time curves were obtained by artificially labeling gas probes as being on the “left”, zi (t0 ) < −L/2, or on the “right”, zi (t0 ) > L/2 at any particular time origin t0 and then subsequently following the evolution of the distribution of these t0 -labeled atoms. In order to improve the resolution, all possible time origins were used and the distributions for “left” and “right” molecules were subsequently symmetrized. They are displayed in Fig. 8, where the ˚ has been used, largest of the possible definitions of L/2 of 66.9 A along with fits to the form of Eq. (7) for each given time-interval under study. While all the fits were obtained with a constant diffusion coefficient of D = Dmembrane = 1.9 × 10−4 cm2 s−1 , they ˚ onwards, which corresponds are only valid from about z = 6.5 A ˚ As noted before, this is actually the region where to z = ±60.5 A. the density approaches that of the middle region of the film. The erfc fits are thus not valid in the smoothly decreasing ρ(z) part of the interfacial region. Fig. 9 characterizes diffusion in this lowdensity part of the polyimide by selecting those gas probes which ˚ are initially in the pure gas phase, i.e. situated at |zi (t0 )| > 72 A at time origin t0 . Their concentration versus time profile is averaged over all time origins separated by 1 ps on a total production time of 500 ps. Within 1–2 ps, the interfacial structure is clearly almost entirely filled up with gas, which means that, despite an intermediate behaviour between diffusion in the pure gas phase and that through the dense membrane, it is a lot closer to the former case. Indeed, the polyimide has a very irregular structure in the immediate vicinity of the interface (see Fig. 2) with sufficiently large and connected holes for the gas probes to move easily. This is progressively reduced as the polymer density increases but still remains a lot faster than transport in both

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

525

QL (0), was determined from the integral of the real average gas ˚ and z > 60.5 A. ˚ This distribution symmetrized over z < −60.5 A region includes both the true gas reservoir and the low-density part of the interface. At each δt, the amount of gas either remaining in the original left region, QL (t), having diffused into the membrane, QM (t), or having exited through the opposite right ˚ QR (t), was evaluated. Since diffusion interface (zmax = 60.5 A), in the interfacial part of the membrane is a lot faster than that in the membrane, the concentrations at the edges of the membrane could be safely reset at each integration loop using: c(zmin , t) = c(zmin , 0)

QL (t) QL (0)

c(zmax , t) = c(zmin , 0) Fig. 9. The evolution of the average gas density distribution as a function of the z-position in the low-density interfacial part of the membrane at time-intervals of 1, 2 and 5 ps. All curves have been averaged over as many time origins as possible and over both sides of the membrane. The 0 ps reference labels those ˚ or at z > 72 A ˚ at a given time origin. gas probes which are situated at z < −72 A

˚ 2 ps−1 ). This confirms that D(z) in bulk models (Dbulk = 1–2 A the low-density part of the interface is close to that of the ideal gas and only drops to a limiting lower value in the more dense interior regions of the membrane. As before [10], we can see in Fig. 8 that the c0 coefficient in Eq. (7) slightly decreases with increasing t. This has been linked to the fact that the stationary source boundary condition, c(z = 0, t > 0) = c0 , is not fully consistent with the conditions used here, i.e. the partial pressure of the necessarily limited number of labeled gas probes in the (fixed volume) reservoir gradually diminishes as gas diffuses into the membrane. However, this was found to be of no consequence on the value of D for the confined membrane [10]. In order to check that it is also the case here and to avoid the problem of the inappropriate boundary conditions, a second approach [10] in which the solution of the one-dimensional diffusion equation is obtained numerically [83] was attempted. The range of the membrane in the z dimension ˚ and the range in t was divided into equal intervals of δz = 1 A into intervals of δt = 0.1 ps. With finite difference methods, it is possible to obtain the concentration at the middle of a slab z and at time t + δt, c(z, t + δt) with good precision from the following numerical integration algorithm: δt (c(z + δz, t) δz2 + c(z − δz, t) − 2c(z, t))

QR (t) QL (0)

(9a) (9b)

Numerical solutions to the form of Eq. (8) obtained using the diffusion coefficient evaluated from the erfc fits, i.e. D = Dmembrane = 1.9 × 10−4 cm2 s−1 , are displayed in Fig. 10 as a function of the z-position and at different time-intervals, t. For convenience, the concentrations have been converted to zdependent mass densities so as to be compared with the actual average mass density distributions. They confirm that the value given by Eq. (7) is a rather good approximation (within statistical uncertainties) and that the non-constant c0 of Fig. 8 does not affect the diffusion coefficient. In order to further check the parameters used, the numerical integration procedure was carried out up to 1.1 ns and the predicted uptake of gas probes was found to be superimposable to the actual uptake ˚ shown in Fig. 7. Furthermore, as mentioned with L/2 = 60.5 A before, carrying on the numerical integration procedure up to 10 ns with the same parameters revealed that the steady state will probably not be reached before 3.5 ns under the conditions of the present model. A third way to estimate the diffusion coefficient in the membrane would have been to perform a time-lag analysis [66,84] by following the quantity of penetrants that have crossed the entire membrane as a function of time. Although this has been done for the confined membrane [10], the basic assumption behind this

c(z, t + δt) = c(z, t) + D

(8)

provided that (D δt/δz2 ) < 1/2 [83]. As has been found with the erfc analysis, the left “edge” of the membrane has to be situated ˚ from the COM in order to have a description at zmin = −60.5 A of the concentration versus time profiles consistent with a constant D(z) in the membrane. It can also be seen from Fig. 8 that a good choice for the initial gas concentrations at the edges of the membrane, c(zmin , 0), is ∼10 kg m−3 , which corresponds to the inflexion point for all the concentration curves when the probes enter the membrane. For the numerical solution, the initial amount of gas situated on the left side of the membrane,

Fig. 10. Same as Fig. 8 except that the lines are now from numerical solutions using the algorithm given in Eq. (8) with D = 1.9 × 10−4 cm2 s−1 and the mass density distribution is displayed as a function of z.

526

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

approach is that statistics are accumulated in the steady state of the permeation process. Consequently, we cannot apply it to the current set of non-equilibrium data. However, this former analysis only confirmed the results of both the analytical erfc and numerical approaches in the confined membrane [10], thus suggesting that the value of Dmembrane = 1.9 × 10−4 cm2 s−1 for the gas probe diffusion coefficient in the free-standing polyimide model is reliable. Dmembrane is equal to twice Dbulk for the normal-density bulk system, which shows that the presence of a free-standing interface and flattened configurations in the vicinity of a glassy polyimide film will not act as a retardant for diffusion. This is primarily governed by lower density since Dmembrane is of a similar order of magnitude than Dbulk = 2.2 × 10−4 cm2 s−1 obtained for the 95%-density bulk model. 3.2.4. Solubility in the free-standing membrane For small penetrants, such as the gas probe under study, solubility in a polyimide membrane is known to follow Henry’s law in that it is linearly proportional to the pressure, S = Sc P, where Sc is the solubility coefficient expressed in appropriate units [18,65]. It is, of course, directly linked to the average probability of insertion for the gas into the polymer, pip , which has been given in Eq. (3a). Using Widom’s test insertion technique [64,65] and 1 million attempted trials per configuration, pip  was found for the gas probe used here to be equal to 0.16 ± 0.01 in the normal-density bulk + gas and 0.20 ± 0.01 in the 95%-density bulk + gas systems. Similar analyses on the pure bulk systems gave 0.14 ± 0.01 and 0.19 ± 0.01, respectively. Although the size of the probe is not realistic, the bulk + gas pip  would correspond roughly to 0.14 × 10−5 cm3 (STP) cm−3 Pa−1 (0.14 cm3 (STP) cm−3 bar−1 ) for the normaldensity model and 0.18 × 10−5 cm3 (STP) cm−3 Pa−1 (0.18 cm3 (STP) cm−3 bar−1 ) for the 95%-density model. The related excess chemical potentials µex of the gas in the polyimide (Eq. (3b)) are thus of the order of 4.6 kJ mol−1 for the normal-density model and 4.0 kJ mol−1 for the 95%-density model. As for the confined membrane [10], Widom’s insertion technique was also applied to the free-standing membrane model as a function of ˚ and the z-position of the inserted probe with a resolution of 1 A a total of 10 million insertions. The results presented in Fig. 11 display the symmetrized pip (z) for the membrane along with the pip  for both bulk models. Superimposing the polyimide mass density curve shows that the pip (z), albeit higher than both bulk values, are fairly constant from the membrane COM up to the start of the lowdensity interface. The lower density and the flattening of the chains associated with the compression step in the membrane preparation procedure lead to a general increase of solubility, which is twice that of the normal-density polyimide bulk and 1.5 time higher than the 95%-density bulk. The associated excess chemical potential µex decreases to ∼3.1 kJ mol−1 . However, as the polyimide density decreases, the solubility then increases in a sigmo¨ıdal way, reaching even values where pip (z) > 1. The latter refers to a further decrease of the excess chemical potential and the presence of very favourable interactions between the penetrants and the polymer chains at the interface. Such inter-

Fig. 11. The average probability of insertion for the gas probes into the freestanding membrane as a function of z, pip (z), obtained using Widom insertion’s technique [64,65] (squares with a maximum standard error of 0.02). pip  in both the normal-density (medium dash) and the 95%-density bulk systems (short dash) are shown for comparison, as well as the average slab mass density for the polyimide (line). All curves have been symmetrized about the membrane COM ˚ and the slab width is 1 A.

actions get progressively less easy as the gas progresses into the denser parts of the matrix, which illustrates the fact that higher density leads to less space available for weakly interacting penetrants, i.e. lower solubility [85]. Thus in addition to faster diffusion, the low-density interface is also associated with higher solubility, which will even increase permeation. It should be noted that results presented here are for a specific model gas probe which has a higher solubility into the polyimide than, for example, helium [14]. However, this increased diffusion + solubility combination allowed for significant penetration into the glassy rigid membrane model to occur in about ∼1 ns. 4. Conclusions The purpose of this study was to allow a very large-scale initially confined glassy polyimide fully atomistic membrane model [10] to adjust itself to its natural density by removing the walls which held it in place, and test whether this relaxation had an effect on gas permeation properties. To stabilize the free-standing surfaces, realistic Lennard–Jones and Coulombic interactions had to be considered in place of the purely repulsive non-bonded potential used in the confined membrane [10]. This required the implementation of a parallelized particle-mesh technique using an iterative real-space solution of the Poisson equation since it is a lot more efficient for large systems than the commonly used Ewald sum [45,46] and can be applied to periodicity in two dimensions [47]. Despite this improvement, the timescale available to the MD simulation remained necessarily limited due to the considerable size of the membrane model. As expected, the free-standing membrane retained in a more diffuse manner the main features of its confined counterpart which had been created using a strategy designed to mimick the experimental solvent casting process [10]. Chains were still flattened and stacked parallel to the interfaces but the density in the middle of the membrane adjusted itself to ∼95% of the ODPA–ODA experimental value. At the free-standing surfaces,

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

the density profile was found to fall smoothly to zero, and this was associated with surface roughness, such as that displayed in Fig. 2a, and a relative increase in atomic mobility. Within this context, the definition of interfacial thickness was not obvious and could be estimated in several ways, but this had no effect on the permeation properties. Results for the diffusion and solubility of the gas probe in the membrane model show clearly that the low-density polyimide interface is highly permeable with both very high diffusion and high solubility. Transport is only slowed down as the gas probes enter the dense parts of the ODPA–ODA glassy membrane. In that case, permeation is mostly governed by the lower density prevailing due to the flattened chain configurations. However, although Dmembrane ≈ Dbulk for the 95%-density bulk model, the anisotropy of the membrane leads to a slight increase in solubility, thus suggesting once again a faster permeability coefficient in such structures. There are still some specific features in this polyimide membrane which need to be addressed at a later stage. Although it was designed with respect to the experimental solvent casting process, the preparation procedure of such a fully atomistic large-scale glassy model remains non-trivial. In addition, solvent evaporation was not modeled explicitly and this could influence the chain configurations. Another open question is the size of the gas probe. As explained in the text, it was kept very small in an attempt to follow gas motion through the membrane with sufficient statistics. Of course, it can be argued that larger and more realistic probes, such as O2 or N2 , would maybe modify to some small extent the structure of the rigid matrix [86–89]. With our current computational resources, this is presently a difficult problem to tackle if full electrostatic and excluded-volume are to be considered in a system of this size. A last point for discussion is that, with hindsight, considering the way in which the density and the chain alignment (Figs. 3 and 4) rise up to their limiting values within a relatively narrow range, it is tempting to think ˚ of the model film could have been disthat the central ∼60 A carded. This would certainly represent a saving of about a factor of 2 in terms of CPU time but it is not entirely obvious whether the results would have been the same. Such a question could only be answered by a study of a smaller system. In addition, with a thinner system alone, there would have been lingering doubts as to whether the true plateau had been reached. It is also worth pointing out that the actual membrane model length of ˚ remains small with respect to the typical experimental ∼140 A dimensions of a few microns for those thin-films [14]. Nevertheless, this work does answer the questions related to the effect of confinement and to the use of a purely repulsive excluded-volume potential in our initially confined ODPA–ODA polyimide membrane [10]. Acknowledgments The IDRIS (Orsay, France) and CINES (Montpellier, France) supercomputing centres and the University of Savoie are acknowledged for computer time. The Rhˆone-Alpes region is thanked for funds dedicated to computer hardware and a doctoral grant for AD.

527

References [1] H. Ohya, V.V. Kudryavtsev, S.I. Semenova, Polyimide Membranes— Applications, Fabrication and Properties, Kodansha Ltd. and Gordon and Breach Science Publishers S.A., Tokyo and Amsterdam, 1996. [2] D.R. Paul, Y.P. Yampolskii (Eds.), Polymeric Gas Separation Membranes, CRC Press, Boca Raton, FL, USA, 1994. [3] Polyimides: Fundamentals and Applications, Marcel Dekker Inc., New York, 1996. [4] G.C. Eastmond, P.C.B. Page, J. Paprotny, R.E. Richards, R. Shaunak, Molecular weight dependence of gas permeability and selectivity in copolyimides, Polymer 34 (1993) 667. [5] C. Joly, D. Le Cerf, C. Chappey, D. Langevin, G. Muller, Residual solvent effect on the permeation properties of fluorinated polyimide films, Sep. Purif. Technol. 16 (1999) 47. [6] H. Kawakami, M. Mikawa, S. Nagaoka, Gas transport properties in thermally cured aromatic polyimide membranes, J. Membr. Sci. 118 (1996) 223. [7] K.C. O’Brien, W.J. Koros, G.R. Husk, Influence of casting and curing conditions on gas sorption and transport in polyimide films, Polym. Eng. Sci. 27 (1987) 211. [8] P. Pandey, R.S. Chauhan, Membranes for gas separation, Prog. Polym. Sci. 26 (2001) 853. [9] E. Pinel, M.-F. Barthe, J. De Baerdemaeker, R. Mercier, S. Neyertz, N.D. Alb´erola, C. Bas, Copolyimides containing alicyclic and fluorinated groups: characterization of the film microstructure, J. Polym. Sci.: Part B: Polym. Phys. 41 (2003) 2998. [10] S. Neyertz, A. Douanne, D. Brown, Effect of interfacial structure on permeation properties of glassy polymers, Macromolecules 38 (2005) 10286. [11] A.M. Shishatskii, Y.P. Yampolskii, K.-V. Peinemann, Effects of film thickness on density and gas permeation parameters of glassy polymers, J. Membr. Sci. 112 (1996) 275. [12] H. Kawakami, M. Mikawa, S. Nagaoka, Gas transport properties of asymmetric polyimide membrane with an ultrathin surface skin layer, Macromolecules 31 (1998) 6636. [13] C.D. Dimitrakopoulos, S.P. Kowalczyk, Scanning force microscopy of polyimide surfaces, Thin Solid Films 295 (1997) 162. [14] E. Pinel, Relations entre architecture a` diff´erentes e´ chelles et propri´et´es physiques de membranes copolyimides. Application a` la s´eparation des gaz, Ph.D. Thesis, University of Savoie, France, 2001. [15] C. Bas, R. Mercier, J. Sanchez-Marcano, S. Neyertz, N.D. Alb´erola, E. Pinel, Copolyimides containing alicyclic and fluorinated groups: solubility and gas separation properties, J. Polym. Sci.: Part B: Polym. Phys. 43 (2005) 2413. [16] S. Neyertz, D. Brown, Preparation of bulk melt chain configurations of polycyclic polymers, J. Chem. Phys. 115 (2001) 708. [17] E. Pinel, D. Brown, C. Bas, R. Mercier, N.D. Alb´erola, S. Neyertz, Chemical influence of the dianhydride and the diamine structure on a series of copolyimides studied by molecular dynamics simulations, Macromolecules 35 (2002) 10198. [18] S. Neyertz, D. Brown, Influence of system size in molecular dynamics simulations of gas permeation in glassy polymers, Macromolecules 37 (2004) 10109. [19] B.L. Sch¨urmann, U. Niebergall, N. Severin, C. Burger, W. Stocker, J.P. Rabe, Polyethylene (PEHD)/polypropylene (iPP) blends: mechanical properties, structure and morphology, Polymer 39 (1998) 5283. [20] O. Okada, K. Oka, S. Kuwajima, S. Toyoda, K. Tanabe, Molecular simulation of an amorphous poly(methyl methacrylate)–poly(tetrafluoroethylene) interface, Comput. Theor. Polym. Sci. 10 (2000) 371. [21] J.H. Jang, W.L. Mattice, Time scales for three processes in the interdiffusion across interfaces, Polymer 40 (1999) 1911. [22] T.C. Clancy, W.L. Mattice, Computer simulation of polyolefin interfaces, Comput. Theor. Polym. Sci. 9 (1999) 261. [23] M. Deng, V.B.C. Tan, T.E. Tay, Atomistic modeling: interfacial diffusion and adhesion of polycarbonate and silanes, Polymer 45 (2004) 6399.

528

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529

[24] W.G. Madden, Monte Carlo studies of the melt–vacuum interface of a lattice polymer, J. Chem. Phys. 87 (1987) 1405. [25] P. Doruker, W.L. Mattice, Simulation of polyethylene thin films on a high coordination lattice, Macromolecules 31 (1998) 1418. [26] J. Chang, J. Han, L. Yang, R.L. Jaffe, D.Y. Yoon, Structure and properties of polymethylene melt surfaces from molecular dynamics simulations, J. Chem. Phys. 115 (2001) 2831. [27] K.F. Mansfield, D.N. Theodorou, Atomistic simulation of a glassy polymer surface, Macromolecules 23 (1990) 4430. [28] K.F. Mansfield, D.N. Theodorou, Molecular dynamics simulation of a glassy polymer surface, Macromolecules 24 (1991) 6283. [29] I. Bitsanis, G. Hadziioannou, Molecular dynamics simulations of the structure and dynamics of confined polymer melts, J. Chem. Phys. 92 (1990) 3827. [30] J.H. Jang, W.L. Mattice, The effect of solid wall interaction on an amorphous polyethylene thin film, using a Monte Carlo simulation on a high coordination lattice, Polymer 40 (1999) 4685. [31] J.H. Jang, R. Ozisik, W.L. Mattice, A Monte Carlo simulation of the effects of chain end modification on freely standing thin films of amorphous polyethylene melts, Macromolecules 33 (2000) 7663. [32] J.A. Torres, P.F. Nealey, J.J. de Pablo, Molecular simulation of ultrathin polymeric films near the glass transition, Phys. Rev. Lett. 85 (2000) 3221. [33] T.S. Jain, J.J. de Pablo, Monte Carlo simulation of free-standing polymer films near the glass transition temperature, Macromolecules 35 (2002) 2167. [34] G. Xu, W.L. Mattice, Monte Carlo simulation on the glass transition of free-standing atactic polypropylene thin films on a high coordination lattice, J. Chem. Phys. 118 (2003) 5241. [35] G. Xu, W.L. Mattice, Monte Carlo simulation of the crystallization and annealing of a freestanding thin film of n-tetracontane, J. Chem. Phys. 116 (2002) 2277. [36] G. Xu, H. Lin, W.L. Mattice, Configuration selection in the simulations of the crystallization of short polyethylene chains in a free-standing thin film, J. Chem. Phys. 119 (2003) 6736. [37] U. Natarajan, W.L. Mattice, Interaction of toluene, hexadecane, and water with the surfaces of random copolymers of styrene and butadiene, J. Membr. Sci. 146 (1998) 135. [38] C. Ayyagari, D. Bedrov, A molecular dynamics simulation study of the influence of free surfaces on the morphology of self-associating polymers, Polymer 45 (2004) 4549. [39] M. Ito, M. Matsumoto, M. Doi, Molecular dynamics simulation of polymer film, Fluid Phase Equilib. 144 (1998) 395. [40] H. Kikuchi, M. Fukura, Evaluation of transport properties of gases in rubbery polymer membrane by molecular dynamics simulations, KGK Kautschuk Gummi Kunststoffe 57 (2004) 416. [41] M. Heuchel, D. Hofmann, P. Pullumbi, Molecular modeling of smallmolecule permeation in polyimides and its correlation to free-volume distributions, Macromolecules 37 (2004) 201. [42] E. Pinel, C. Bas, S. Neyertz, N.D. Alb´erola, R. Petiaud, R. Mercier, Copolyimides with trifluoromethyl or methoxy substituents. NMR characterization, Polymer 43 (2002) 1983. [43] D. Brown, The gmq User Manual Version 3, 1999. Available at http://www.lmops.univ-savoie.fr/brown/gmq.html. [44] D. Brown, H. Minoux, B. Maigret, A domain decomposition parallel processing algorithm for molecular dynamics simulations of systems of arbitrary connectivity, Comput. Phys. Commun. 103 (1997) 170. [45] P.P. Ewald, Die berechnung optischer und elektrostatishe gitterpotentiale, Ann. Phys. 64 (1921) 253. [46] W. Smith, A replicated data molecular dynamics strategy for the parallel Ewald sum, Comput. Phys. Commun. 67 (1992) 392. [47] C. Sagui, T. Darden, Multigrid methods for classical molecular dynamics simulations of biomolecules, J. Chem. Phys. 114 (2001) 6578. [48] J.V.L. Beckers, C.P. Lowe, S.W. de Leeuw, An iterative PPPM method for simulating coulombic systems on distributed memory parallel computers, Mol. Simul. 20 (1998) 369. [49] A. Brandt, Multi-level adaptive solutions to boundary value problems, Math. Comput. 31 (1977) 333.

[50] D. Brown, J.H.R. Clarke, M. Okuda, T. Yamazaki, A molecular dynamics study of chain configurations in n-alkane liquids, J. Chem. Phys. 100 (1994) 1684. [51] D. Brown, J.H.R. Clarke, M. Okuda, T. Yamazaki, The preparation of polymer melt samples for computer simulation studies, J. Chem. Phys. 100 (1994) 6011. [52] D. Brown, J.H.R. Clarke, M. Okuda, T. Yamazaki, A large scale molecular dynamics study of chain configurations in the n = 100 alkane melt, J. Chem. Phys. 104 (1996) 2078. [53] S. Neyertz, D. Brown, A computer simulation study of the chain configurations in poly(ethylene oxide)-homologue melts, J. Chem. Phys. 102 (1995) 9725. [54] S. Neyertz, D. Brown, Erratum: a computer simulation study of the chain configurations in poly(ethylene oxide)-homologue melts, J. Chem. Phys. 104 (1996) 10063. [55] S. Neyertz, D. Brown, J.H.R. Clarke, The local energy approximation and the predictability of chain configurations in polymer melts, J. Chem. Phys. 105 (1996) 2076. [56] G. Ten Brinke, D. Ausserre, G. Hadziioannou, Interactions between plates in a polymer melt, J. Chem. Phys. 89 (1988) 4374. [57] M. Vacatello, Monte Carlo simulations of the interface between polymer melts and solids. Effects of chain stiffness, Macromol. Theory Simul. 10 (2001) 187. [58] S.K. Kumar, M. Vacatello, D.Y. Yoon, Off-lattice Monte Carlo simulations of polymer melts confined between two plates. 2. Effects of chain length and plate separation, Macromolecules 23 (1990) 2189. [59] M. Vacatello, Ordered arrangement of semiflexible polymers at the interface with solids, Macromol. Theory Simul. 11 (2002) 53. [60] C. Mischler, J. Baschnagel, K. Binder, Polymer films in the normalliquid and supercooled state: a review of recent Monte Carlo simulation results, Adv. Colloid Interface Sci. 94 (2001) 197. [61] K.F. Mansfield, D.N. Theodorou, Atomistic simulation of a glassy polymer/graphite interface, Macromolecules 24 (1991) 4295. [62] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684. [63] S. Neyertz, D. Brown, A. Douanne, C. Bas, N.D. Alb´erola, The molecular structure and dynamics of short oligomers of PMDA–ODA and BCDA–ODA polyimides in the absence and presence of water, J. Phys. Chem. B 106 (2002) 4617. [64] B. Widom, Some topics in theory of fluids, J. Chem. Phys. 39 (1963) 2808. [65] F. M¨uller-Plathe, Permeation of polymers. A computational approach, Acta Polym. 45 (1994) 259. [66] B. Flaconn`eche, J. Martin, M.H. Klopffer, Transport properties of gases in polymers: experimental methods, Oil Gas Sci. Technol. - Rev. IFP 56 (2001) 245. [67] H. Kikuchi, S. Kuwajima, M. Fukuda, Novel method to estimate solubility of small molecules in cis-polyisoprene by molecular dynamics simulations, J. Chem. Phys. 115 (2001) 6258. [68] H. Kikuchi, S. Kuwajima, M. Fukuda, Permeability of gases in rubbery polymer membrane: application of a pseudo-nonequilibrium molecular dynamics simulation, Chem. Phys. Lett. 358 (2002) 466. [69] T.C. Clancy, J.H. Jang, A. Dhinojwala, W.L. Mattice, Orientation of phenyl rings and methylene bisectors at the free surface of atactic polystyrene, J. Phys. Chem. B 105 (2001) 11493. [70] S. Queyroy, Simulations mol´eculaires dynamiques de surfaces de polym`ere amorphe: cas de la cellulose, Ph.D. Thesis, University of Savoie, France, 2004. [71] S. Queyroy, S. Neyertz, D. Brown, F. M¨uller-Plathe, Preparing relaxed systems of amorphous polymers by multiscale simulation: application to cellulose, Macromolecules 37 (2004) 7338. [72] P. Doruker, W.L. Mattice, Effect of surface roughness on structure and dynamics in thin films, Macromol. Theory Simul. 10 (2001) 363. [73] E. Helfand, Y. Tagami, Theory of the interface between immiscible polymers, J. Chem. Phys. 56 (1972) 3592. [74] E. Helfand, Y. Tagami, Theory of the interface between immiscible polymers. II, J. Chem. Phys. 57 (1972) 1812.

S. Neyertz et al. / Journal of Membrane Science 280 (2006) 517–529 [75] D. Barbier, D. Brown, A.-C. Grillet, S. Neyertz, Interface between end-functionalized PEO oligomers and a silica nanoparticle studied by molecular dynamics simulations, Macromolecules 37 (2004) 4695. [76] A.A. Gusev, F. M¨uller-Plathe, W.F. van Gunsteren, U.W. Suter, Dynamics of small molecules in bulk polymers, Adv. Polym. Sci. 116 (1994) 207. [77] A.A. Gusev, U.W. Suter, D.J. Moll, Relationship between helium transport and molecular motions in a glassy polycarbonate, Macromolecules 28 (1995) 2582. [78] H.L. Tepper, W.J. Briels, Comments on the use of the Einstein equation for transport diffusion: application to argon in AlPO4 -5, J. Chem. Phys. 116 (2002) 9464. [79] M. Tsige, G.S. Grest, Molecular dynamics simulation of solvent-polymer interdiffusion: fickian diffusion, J. Chem. Phys. 120 (2004) 2989. [80] S. Neyertz, A. Douanne, C. Bas, N.D. Alb´erola, Comparative MD simulations of PMDA–ODA and BCDA–ODA polyimides, in: M.J.M. Abadie, B. Sillion (Eds.), Proceedings of the 5th European Technical Symposium on Polyimides and High Performance Functional Polymers (STEPI 5), Montpellier, France, May 3–5, 1999, 2001, p. 133. [81] T.R. Cuthbert, N.J. Wagner, M.E. Paulaitis, G. Murgia, B. D’Aguanno, Molecular dynamics simulation of penetrant diffusion in amorphous polypropylene: diffusion mechanisms and simulation size effects, Macromolecules 32 (1999) 5017. [82] J.H.D. Boshoff, R.F. Lobo, N.J. Wagner, Influence of polymer motion, topology and simulation size on penetrant diffusion in amorphous, glassy

[83] [84]

[85]

[86]

[87]

[88] [89]

[90]

529

polymers: diffusion of helium in polypropylene, Macromolecules 34 (2001) 6107. J.C. Crank, The Mathematics of Diffusion, second ed., Oxford University Press, Oxford, 1979. M.H. Klopffer, B. Flaconn`eche, Transport properties of gases in polymers: bibliographic review, Oil Gas Sci. Technol. - Rev. IFP 56 (2001) 223. N.F.A. van der Vegt, W.J. Briels, M. Wessling, H. Strathmann, Free energy calculations of small molecules in dense amorphous polymers. Effect of the initial guess configuration in molecular dynamics studies, J. Chem. Phys. 105 (1996) 8849. H. Yamamoto, Y. Mi, S.A. Stern, A.K. St-Clair, Structure/permeability relationships of polyimide membranes. II, J. Polym. Sci.: Part B: Polym. Phys. 28 (1990) 2291. L. Jia, J. Xu, A simple method for prediction of gas permeability of polymers from their molecular structure, Polym. J. 23 (1991) 417. S.A. Stern, Polymers for gas separations: the next decade, J. Membr. Sci. 94 (1994) 1. Y. Li, M. Ding, J. Xu, Relationship between structure and gas permeation properties of polyimides prepared from oxydiphtalic dianhydride, Macromol. Chem. Phys. 198 (1997) 2769. W. Humphrey, A. Dalke, K. Schulten, VMD - Visual Molecular Dynamics, J. Mol. Graphics 14 (1996) 33.