Integral Equation Theory for Gas Sorption and Swelling of Glassy

4914

Ind. Eng. Chem. Res. 2010, 49, 4914–4922

Integral Equation Theory for Gas Sorption and Swelling of Glassy Atactic Polystyrene Qinzhi Xu, Jianguo Mi,* and Chongli Zhong Laboratory of Computational Chemistry, Department of Chemical Engineering, Beijing UniVersity of Chemical Technology, Beijing 100029, China

On the basis of the polymer reference interaction site model (PRISM), a theoretical approach was proposed to describe gas solubilities in atactic polystyrene (aPS) melt. In the PRISM, our recently developed multisite semiflexible model was applied to derive intramolecular correlation functions. The intermolecular correlation functions for gas-gas and gas-aPS were obtained and further applied to calculate the excess chemical potential of gas-aPS system. The solubility coefficients of a wide range of gases in aPS were predicted at low pressure. These predictions are in good agreement with the experimental data. Particularly, the sorption curves of CO2 in aPS under high pressures were calculated, and the reasonable results were obtained. These results demonstrate that the theoretical model is reliable to describe the sorption properties of small molecules in aPS chain. On the basis of the CO2 sorption isotherms, the swelling of aPS under high pressure were also discussed. 1. Introduction Gases dissolved in amorphous polymers affect the thermophysical and mechanical properties of the polymer. Gas solubility parameter plays an important role in applications such as designing polymer barrier materials for packaging applications, developing membranes for gas separations, foaming, and plasticization.1-4 Knowledge of the solubility of small molecules in polymers is essential for design and operation of polymers and would help in understanding, predicting, and ultimately controlling the relations between structure, properties, processing, and performance of bulk phases, as well as of films, blends, nanocomposites, and other products based on them.5-7 To obtain solubility data, experimental measurements are still the dominant method to date. A great number of experiments have been undertaken in the past decades to investigate the sorption of various volatiles in polymers. For example, supercritical CO2 is widely used as a plasticizer due to its mild critical properties, whereas atactic polystyrene (aPS), which is glassy under ambient conditions, is one of the most widely used amorphous polymers. A very recent investigation is that of Panayiotou and collaborators,8 who measured the sorption of CO2 in aPS by using two experimental techniques, that is, the quartz crystal microbalance and the mass-loss analysis, over a wide spectrum of temperatures and up to very high pressures. Thus obtained solubility data are important to the design and operation of polymer plants, either for polymer plasticization, or for the removal of residual monomers, oligomers, and polymerization solvents from the polymer products. Experiments provide the practical data directly but cannot predict new data in varied conditions. Molecular simulation has become a very useful and accurate method to estimate sorption properties of a wide range of materials. During the past few years, modeling at the atomistic level, either molecular dynamic9-12 or Monte Carlo simulation,12,13 has been used to study mainly structural and some short time local dynamical aspects of gas solubility. Both the all-atom and united-atom simulations usually provide exact results in modeling to compare with the experimental data. However, they are quite time-consuming due to the complicated structure of * To whom correspondence should be addressed. E-mail: [email protected].

polymer systems. Recently, some more efficient coarse-grained (CG) models have been presented in the literature.14,15 These models vary in the degree of coarse graining as well as in the procedure for obtaining the effective interactions between the CG beads. Very recently, van der Vegt et al.16 proposed a CG model for aPS in which each aPS monomer is represented by two superatoms. The CG model has been further developed by employing a structure-based CG methodology that combines atomistic and CG simulations. The model can describe aPS sequences with varying tacticities and has been tested and validated for a number of structural properties of aPS. Furthermore, the CG model was used to predict melt properties, including the melt packing and the density.17 Coarse-grained models permit fast equilibration of polymer conformations in the melt as opposed to detailed-atomistic models. However, to compute gas solubility, CG models are not suitable. But this problem can be overcome by a (computationally cheap) inversemapping procedure that reintroduces atomistic degrees of freedom in the equilibrated, CG melt configurations.18 Theoretical investigation is always necessary for mechanism explanation, property prediction, and practical application. Previous workers applied the well-known Flory-Huggins theory to the sorption of gases in polymers.19 The theory reproduces qualitative features such as the shape of the sorption isotherm. Good fits to experimental solubility data at a given temperature can be achieved by the proper choice of an empirical Flory-Huggins parameter. However, the theory neglects the effect of polymer chain architecture on the intermolecular packing, therefore limiting its usefulness as a predictive tool.20,21 To describe gas solubility in polymers, the statistical mechanicsbased models are often used to calculate the sorption and swelling properties in the rubbery and glassy sate, such as the perturbed-hard chain theory,22 the lattice-fluid theory,23-25 the nonrandom hydrogen-bonding (NRHB) model,8 and the polymer reference interaction site model (PRISM).26,27 Compared with the three former theories, the PRISM can also provide the structure information of the system. The solubility coefficients of a wide range of monatomic gases in polyethylene were calculated in polymers with more complex structure.28 As pointed out by Curro,29 a related application of PRISM theory is to consider the solubility of a cavity particle in a polymer liquid. The cavity particle interacts with the polymer but has

