J. Phys. Chem. B 1999, 103, 3539-3544
3539
Solubility of Small Molecules and Their Mixtures in Polyethylene Shyamal K. Nath and Juan J. de Pablo* Department of Chemical Engineering, UniVersity of Wisconsin-Madison, Madison, Wisconsin 53706 ReceiVed: October 29, 1998; In Final Form: February 26, 1999
The solubility of small molecules and their mixtures in polyethylene is studied using a recently proposed united atom force field.1 By comparing the simulation results with experimental data on solubility of small molecules in semicrystalline polyethylene, it is found that crystallinity can severely effect the solubility of a solute, especially if it is highly soluble in the polymer of interest. It is also observed that small molecules form aggregates within the polymer in the liquid phase.
1. Introduction
TABLE 1: Intramolecular Potential Energy Functions for Alkanes
In the manufacture of polyolefins, knowledge of the levels of gaseous reactants and inerts (monomer, comonomer, etc.) sorbed by polymer particles at reactor conditions is essential for developing a sound understanding of the kinetics of the polymerization process and, hence, the quality of the product resin. Solubility data are required for proper design of polymerization reactors, as well as devices for safe handling of reaction byproducts. Experimental studies of sorption of multicomponent gas mixtures in polyolefins under reactor conditions are expensive and time consuming. Further, safety considerations can add significantly to the cost of such work. For most practical industrial applications, empirical correlations and semitheoretical equations of state are often used to obtain solubility data of solutes and their mixtures in the relevant polymeric systems. In a previous study,2 we have shown that several widely used equations of state fail to provide high-accuracy predictions of phase equilibria for binary polyolefin mixtures (particularly mixtures of highly assymetric components such as ethyleneC40). In contrast, when a suitable molecular force field is available, molecular simulations have been shown to constitute a valuable predicting tool for such calculations. With the advent of novel Monte Carlo techniques such as the Gibbs ensemble,3 configurational bias methods,4 and expanded ensemble methods,5 it is now possible to use molecular simulations to generate phase diagrams for mixtures of moderately long molecules with relative ease. Interestingly, even though the extension of the Gibbs ensemble technique to simulate ternary mixtures is straightforward, only a few simulation studies exist on phase equilibria for ternary systems.6 Further, these previous studies were limited to mixtures of spherical particles. Panagiotopoulos et al.,6 for example, simulated phase diagrams for ternary Lennard-Jones spheres of various size and interaction strength assymetries. To the best of our knowledge, phase equilibria simulations of ternary chain molecular systems have not been conducted before. In this study, we present results for the phase behavior of binary and ternary mixtures of small molecules and polyethylene. Knowledge of ternary phase equilibria for alkane systems is important in both the polymer and the petrochemical industries. Considering the chemical simplicity that characterizes such systems, molecular simulations could be a useful tool in * Author to whom correspondence should be addressed.
Bond Stretching Potential
V(r)/kB ) Kr ) 96500 K/Å2
Kr (r - beq)2 2
beq ) 1.54 Å (alkanes) beq ) 1.34 Å (ethylene)
Bond Bending Potential
V(θ)/kB ) Kθ ) 62500 K/rad2
Kθ (θ - θeq)2 2
θeq ) 114.0°
Torsional Potential V(φ)/kB ) V0 + V1(1 + cos φ) + V2/(1 - cos 2φ) + V3(1 + cos 3φ) V0 ) 0 V1 ) 355.04 K V2 ) -68.19 K V3 ) 701.32 K
understanding and determining the phase behavior of these systems. Furthermore, simulations could also provide molecular level insights into the complicated phase behavior of these systems. This work, in which we study the solubility of gaseous mixtures in a polymer for several systems of engineering importance, is a first step towards that end. The paper is organized as follows. We begin with a brief description of the model and simulation methods employed to generate phase diagrams. We then present solubility results for small alkanes and their mixtures in polyethylene. Our simulation results are compared to available experimental data. We analyze solute-solute and polymer-solute intermolecular correlation functions, which are shown to exhibit some unexpected features. We conclude with a few remarks concerning the solubility of small molecules in semicrystalline polymers and the possibilities of molecular simulations in regards to the study of realistic mixtures of industrial interest. 2. Simulation Details A united-atom representation of the alkanes is adopted throughout this work. For simple n-alkanes and ethylene, we use the recently proposed NERD force field.1 The NERD force field has been shown to provide good agreement with experimental phase equilibria data of pure alkanes and their binary mixtures. A complete list of available intra- and intermolecular NERD force field parameters is provided in Tables 1 and 2. A
10.1021/jp9842796 CCC: $18.00 © 1999 American Chemical Society Published on Web 04/15/1999
3540 J. Phys. Chem. B, Vol. 103, No. 18, 1999
Nath and de Pablo
TABLE 2: Intermolecular Potential Energy Functions for Alkanes Nonbonded Interaction Potential
[(σr ) - (σr ) ]
V(r) ) 4
12
6
All Molecules, when CHn Is Not a Terminal Group σCH2 ) 3.93 Å CH2 ) 45.8 K σCH2 ) 3.79 Å
Ethylene
CH2 ) 84.7 K
Alkanes Longer than Propane CH3 ) 104.0 K σCH3 ) 3.91 Å
Lennard-Jones potential energy function is adopted to describe site-site interactions both for sites located more than three bonds apart (on the same molecule) and sites located on different molecules. For alkanes, a torsional potential energy function is imposed on rotations about carbon-carbon bonds.7 Both bond stretching and bond angle bending are manipulated by means of a harmonic potential. For nonbonded, unlike pair interactions, we use Lorentz-Berthelot combining rules, which in previous work have been shown to be suitable for alkane mixtures.2 In this work, nitrogen is modeled as a dumbbell; the proposed force field parameters for nitrogen are,8 /kB ) 36.0 K, σ ) 3.31 A, with an equilibrium bond length of 1.0897 Å and a bond spring constant of Kr ) 96500 K/Å2. As with the NERD force field for alkanes, parameters for nitrogen were arrived at by simultaneous regression of second virial coefficient and phase coexistence data. In all calculations, a cutoff radius of 10.0 Å was employed for Lennard-Jones interactions and standard tail corrections were implemented.9 All phase equilibria simulations were conducted in the Gibbs ensemble at constant pressure and temperature.3 The transfer of long alkane molecules between coexisting phases was made possible by the use of an expanded Gibbs ensemble method.5 In the expanded Gibbs ensemble method, the required transfer moves of long-chain molecules between coexisting phases are attemped by introducing intermediate states of a "tagged" chain molecule. Thus, the difficulty of inserting a whole molecule into the system at once (whose acceptance probability is prohibitively low for large molecules at intermediate to high densities) is reduced to that of transfering a few segments of the tagged chain from one coexisting phase to the other in a stepwise manner. All of the simulation techniques employed here are well established and have been published elsewhere;3-5 readers are therefore referred to the literature for additional details. For a polymer-small molecule mixture at a low to moderate pressure, the vapor phase is polymer free. To increase efficiency, we reject trial attempts to transfer the polymer molecules from the liquid phase to the vapor phase during those simulations. To verify that the restriction on the presence of the polymer in the vapor phase does not affect the simulation results appreciably, we have also conducted a few simulations where this constraint was not imposed. In order to facilitate comparison to available experimental data, simulations of solubility must be carried out at a specified gas-phase composition. Simulations of that nature, however, are difficult without prior knowledge of the relative composition of the solvents in the liquid phase. To obtain a rough estimate of the polymer-phase composition, short preliminary simulations were performed before starting an actual production run. In all our simulations, molecules were displaced by means of a hybrid Monte Carlo (HMC) procedure.10,11 Within the HMC procedure, five molecular dynamics steps were used to generate
Figure 1. Pressure vs composition diagram for a mixture of nitrogenbutane at T ) 66.3 °C. The filled symbols correspond to our simulations and the open symbols correpond to experimental data.12 The lines are provided as a guide to the eye for the simulation data.
a global trial Monte Carlo move. As part of the HMC procedure, a time step of 2.345 ps was used for the molecular dynamics trial moves. In this study, we model polyethylene as linear alkane chains of 70 carbon units. Note that, beyond this molecular weight (C70), the density of long alkane oligomers does not change significantly and we do not expect the solubility of small molecules to exhibit appreciable differences in systems of longer chains. In conducting all the simulations, we restricted our system to 20 polyethylene chains (of 70 segments each). The number of solvent particles was adjusted so that the simulation box size of the gas phase was always slightly larger than that of the liquid phase. The average number of small molecules present in the liquid phase in a binary methane-polyethylene mixture at 25 °C and 600 psig is in the range 12-18, and that in a binary ethylene-polyethylene mixture is in the range 105115. For ternary mixtures of methane, ethylene, and polyethylene, the number of solvent molecules in the polymer-rich phase ranges from 45 to 55, at 25 °C and 600 psig. Equilibrium averages were collected for about 2 × 107 simulation steps, of which 5% were volume moves, 5-20% were HMC moves, and the rest were transfer moves. 3. Results and Discussion 3.1. Solubility. Before presenting results for the solubility of small molecules in polyethylene, we show several simulated binary phase diagrams for which direct comparison to experimental data can be made. Figure 1 shows a pressure vs composition diagram for a binary nitrogen-butane mixture at T ) 66.3 °C. The filled symbols correspond to results of simulations, and the open symbols correspond to experimental data.12 Despite the fact that simulations were conducted without the use of any binary parameters, agreement with experiment is satisfactory. At high pressures, simulations tend to slightly overpredict both the vapor and liquid concentrations of nitrogen. However, the uncertainty of simulations is also higher near the critical point, which precludes us from making more definite statements regarding the accuracy of the force field used. Figure 2 shows simulated and experimental13 phase diagrams for a binary ethylene-tetracontane system at P ) 96.5 bar. Agreement between simulated and experimental data is remarkable; simulations provide quantitative results and are even able to describe the LCST phenomena observed for this mixture.
Solubility of Small Molecules in Polyethylene
Figure 2. Temperature vs composition diagram for a mixture of ethylene-tetracontane at P ) 96.5 bar. The meaning of the symbols and lines is the same as in Figure 1. Experimental results are from de Loos et al.13
Figure 3. Solubility of small molecules in polyethylene at 25 °C as a function of pressure. The meaning of symbols and lines is the same as in Figure 1. The squares correspond to results for methane and the circles to results for nitrogen. Experimental results are from Li and Long.14
Our results for the nitrogen-butane and ethylene-tertacontane mixtures lend credence to the hypothesis that the force field employed in this work is capable of describing linear alkanes and their mixtures with nitrogen and ethylene with relatively high accuracy. More importantly, given the quality of binary predictions, we do not foresee a need to introduce additional parameters to describe ternary mixtures of these components. Experimental solubility data are available for methane and nitrogen in polyethylene at room temperature (25 °C) and at pressures up to 1200 psig. These data are shown in Figure 3, along with the results of our simulations. The open symbols correspond to experiment, and the filled symbols correspond to the results of Monte Carlo simulations. Note that the experimental data were collected from semicrystalline polyethylene samples and corrected for crystallinity according to
100% amorphous polymer ) solubility in semicrystalline polymer (1) 1-R where, R is the fraction of crystalline polymer (55%, for the data in Figure 3). This correction for crystallinity is only valid
J. Phys. Chem. B, Vol. 103, No. 18, 1999 3541
Figure 4. Solubility of binary methane-nitrogen mixtures in polyethylene at 25 °C as a function of pressure. The meaning of the lines and symbols is the same as in Figure 1. Experimental results are from Li and Long.14 The circles are results for pure nitrogen, and the squares are results for an equimolar mixture of nitrogen and methane.
if the crystalline part of the polymer acts as a strictly inert zone which does not influence solubility in any way. This is only true if the crystallinity of the polymer is low; unfortunately, the samples on which Li and Long conducted their experiments were highly (55%) crystalline. In contrast, our simulations were carried out directly on amorphous samples. The simulated solubility of nitrogen at low pressures is slightly below the experimental data of Li and Long,14 whereas it is slightly above at high pressures. For methane, we find that simulated results are always slightly higher than the experimental data. If the uncertainity of experiments and simulations (as shown in the figure) is taken into account, however, the agreement between both can be considered to be satisfactory. Figure 4 shows the solubility of a binary methane-nitrogen mixture (50/50 mole fraction) in polyethylene at various pressures and at 25 °C. Simulated solubilities are slightly below experiment at all pressures. One of the most interesting aspects of the results shown in Figure 4 is the synergism of nitrogen and methane; the solubility of the methane-nitrogen binary mixture in polyethylene is higher than the solubilities of either pure methane or pure nitrogen. This synergism is also apparent in the experimental data, and it is captured by our simulations. One way in which crystalline domains could influence solubility would be by becoming partially amorphous in response to the presence of a solute. In that event, corrections to solubility data in the form of eq 1 would predict a solubility of the gas above that for a completely amorphous polymer sample. Such a mechanism could explain the slight discrepancy we observe in Figure 4 for solubility of the binary methanenitrogen mixture. The crystalline domains of polyethylene can also be viewed as a network, which would give rise to elastic forces that would oppose swelling. A schematic representation of the swelling ability of a semicrystalline network and a pure amorphous polymer is shown in Figure 5. In this case, solubility data (on the basis of amorphous material) would be lower than those corresponding to a completely amorphous or less crystalline material. Furthermore, this effect would become more important as the swelling of the polymer increases, i.e. for more soluble fluids. Figure 6 shows the solubility of ethylene in polyethylene from both simulation and experiment. The severe effects of constrained swelling can be observed on this mixture; simulated
3542 J. Phys. Chem. B, Vol. 103, No. 18, 1999
Figure 5. Schematic diagram illustrating the difference between restricted and unrestricted swelling of a semicrystalline polymer. The small solvent molecules do not dissolve in the crystalline regions of a semicrystalline polymer; such regions remain unperturbed by the presence of a swelling agent. Since the crystalline domains are rigid, and since they are connected to the amorphous regions, their presence places constraints on the swelling of the amorphous parts of the polymer. The presence of restrictions can reduce significantly the swelling ability of a semicrystalline polymer.
Nath and de Pablo
Figure 7. Solubility of binary methane-ethylene mixtures in polyethylene at 25 °C as a function of pressure. The meaning of the lines and symbols is the same as in Figure 1. Experimental results are from Li and Long.14
Figure 6. Solubility of ethylene in polyethylene at 25 °C as a function of pressure. The meaning of the lines and symbols is the same as in Figure 1. Experimental results are from Li and Long.14
results are well above the corresponding experimental data. Since molecules such as methane and nitrogen have a much smaller solubility (in polyethylene) than ethylene, the presence of crystalline domains does not effect their solubility as much; this would explain why constrained swelling effects are not observed in Figures 3 and 4. The idea of constrained swelling is not new; Rogers et al.15 were the first to consider the crystalline regions of a polymer as crosslinks that restrict the swelling of a network. They found that the swelling of semicrystalline polyethylene in a poor solvent (methylbromide) could be easily explained using the Flory-Huggins theory, with the assumption that swelling only occured in the amorphous region. However, the same theory failed to explain the swelling of a semicrystalline sample by a good solvent (hexane); the authors then postulated that this was due to the existence of a network. Since the work of Rogers et al., this view has been adopted by many authors. Most recently, Doong and Ho16 have shown that, by considering an elastic factor arising from the network formed by crystalline regions, regular solution theories can explain the solubility of highly soluble penetrants in semicrystalline polymers. Figure 7 compares simulated solubilities for a binary mixture of methane-ethylene (50/50 by mole fraction) in polyethylene to experimental data. As expected, experimental data are consistently below the simulation results. However, the effects
Figure 8. Site-site intermolecular correlation functions for a binary mixture of methane and polyethylene. The solid line corresponds to the methane-methane correlation function, the dashed line corresponds to the methane-polyethylene correlation function, and the small-dashed line correspond to the methane-polyethylene end group correlation function.
of constrained swelling for the methane-ethylene system are much smaller than they are for pure ethylene and agreement with experimental binary data is reasonable. 3.2. Pair Correlation Functions. Some interesting features are revealed by the solute-solute and solute-polymer intermolecular pair correlation functions of the methane, ethylene, and polyethylene mixtures. Figure 8 shows intermolecular pair correlation functions for the methane-polyethylene binary mixture at 25 °C and 600 psig. All the calculations for correlation functions are performed at the equilibrium mixture compositions reported in the previous section. The solid line corresponds to the methane-methane correlation, the dashed line corresponds to the methane-polyethylene correlation, and the small dashed line corresponds to the methane-polyethylene
Solubility of Small Molecules in Polyethylene
Figure 9. Site-site intermolecular correlation functions for a binary mixture of ethylene and polyethylene. The meaning of the lines is the same as in Figure 8.
end-group correlation. The methane-polyethylene function represents correlations between any methane particle and any site of a polyethylene molecule, whereas the methanepolyethylene end-group function only includes correlations between any methane particle and an end site of any polyethylene molecule. The first peak of the methane-methane pair correlation function is large and has a much higher value than that of the methane-polyethylene function. A stronger first peak for the methane-methane interaction indicates that methane particles prefer to stay close to each other and form “aggregates” when dissolved in polyethylene. Another interesting observation from Figure 8 pertains to the relative strength of the first peak of the methane-polyethylene and the methane-polyethylene end-group correlations. The first peak of the methane-polyethylene end-group function is higher than that of the methane-polyethylene, indicating that methane aggregates prefer to reside near the end sites of the polymer chains. Of course, it must be noted that, since we model polyethylene molecules as alkane chains of 70 units, we have a higher concentration of end groups than would be found in a truly polymeric system. However, the concentration of end groups in the melt should not change the observation that methane particles prefer the end groups of chain molecules. As shown in Figure 8, the first peak of the methanepolyethylene function is much wider than that of the methanemethane or the methane-polyethylene end-group correlation. The wider first peak is due to the fact that interaction sites of polyethylene chains are of the fused-sphere type, the bond length connecting two methylene sites being smaller than the radius of a methylene site. Figure 9 shows ethylene-ethylene, ethylene-polyethylene, and ethylene-polyethylene end-group site-site intermolecular pair correlation functions for a binary ethylene-polyethylene mixture at 25 °C and 600 psig. The general features of these functions are similar to those for the methane-polyethylene system. The first peak of the ethylene-ethylene correlation function is the highest, that of the ethylene-polyethylene correlation function is the weakest, and that of the ethylene-
J. Phys. Chem. B, Vol. 103, No. 18, 1999 3543 polyethylene end-group correlation function is in between. However, the relative differences between the peaks for the ethylene-polyethylene system are much smaller than those for the methane-polyethylene mixture. It should be pointed out that ethylene has a higher solubility in polyethylene than methane, which partly explains why ethylene has a lower tendency to form aggregates in polyethylene. As should be expected, the first peaks in Figure 9 are wider than those in Figure 8, due to the dumbbell shape of an ethylene molecule. Figure 10 shows pair correlation functions for the ternary methane-ethylene-polyethylene mixture at 25 °C and 600 psig. From Figure 10a, we observe that the first peaks of both the methane-methane and the ethylene-ethylene correlation functions are more pronounced than those of the methanepolyethylene or ethylene-polyethylene correlation functions. The first peak of the methane-ethylene correlation function falls between that of the methane-methane and the ethyleneethylene correlation functions. From Figure 10b we observe that, as for the binary mixtures, the first peak of the solutepolyethylene correlation function is weaker than that of the solute-polyethylene end-group correlation function. From the relative strengths of the first peaks of the correlation functions, we can argue that, as for the binary mixtures, the solutes tend to form aggregates near the end sites of the polymer chains in the ternary methane-ethylene-polyethylene system. A schematic representation of this tendency is shown in Figure 11. The solute particles reside around the end group of a polyethylene molecule. 4. Conclusions Monte Carlo simulations with an expanded Gibbs ensemble method have been conducted to study the solubility of small molecules and their mixtures in amorphous polyethylene. A recently proposed united-atom force field was used for this study. Results of simulations have been compared to available experimental data; the systems studied here included various binary and ternary mixtures of polyethylene with methane, ethylene, and nitrogen. In general, good agreement is obtained between simulation and experiment for binary mixtures of small molecules and long alkanes. These results serve to validate the force field used here. For sparingly soluble gases, agreement between available experimental data and simulation results of solubility of gases in long-chain polyethylene is satisfactory. However, for highly soluble gases, simulated solubilities are higher than the corresponding experimental values. We attribute these discrepancies to constrained swelling effects. An interesting synergism has been observed for mixtures of nonpolar gases such as methane and nitrogen. In agreement with experimental data, the simulated solubility of such mixtures is greater than that for those of the individual components. Interesting observations have also been made regarding the effects of crystallinity on the solubility of gaseous solutes. Available experimental data for solubility in polyethylene correspond to highly crystalline samples. It is found that for sparingly soluble gases (such as methane, nitrogen, or their mixtures), the effects of crystallinity are not very significant. To obtain the solubility of such solutes in a 100% amorphous polyethylene, a simple correction based on the assumption that crystalline domains behave as completely inert zones is sufficient. For highly soluble gases, however, crystallinity does affect solubility by constraining the ability of the polymer to swell. These effects will be addressed in a future publication
3544 J. Phys. Chem. B, Vol. 103, No. 18, 1999
Nath and de Pablo
Figure 10. Site-site intermolecular correlation functions for a ternary mixture of methane, ethylene, and polyethylene. (a) The dashed-dotted line corresponds to methane-methane correlations, the dashed line corresponds to methane-ethylene correlations, the solid line corresponds to ethyleneethylene correlations, the dotted line corresponds to methane-polyethylene correlations, and the small-dashed line corresponds to ethylenepolyethylene correlations. (b) The dotted line corresponds to the methane-polyethylene end group correlation function, the dashed line corresponds to the ethylene-polyethylene end group correlation function, the small-dashed line corresponds to the methane-polyethylene correlation function, and the solid line corresponds to the ethylene-polyethylene correlation function.
References and Notes
Figure 11. A schematic diagram depicting the aggregation of solute molecules in a ternary methane, ethylene, and polyethylene mixture.
by examining the solubility of other systems and by generating the necessary experimental data in our laboratories. Several solute-solute and solute-polymer site-site pair correlation functions have been determined. Our results indicate that solute molecules tend to form aggregates of their own, and these exhibit a tendency to reside in the vicinity of the end groups of long alkane molecules. Acknowledgment. The authors are grateful to Union Carbide for partial support of this work. This work was also supported by a PECASE award from the NSF to Juan de Pablo.
(1) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1998, 108, 9905. (2) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J.; Patramai, I. Ind. Eng. Chem. Res. 1998, 37, 3195. Note that, accidentally, the bond length for ethylene was not mentioned in Table 1 of this paper. The equilibrium bond length for ethylene is 1.34 Å. (3) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (4) Laso, M.; de Pablo, J. J.; Suter, U. W. J. Chem. Phys. 1992, 97, 2817. (5) Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1996, 105, 4391. (6) Tsang, P. C.; White, O. N.; Perigard, B. Y.; Vega, L. F.; Panagiotopoulos, A. Z. Fluid Phase Equilibr. 1995, 107, 31. (7) Jorgensen, W. I.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638. (8) Nath, S. K.; Rivera, J. L.; Alejandre, J.; de Pablo, J. J. Mol. Phys., submitted. (9) Note that, in this work, we use a cutoff radius of 10 Å instead of the longer cutoff radius (13.5 Å) used in our earlier works. We have found that the results of simulation are independent of the cutoff radius used. Use of a smaller cutoff radius helps to increase the efficiency of our simulations. (10) Duane, S.; Kennedy, A. D.; Pendelton, B. J.; Roweth, D. Phys. Lett. B 1987, 195, 216. (11) Nath, S. K.; de Pablo, J. J.; DeBellis, A. D. J. Am. Chem. Soc., in press. (12) Malewski, M. K. F.; Sandler, S. I. J. Chem. Eng. Data 1989, 34, 424. (13) de Loos, Th. W.; Poot, W.; Lichtenthaler, R. N. Ber. Bunsen-Ges. Phys. Chem. 1984, 88, 855. (14) Li, N. N.; Long, R. B. AIChE J. 1969, 15, 73. (15) Rogers, C. E.; Stannett, V.; Szwarc, M. J. Phys. Chem. 1959, 63, 1406. (16) Doong, S. J.; Ho, W. S. W. Ind. Eng. Chem. Res. 1991, 30, 1351.