10.1021/ie901728p  2010 American Chemical Society Published on Web 04/26/2010

Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010

no interaction with other cavity particles. From the chemical potential of the cavity particle as a function of the particle diameter, a measure of the free volume distribution of the polymer liquid can be deduced.30 Thus, the PRISM can be employed to probe the free volume distribution of various polymers as a function of chain architecture and backbone stiffness. However, this spherical approximation is only reasonable for spherical molecules such as argon and methane, whereas for nonspherical molecules, such as CO2 and chainlike alkanes, the overlap of atoms and the formation of molecule chain were neglected. To investigate more complicated polymers and gas-polymer systems, the intramolecular correlation functions are usually described by the flexible chain model or by molecular simulation. The flexible chain model neglects the fine architecture of polymer chains, such as bond lengths, bond angles, dihedral angles, and bond potentials, etc. Thus obtained intramolecular correlation functions deviate from real values systematically and can only give a qualitative description of chains. If they are applied to analyze phase equilibria, the deviations are inevitably exaggerated in Helmholtz free energy and chemical potential calculations. Although molecular simulation can produce accurate intramolecular correlation functions, the external input to the PRISM deteriorates the consistency of the theory itself. An exact theoretical method that can predict both the structure and the solubility of gas dissolved in complicated polymer systems is still worthwhile to explore. In this work, we extended our recently developed multisite semiflexible chain model31 for PRISM to investigate the solubility coefficients of spherical and nonspherical gases in aPS at low pressure, and sorption isotherms of carbon dioxide under high pressures. To establish a relatively accurate model to calculate the excess Helmholtz free energy and chemical potential of the gas-aPS systems, the real structure information of both gas and polymer molecules, including bond lengths, bond angles, and bond potentials, is included in the derivation of the intramolecular correlation functions. Furthermore, by introducing the potential of mean force (PMF), the original Lennard-Jones (LJ) force field for nonbonded interaction is corrected. With the two improvements, the intermolecular correlation functions of aPS were proven to be much more precise than those from the original PRISM with a flexible chain as input.31 The similar intramolecular correlation functions can also be obtained for nonspherical gas molecules. Accordingly, the energy and chemical potential of the system can be calculated. In addition, the swelling property of aPS caused by high solubility of carbon dioxide under high pressures was also discussed. 2. Theory 2.1. PRISM. For the purpose of the calculation of gas solubilities, the distribution function is the key issue and can be obtained by the integral equation theory of molecular liquids. Here we employ the PRISM to describe the intermolecular packing of the gas-polymer system. Since it has been extensively discussed in previous publications,28,29 we just briefly summarize the formalism as applied to aPS and gas solute with the multisite model. As for the system we work on, it contains small molecules in equilibrium with aPS having repeat units per chain. Each monomer structure of aPS is made up of eight interaction sites. For gas solute, taking CO2 as an example, each molecule contains three sites. The schematic structure of aPS and CO2 is shown in Figure 1. In part (a), A and B represent the CH2 and CH in the backbone, and C, D, E, F, G, and H correspond to the C atom and CH in the side phenyl ring group,

4915

Figure 1. Schematic representation of aPS (a) and CO2 (b)

respectively. In part (b), the three sites denote O, C, and O atoms, respectively. The distribution function of primary interest is the intermolecular correlation functions gRMγM′(r) for gas-gas or gas-polymer. In general, for molecular fluids, it can be used to characterize the intermolecular packing, which is defined as follows32 FRMFβM′gRMγM′(r) ) 〈

∑ δ(br

RM r i )δ(b

-b r γj M′)〉

(1)

i*j)1

where FRM is the site density of type R in molecules of species M, and b riRM is the position vector of site R on chain i in molecule M. Following Chandler and Andersen’s definition28,29 of integral equations to multicomponent systems, the relationship between intermolecular correlation functions and direct correlation functions can be further expressed as the following matrix form in Fourier space ˆ (k)C ˆ (k)[Ω ˆ (k) + G˜ hˆ(k)] hˆ(k) ) Ω

(2)

where G˜ denotes the density matrix with elements FRMγM′ ) δRMγM′FM, FM is the monomer number density of molecular ˆ (k) accounts for species M, δRMγM′ is the kronecker delta, C the matrix of direct correlation functions, hˆ(k) stands the matrix of total correlation functions and its element hRMγM′(r) in r space is defined according to hRMγM′(r) ) gRMγM′(r) - 1, ˆ (k) represents intramolecular correlation functions, and Ω which can be given by the following Fourier transform ˆ R γ (k) ) Ω M M′

δRMγM′ NRM

∑ ∑ ωˆ (k) ij

(3)

i∈RM j∈γM′

where ω ˆ RMγM′(k) is the normalized probability density between RM and γM′ sites on a single chain or molecule of type RM. To describe different properties of units clearly, the matrices in eq 2 can be further divided into four submatrices representing of gas-gas, gas-polymer and polymer-polymer interactions, respectively hˆ(k) )

ˆ (k) ) Ω

ˆ (k) ) C

( ( (

)

(4a)

ˆ 11(k) (Ω 0 ˆ 22(k) Ω 0

(4b)

ˆ 11(k) C ˆ 12(k) C ˆ 21(k) C ˆ 22(k) C

(4c)

hˆ11(k) hˆ12(k) hˆ21(k) hˆ22(k)

G˜ )

(

G˜ 11 0 0 G˜ 22

)

) )

(4d)

Here the subscript 11 represents the information of the monomer sites of polymer, 22 corresponds to gas molecule sites,

4916

Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010

and the intersect terms 12 and 21 describe the relationship between gas and polymer. When applied to gas-aPS system, ˆ 11 are an 8 × 8 matrix for aPS; hˆ22, Ω ˆ 22, and ˆ 11, and C hˆ11, Ω ˆ 22 are a 3 × 3 matrix for CO2 or propane, and can be simplified C to a 1 × 1 matrix for a spherical gas solute. To solve eq 2 for the considered system, the Percus-Yevicklike closure33 has been used

a previous successful run with similar parameters. For simplicity, the cutoffs for the LJ potential and the PMF are the same, 2.25σ. 2.2. Molecular Model of Gas Solubility. It has been demonstrated that the PRISM can be used to analyze the solubility of monatomic gases in polymer28,29 by calculating the chemical potentials of both gas and gas in polymer. For pure gas, the chemical potential can be easily obtained by

CRMγM′(r) ) gRMγM′(r)[1 - eβuRMγM′(r)]

βµAG ) ln[βPΛ3]

(5)

in which β ) 1/kBT, uRMγM′(r) is the LJ force field to describe the interactions between nonbonded sites

[( ) ( ) ] σRMγM′

uRMγM′(r) ) 4εRMγM′

12

-

r

σRMγM′ r

6

(6)

where εRMγM′ and σRMγM′ are the site parameters, and the interactions between polymer and gas are determined by Lorentz-Berthelot mixing rules: σRMγM′(r) ) (σRMRM′ + σγMγM′)/2

(7)

εRMγM′(r) ) (εRMRM′εγMγM′)1/2

(8)

According to Schweizer,34 the PMF W(r) is usually introduced to modify the correlation functions. W(r) is a medium-induced potential that reflects the effects of the intermolecular chains in the system, and can be obtained by35 βWRγ(k) ) -

∑ Cˆ

ˆ

ˆ

˜ Sij(k)Cjγ(k) Ri(k)F

Here Sˆ(k) is the structure factor and defined in elsewhere.36 W(r) corresponds to both intra- and intermolecular interactions, and combines them together, leading to a progress on the theoretical self-consistence. Following Li et al.,35 we apply W(r) to modify the LJ potential u(r) ) uLJ(r) + W(r)

where P is the pressure and Λ is the de Broglie wavelength. For gas A dissolved in polymer B at low pressure, the chemical potential of A can be decomposed into three distinct parts: the ideal term µ0, the hard-sphere reference term µref, and the attractive contribution µatt; therefore, an equation of state can be written as µAP ) µ0 + µref + µatt

(10)

In the whole calculation of the PRISM, the intramolecular correlation functions of aPS are the most important input. In our previous work, they were formulated with a six-site semiflexible chain model by combining the generator matrix method38 and Koyama distribution function,37 in which the two sites G and H are approximately represented by E and F due to the symmetry structure of the benzene. Therefore, 21 independent correlation functions were derived out to stand for eightsite interactions.31 For gas solutes, the corresponding intramolecular correlation functions can be easily obtained by the same way. With these functions, eq 2 can be solved numerically by using simple Picard iteration method. In the process of the PRISM calculation, we first propose the initial guess for the indirect correlation function γ ) h C, and then iterate on eq 2 using the damped iteration scheme γ )(1 - R)γnew + Rγold, where R