Density-Functional Theory for PolymerâCarbon Dioxide Mixtures

Article pubs.acs.org/IECR

Density-Functional Theory for Polymer−Carbon Dioxide Mixtures Xiaofei Xu,† Diego E. Cristancho,‡ Stéphane Costeux,‡ and Zhen-Gang Wang†,* †

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, United States The Dow Chemical Company, Midland, Michigan 48674, United States

‡

ABSTRACT: We propose a new density-functional theory (DFT) describing inhomogeneous polymer−carbon dioxide (CO2) mixtures. The theory is constructed by combining the bulk Peng−Robinson equation of state (PR-EOS) with the statistical associating fluid theory (SAFT) and the fundamental measure theory (FMT). The weight density functions from FMT are used to extend the bulk excess Helmholtz free energy of PR-EOS to the inhomogeneous case, while the SAFT is used to describe correlations due to polymer chain connectivity and short-range forces due to weakly polar or association interactions. The additional long-range dispersion contributions are included using a mean-field approach. We apply our DFT to the interfacial properties of polystyrene−CO2 and poly(methyl methacrylate)−CO2 systems. The calculated interfacial tension values are in good agreement with experimental data at low to intermediate pressures. The inclusion of association energy for CO2 is shown to have a significant effect. We also point out the limitation of the PR-EOS for describing polymer−CO2 mixtures at high pressures (P > 35 MPa).

I. INTRODUCTION The phase behavior of polymers in condensed carbon dioxide (CO2) is of tremendous practical importance in industry due to the widespread use of CO2 as an environmentally benign solvent.1−4 For instance, the fabrication of novel foam materials with nanoscale pores requires a systematic and quantitative study of phase separation.5 Despite numerous advances on the experimental front,6−9 the modeling of polymer and CO2 mixtures is still a challenge owing to the lack of theoretical models that can accurately describe the binary system at the molecular level.4,10 For practical applications, CO2 is usually in the supercritical state. Therefore, a small change in control parameters (such as pressure or temperature) results in a large variation of the physical properties (such as density) of the system. Compressibility of the system must be properly taken into account. The appearance of both liquid−liquid and vapor−liquid coexistence, together with the large disparity in the molecular weight between the polymer and CO2, makes the theoretical description of polymer−CO2 mixtures highly nontrivial. Previous efforts are primarily based on molecular simulation and self-consistent field theory (SCFT). In the classification of van Konynenburg and Scott,11 the phase behavior of a short polymer chain solution is usually of Type I or Type II. In this case, the critical points of pure polymer and solvent are connected by the liquid−vapor critical line in the pressure− temperature plane. For long polymer chains, the phase behavior is usually of Type III, where not only liquid−vapor but also liquid−liquid phase separation occurs. Computer simulations have been primarily used to study Type I and II phase behavior,12−15 while SCFTs have been used for Type III.16,17 By using SCFT, Müller et al.16 studied the phase separation of Type III. At low pressures, the interfacial properties and nucleation behavior are similar to a one-component fluid. At higher pressures, there is a triple point at which the polymer-rich phase coexists with a vapor of the solvent and a solvent-rich liquid. Park et al.17 studied the effect of temperature and pressure © 2012 American Chemical Society

on the surface tension of polystyrene in supercritical CO2 by SCFT. They found that the surface tension varies linearly with temperature and with pressure. Virnau et al.12 pointed out that the presence of the triple point has a profound effect on the interfacial properties and the early stage of nucleation behavior. In particular, there is a discontinuity in the slope of the interfacial tension versus pressure across the triple point. MacDowell et al.13 investigated the critical lines by thermodynamic perturbation theory of first order (TPT1) and computer simulation. Their results indicate that TPT1 works well at low pressures, while at higher pressures it overestimates the pressure in critical lines for very asymmetric mixtures. Virnau et al.14 studied the phase behavior of n-alkanes in supercritical solution by using Monte Carlo simulation. Their results show that a coarse grained model could reproduce experimental results quite accurately. Mognetti et al.15 investigated the effect of dipolar interaction on the phase diagram of polymer solutions by using a spherically averaged quadrupolar potential. These investigations reveal that it is possible to develop a simplified and efficient coarse-grained model for polymer− CO2 mixtures. However, most of these studies are limited to polymers with rather short chain lengths of 5 to 10, which is far shorter than typical polymer lengths in practice of application. Although computer simulation in principle can give accurate results, the computation time becomes prohibitively long for long polymers.18 SCFT is less demanding computationally, but it does not provide structural information at the molecular level. DFT is an attractive alternative as it can successfully capture the necessary molecular level details while requiring less computational effort. The free energy in DFT of a given system is represented as a functional of components densities. Received: Revised: Accepted: Published: 3832

December 13, 2011 February 7, 2012 February 9, 2012 February 9, 2012 dx.doi.org/10.1021/ie2029267 | Ind. Eng. Chem. Res. 2012, 51, 3832−3840

Industrial & Engineering Chemistry Research

Article

atom contributes to one attractive site, and each oxygen atom contributes to one attractive site. We note that the interaction between the CO2 molecules can alternatively be modeled by quadrupoles.23 We adopt the association model of Shin and Wu because we use the parametrization of the other parameters in our DFT as in their equation of state. A polymer is modeled as a freely jointed chain with N identical segments of diameter σ2 tangentially connected. Mathematically, chain connectivity is enforced by the bonding potential between nearest-neighbor segments,

Equilibrium density profiles are then obtained by minimizing the grand potential, which yields the equilibrium properties of the system such as interfacial tension and surface excess. Talreja et al.19 investigated the interfacial properties of thin polymer films in the presence of CO2 by DFT. However, their functional is based on a crude weight function and hence the resulting DFT does not provide an accurate description of the structural details at the molecular level. The structural properties are usually quite sensitive to the form of the weight function.20 Li and Firoozabadi21 developed a DFT to study the interfacial tension of binary mixtures of small molecules based on the Peng−Robinson equation of state (PR-EOS). Although their theory yields good results for the interfacial tension over a wide range of conditions, the theory is limited to small molecules with no chain connectivity. Furthermore, as will be demonstrated in our work, the association interaction between CO2 molecules plays a significant role. The purpose of this work is to develop an accurate and an efficient DFT for polymer−CO2 mixtures that includes most of the essential interaction and correlation effects and is valid for long polymers having thousands of monomers. We develop our DFT from the PR-SAFT EOS of Shin and Wu;22 this EOS has been shown to accurately describe the bulk thermodynamic properties of polymer−CO2 mixtures at low to moderately high pressures. The DFT is constructed by combining the bulk PR-EOS with the statistical associating fluid theory (SAFT) and the fundamental measure theory (FMT). We then apply our DFT to examine the phase behavior and interfacial properties of polystyrene (PS)−CO2 and poly(methyl methacrylate) (PMMA)−CO2 mixtures.

exp[−βVB(r N)] =

N−1

δ(|ri + 1 − ri| − σ2) 4πσ2 2 i=1

∏

(1)

where r = (r1, r2, ..., rN) and δ is the Dirac delta function, and σ2 is the diameter of polymer segments. β−1 = kBT stands for the temperature multiplied by the Boltzmann constant. The interaction between two arbitrary species (i.e., polymer segments or CO2) is described by N

⎧∞ r < σij ⎪ uij(r ) = ⎨ ⎪−εij(σij /r )6 r ≥ σij ⎩

(2)

The constants εij and σij (i, j = 1,2, with 1 denoting CO2 and 2 denoting the monomer of polymer) are energy parameters and the effective diameter of interaction, respectively. The effective diameter of a monomer is related to the size parameter bi as σi3 = 3bi/(2πNA), where NA is the Avogadro number. The energy parameter is determined by εi = 3ai/(2πσi3NA2). The parameters ai and bi are obtained by fitting the volumetric properties of pure polymer and the binary vapor−liquid equilibrium data for CO2−polymer mixtures.22 Because a in the PR-EOS is an empirical parameter obtained by best fitting, it is generally temperature-dependent, thus making ε also temperaturedependent.10,22 We therefore take it to be the average value in the range of T = 260 K to T = 450 K. The overall variation is on the order of 1% or less, which leads to a corresponding range of variation for the results. The Lorentz− Berthelot mixing rule is used to obtain the cross interaction parameters: σij = (σi + σj)/2, εij = (εiεj)1/2. The chain length (N) of the polymer is linearly proportional to the molecular weight (M). The ratio N/M for each polymer is obtained in the work of Shin and Wu22 by fitting the equilibrium EOS data for the polymers studied. In comparing our results with experimental data, we use the same ratios as determined by Shin and Wu and the actual molecular weight M of the polymer to determine the chain length N. In Table 1, we list the value of

II. MODEL AND THEORY A. Molecular Model. We consider a compressible polymer− CO2 mixture, where the molecular units of both components are coarse-grained as spherical particles with a hard core. The CO2 molecule is modeled as a particle of diameter σ1 with three binding sites, accounting for the different attractive sites of carbon and oxygen atoms.22 As is shown in Figure 1, the carbon

Table 1. Parameters for the Species Studied in This Work

Figure 1. Schematic of associating CO2 molecules. The association interaction arises from the positive and negative partial charges on the C and O atoms, respectively. Following ref 22, we capture this interaction by assigning three interaction sites, two corresponding to the O atoms (cyan solid circles) and one to the C atoms (red solid circles). However, it should be emphasized that the location of the association site is not to be taken literally. The association site is used merely to indicate the saturated nature (due to possible shift in partial charge of the C atom upon interacting with an O atom from another molecule) of the association interaction.

species

M [g/mol]

N [−]

ε/kB [K]

σ [A]

CO2 PMMA PS

44 89 230 150 000

1 2422 1915

1035 6698 6674

2.74716 2.92616 2.92709

the molecular weight (M), chain length (N), energy parameter (ε), and monomer diameter (σ) for the species used in our model. B. Helmholtz Free Energy Functional. The Helmholtz energy functional F is expressed as a summation of an ideal term corresponding to the system of an ideal gas of monomers and polymers free of nonbonded interaction, and an excess 3833

dx.doi.org/10.1021/ie2029267 | Ind. Eng. Chem. Res. 2012, 51, 3832−3840


Article

represents the geometrical properties of molecules. Based on these weight functions, the weighted densities are defined as

contribution accounting for the inter- and intramolecular interactions, F[ρ1(r), ρ̂2(r N)] = F id[ρ1(r), ρ̂2(r N)] + F ex [ρ1(r), ρ̂2(r N)]

N

(4)

i=1

∫ dr ρ1(r)[ln(ρ1(r) − 1] + ∫ dr N ρ̂2(r N)[ln ρ̂2(r N) − 1] + β ∫ dr N ρ̂2(r N) VB(r N)

(5)

ex + F ex + F ex + F ex F ex = Fpr assoc disp chain

∫ dr′ ρi(r′) Θ(σi/2 − |r − r′|)

n V, i(r) =

∫ dr′ ρi(r′) δ(σi/2 − |r − r′|) |rr −− rr′|′

aρ 1 + (1 + ln 2 2 b 1 + (1 −

n0, i ζii ln(1 − 4n3) − 1 + 4(1 + 1 + 4(1 −

1 8 2 n3

2 )n3 2 )n3

(9)

with ζij = 1 − nV,i·nV,j and n3 = n3,1 + n3,2.We note that in the DFT for small molecules developed by Li and Firoozabadi,21 which is also based on the PR-EOS, the vector density function nV,i was not included. This function is necessary for capturing the molecular level structure at high densities.24,25 The WDA alone does not sufficiently describe the long-range intermolecular attractions.26 The additional dispersion contributions are included in a mean-field manner by27 ex [ρ (r), ρ (r)] Fdisp 1 2

ϕpr is the reduced excess Helmholtz free energy density and in bulk phase it has the expression ϕpr = − ∑ ρibulk ln(1 − bρbulk ) −

n0, in0, jaij ζij ln

i , j = 1,2

(7)

bulk

∑

×

(6)

∫ dr ϕpr[ρ1(r), ρ2(r)]

∑ i = 1,2

Fprex accounts for the short-range repulsion and partially for the long-range attraction, which is built on the PR-EOS, given by

i

n3, i(r) =

ϕpr[ρ1(r), ρ2(r)] = −

The excess Helmholtz free energy is assumed to be a functional of segmental density distributions. It includes the contribution from short-range interactions, long-range dispersion interactions, and associating interaction among CO2 molecules and correlations due to the intramolecular polymer chain connectivity,

ex βFpr =

∫ dr′ ρi(r′) δ(σi/2 − |r − r′|)

where Θ is the Heaviside step function. n0,i is the number density of monomers, n2,i describes the effect from geodesic curvature of density variation,25 and n3,i is the packing fraction of monomers for component i. The vector function nV,i is needed to describe the density gradient, which vanishes in uniform density profiles. In terms of these weighted densities, ϕpr can be expressed as

The ideal term of Helmholtz free energy is known exactly as βF id[ρ1(r), ρ̂2(r N)] =

∫ dr′ ρi(r′) δ(σi/2 − |r − r′|)

n2, i(r) =

N

∑ ρ2, i(r) = ∑ ∫ dr N δ(r − ri)ρ̂2(r N) i=1

πσi 2

(3)

ρ1(r) is the density profile of CO2. ρ̂2(rN) is the multidimensional density profile of polymer chain, which is related to the segmental densities ρ2,i(r), (i = 1, 2 ,..., N) by ρ2(r) =

1

n0, i(r) =

=

bulk

2 )bρ

bulk

2 )bρ

1 4

∑

∫ ∫ dr dr′gij(|r − r′|)uij(|r − r′|)

i , j = 1,2

× [ρi(r) − ρi(r′)][ρj(r) − ρj(r′)]

(8)

(10)

ex This form ensures that Fdisp vanishes for the uniform bulk state and hence does not over account the contribution included in the bulk EOS. The pair distribution function gij should in principle be evaluated at some reference liquid density.28 However, the choice of a good reference density is not obvious for our system since there is a large change in density and composition across the liquid−vapor interface. Thus we make the simplest mean-field approximation by writing gij(r) = Θ(r − σij), which corresponds to the limit of zero density in the reference fluid. Because of the long-range nature of the dispersion interaction, we expect that this approximation, which essentially amounts to averaging out oscillations in the gij, should be a reasonable one. ex The Helmholtz energy Fassoc due to intermolecular and intramolecular association of CO2 molecules is29,30

where ρbulk = ∑iρi bulk is the total number density in bulk phase; a and b are the energy and volume parameters of the mixture, which can be estimated by applying the mixing rules of a = ∑i,jxixjaij, b = ∑ixibi and aij = (aiaj)1/2. xi, ai, and bi stand for the mole fraction, energy parameter, and volume parameter of pure component i, respectively. The central idea in FMT is to describe the short-range interaction by a set of weight density functions, which arise from a geometrical decomposition of the excluded volume between two hard spheres.24 Here we use the same decomposition, taking the point of view that the PR EOS is an approximate EOS for effective hard spheres with attractive interactions. This approach is not a true first-principles construction of the DFT, but rather a practical way of extending a simple and relatively accurate bulk EOS to describe some coarse-grained features of the inhomogeneous polymer−CO2 mixtures. We adopt weighted density approximation (WDA) to extend ϕpr to the inhomogeneous states. Four weight functions from Rosenfeld FMT24 are introduced, which mathematically

ex βFassoc =

3 ⎡

∫ dr n01(r) ∑ ⎢⎣ln Xl(r) − l=1

Xl(r) 1⎤ + ⎥ 2 2⎦ (11)

3834



Article

reveals that the end effect becomes insignificant as the segment index increases for sufficiently large i, that is, Ii(r) ≈ I(r) becomes independent of the segment index. Therefore, for a long homopolymer chains, we make the approximation Ii(r) ≈ I100(r) for i > 100. The end effect at the other end of chain is described by the symmetric expression of the eq 17. This simplification considerably accelerates the convergence of the iterations and enables us to investigate systems with very long polymer chains. To determine the density profile across an interface ρi(z), accurate values of the coexisting density {ρiL,ρiV|, i = 1, 2} are required, where ρiL and ρiV are the equilibrium density for component i in the liquid and vapor side, respectively. These are determined by searching the chemical potentials such that μi(ρ1L, ρ2L) = μi(ρ1V, ρ2V), i = 1, 2 at given temperature and pressure. Using {ρi L,ρiV} as an input, the equilibrium density profiles of interface are obtained by iterating the Euler− Lagrange equations. The calculations are performed within the range 0 ≤ z ≤ L with the Dirichlet boundary

with ⎡ Xl(r) = ⎢1 + ∑ Δlm(r) ⎢ ⎣ m

∫

⎤−1 δ(|r − r′| − σ1) ⎥ dr′ Xj(r′)ρ1(r′) ⎥ 4πσ12 ⎦

hs Δlm(r) = y11 (σ1, nα)[exp(βelm) − 1]klm hs (σ , n ) = y11 1 α

n2 σ1ζ11 n22 σ12ζ11 1 + + 1 − n3 4(1 − n3)2 72(1 − n3)3

where Xl(r) stands for the mole fraction of molecules at r not bonded at site l; elm is the bonding energy between site l and m, and klm is a parameter characterizing the accessible volume for binding. The value of these parameters can be found in ref 22. ex Finally, Fchain due to chain connectivity of the polymer is given 31 by ex = βFchain

1−N N

∫ dr n02ζ22 ln y22 (σ2, nα)

ρi(0) = ρiV ,

(12)

with hs (σ , n ) = y22 2 α

n2 σ2ζ22 n2 2 σ2 2ζ22 1 + + 1 − n3 4(1 − n3)2 72(1 − n3)3

Ω[ρ1(r), ρ̂2(r)]

∫ ρ1(r) dr − μ2 ∫ ρ̂2(r N)d r N

(13)

where μ1 and μ2 are the chemical potential of CO2 and the polymer chain, respectively. At equilibrium, the grand potential is minimized, δΩ δΩ = =0 δρ1(r) δρ̂2(r N)

i = 1, 2

In most cases, we assumed that L = 20σ1 and used a grid of the size Δz = 0.01σ1 for the computation. However, nearby the critical point, we use L = 30σ1 and Δz = 0.005σ1 to ensure the accuracy. The equations are solved using Picard iteration with a convergence criterion that the deviation of Euclidean distance between two consecutive density profiles is less than 10−6. During the iterations, we may encounter the floating problem for the density profiles.32 This is because the location of the interface is translationally invariant, that is, if {ρi(z)} is a solution of the density profile, then the function of {ρi(z + d)} also is a solution, where d is any arbitrary shift. To avoid the floating problem, we maintain fixed the position of the Gibbs dividing surface in the following manner. At the beginning of the iteration process, the position of the Gibbs dividing surface is at z = L/2 with the starting profile:

C. Euler−Lagrange Equations. The grand potential Ω of the system is related to the Helmholtz free energy functional as

= F − μ1

ρi(L) = ρiL ,

⎧ρV , ⎪ i ρi(z) = ⎨ ⎪ρL , ⎩ i

(14)

Combining eqs 3, 4, 13, and 14 yields the Euler−Lagrange equations

0 ≤ z < L /2 L /2 ≤ z ≤ L

During the iterations, the entire profiles are shifted left or right, depending on the sign of the difference

⎡ δβF ex ⎤ ⎥ ρ1(r) = exp⎢βμ1 − ⎢⎣ δρ1(r) ⎥⎦

(15)

N ⎡ δβF ex ⎤ ⎥ ∑ Ii(r)IN + 1 − i(r) ρ2(r) = exp⎢βμ 2 − ⎢⎣ δρ2(r) ⎥⎦ i=1

(16)

where ρ(z) = ρ1(z) + ρ2(z) is the total density profile. In practice, the position of the Gibbs dividing surface is adjusted every 10 iterations.

(17)

III. RESULTS AND DISCUSSION We first examine the bulk phase behavior, in order to establish the range of applicability of our PR-SAFT EOS based DFT. In typical polymer solutions with supercritical CO2, both the liquid−vapor and liquid−liquid phase equilibria are present in the phase diagram, leading to the occurrence of a triple point. We thus focus the region of the phase diagram near the triplepoint condition. In Figure 2, we present the density−density phase diagram for PS−CO2 mixtures at various temperatures T = 300, 310, 320, and 373.15 K. For T = 320 K and T = 373.15 K, we find coexistence between a nearly pure PS-rich phase and a low density CO2-vapor, which ends at high pressure. For T = 300 and 310 K, there is a triple point at P = 6.45 MPa and

∫0

where the recursive function is given by Ii(r) ⎧1, i=1 ⎪ ⎪ ⎡ δβF ex ⎤ δ(σ − |r − r′|) =⎨ 2 ⎥ Ii − 1(r′), i > 1 ⎪ dr′exp⎢ − ⎪ ⎢⎣ δρ2(r′) ⎥⎦ 4πσ2 2 ⎩

∫

Equations 14, 15, and 16 provide the key equations for calculating the density profiles of polymer and CO2 mixtures. D. Numerical Procedure. Because the propagator function Ii(r) is defined in terms of individual polymer segments, the numerical iteration becomes computationally intensive for a long polymer chain. Solution of the recursive relation (eq 17) 3835

L /2

ρ(z) dz −

L

∫L /2 ρ(z) dz



Article

Figure 2. Density-density phase diagram of PS−CO2 mixtures at (a) T = 300 K, (b) T = 310 K, (c) T = 320 K, and (d) T = 373.15 K. The red lines are the spinodal curves and the blue lines are the binodal curves. The dashed lines are the tie lines of phase coexistence. The triple points are marked by squares.

P = 7.73 MPa, respectively (the triple-point pressure), where a PS-rich phase coexists with an almost pure CO2-vapor and a CO2-rich liquid-like phase. Below the triple pressure, the PS-rich phase coexists with the CO2-rich vapor, while above the triple pressure it coexists with a CO2-rich liquid-like phase. In addition to the equilibrium phase boundary (the binodal), we present the spinodal curve in Figure 2. It can be seen that the binodal curves terminate at an end point at high pressure. Above that pressure, the spinodal curves extend beyond the region enclosed by binodal curves; such a behavior is clearly unphysical. This unphysical behavior is also reported by van Konynenburg and Scott11 for the van der Waals EOS. Some peculiar behaviors of third-order polynomial EOS at high densities have also been reported in the work of Müller et al.16 Our ongoing work indicates that this unphysical behavior does not appear in a phase diagram calculated in using perturbedchain SAFT EOS,38 which better captures the local packing effects at high densities. Moreover, experiments on CO2 and hexadecane show that the binodal and spinodal curves terminate in a critical point.33 We thus conclude that the PR EOS or any EOS represented by a third-order polynomial, and hence DFT, based on these EOS, are not reliable at very high pressures. For our PR EOS based DFT, on the basis of our phase diagram data and data from the work by Shin and Wu, we estimate the pressure beyond which the theory becomes unreliable (for example, the deviation from experimental data on the solubility becomes greater than 10%) to be about 35 MPa. To further investigate the phase behavior, we examine the phase boundaries in the P−T plane. Figure 3 shows the coexistence lines of pure components and the two critical lines for PS−CO2 mixtures in P−T plane. The results reveal that the phase behavior is of Type III in the classification of van Konynenburg and Scott.11 The critical line I emerges from the critical point of pure CO2 and terminates at an end point. There it connects to a triple line, on which a PS-rich phase coexists with a CO2-rich vapor and a CO2-rich liquid-like phase. As the temperature decreases, the triple line merges with the coexistence curve of pure CO2. The critical line II emerges from

Figure 3. Coexistence lines of the pure components and critical and triple lines of PS−CO2 mixtures in the P−T plane.

the critical point of the pure PS and moves up in pressure with decreasing temperature. Because of the inaccuracy of the PR EOS at pressure above 35 MPa, the portion of the critical line above this pressure is not quantitatively accurate, although qualitatively it has the expected behavior. Having obtained a global view of the phase behavior and established the range of validity for our PR-SAFT EOS based DFT, we now apply the theory to investigate the interfacial properties at the liquid−vapor and liquid−liquid coexistence. We first examine the interface at the liquid−vapor coexistence. Figure 4 shows the density profiles between the PS-rich liquid and the CO2-rich vapor at T = 373.15 K and P = 1 MPa. At the interface, there is a sharp change in the PS density and notable surface enhancement of CO2. Examination of the segmental contribution to overall PS density indicates that the middle segments tend to be confined to the interior of the interface (the PS-rich side) while the end segments are segregated closer to the vapor. If we neglected the contribution due to the longrange dispersion interaction, eq 10, the theory would predict significant density oscillations on the liquid side (the dashed lines), 3836



Article

practical applications. The interfacial tension at a planar interface is calculated from γ = (Ω − Ωbulk)/A, where Ωbulk is the grand potential in bulk phase and A is the surface area. Figure 6 shows the interfacial tension of the coexisting phase as

Figure 4. Density profiles across the coexisting interface for PS−CO2 mixtures at T = 373.15 K and P = 1 MPa. The inset shows the magnification of the density profile of CO2 in the liquid side. The dashed lines are the results in absence of long-range dispersion ex . interaction Fdisp

reflecting the local layering of the PS monomers against a nearly hard wall. The density oscillations at an interface have also been reported for pure fluids by computer simulations.36 The amplitude of these oscillations for a general liquid−vapor interface depends on the thermodynamic state and on the strength and range of the attractive interactions.37 Upon inclusion of the additional long-range dispersion interaction these oscillations are greatly reduced. To further explore this phenomenon, we show the pressure dependence of PS density profiles of the coexisting phase in Figure 5. As the pressure increases,

Figure 6. Interfacial tension between polymer-rich phase and solventrich phase for PS−CO2 at T = 310 K. The insets show the density profiles between coexisting phases near the triple pressure at (a) P = 7.72 MPa, (b) P = 7.73 MPa, and (c) P = 7.74 MPa. The triple-point pressure is 7.73 MPa.

a function of pressure at a supercritical temperature T = 310 K. At this temperature, increasing the pressure can take the systems across the triple line; see Figure 2. The overall decrease of the interfacial tension with the pressure is as expected, consistent with our consideration based on density contrast given earlier. At the triple pressure P = 7.73 MPa, there is a discontinuity in the slope of the interfacial tension. For P < 7.73 MPa, the PS-rich phase coexists with a CO2-vapor, while for P > 7.73 MPa it coexists with a CO2-liquid-like phase. As pressure further increases, the interfacial tension decreases and eventually vanishes at the liquid−liquid critical point. In the insets, we show the density profiles across the interface. At P = 7.72 MPa, the density of the CO2 is almost equal in both phases, with an excess of CO2 at the interface. At P = 7.73 and 7.74 MPa, the density of CO2 in the CO2-rich side (left side) increases due to the occurrence of CO2-liquid-like phase. In Figure 7, we compare our DFT predictions with the experimental results34 of interfacial tension for both PS−CO2 and PMMA−CO2 mixtures. At these conditions, the CO2-rich phases are vapor. In general, our theory satisfactorily reproduces the experimental results. In that figure, we also address the effects of including the CO2 association energy. The results show that the inclusion of this energy, which describes the intermolecular association of CO2, considerably improves the agreement with experiments. Although the association energy of CO2 is weak,38 the incorporation of association interaction is necessary in describing compressed polymer−CO2 binary mixtures. We note that the agreement for PMMA is not good as for PS. This may be due to the association-like interaction between CO2 and PMMA that is not considered in this work. It is of interest to study the distribution of CO2 in the interfacial region. At low pressures, the density of CO2 is higher in the PS liquid than in the vapor phase due to the favorable attraction between CO2 and PS. As the density of the (nearly pure CO2) vapor phase increases, the situation becomes

Figure 5. Density profile of PS across the coexisting interface for PS−CO2 mixtures at T = 373.15 K. The dashed lines are the results in ex . absence of long-range dispersion interaction Fdisp

the amplitude of these oscillations is also greatly reduced. When the density in the vapor side is very low, excursion of the PS segments into the vapor phase is energetically very costly, so the interface effectively acts as a hard wall for the PS. At higher pressures, the density contrast between the vapor and liquid becomes less, so fluctuation (including capillary waves) becomes less costly, thus smearing out the oscillations. Interfacial tension is one of the most important physicochemical properties for the polymer−CO2 mixtures for many 3837



Article

Figure 7. Interfacial tensions between polymer-rich phase and solvent-rich phase as a function of pressure. The points are the experimental data.34,35 For the molecular weights of PS and PMMA, we use the weight-averaged molecular weights in the experiments. The polydispersity index for PS and PMMA are Mw/Mn = 1.804 and 2.302, respectively. The results become insensitive to the molecular weights at high enough molecular weights (N > 1000).

with liquid CO2. At T = 320, 373.15, and 423.15 K, the excess adsorption also shows a sharp slope at low pressure. With further increasing pressure, they go through a broad maximum and drop continuously as there is no phase change. These properties are similar to that of CO2 adsorption on octacosane20 and other substances in the experiment.39,40

reversed. However, there is excess adsorption of CO2 at the interface for all pressures. Figure 8 shows the excess adsorption

IV. CONCLUSIONS We have developed a density-functional theory (DFT) for describing inhomogeneous polymer−carbon dioxide (CO2) mixtures. The functional is constructed from a PR-SAFT EOS.22 The weighted density functions from fundamental measure theory24 are adopted to extend the bulk excess Helmholtz free energy to the inhomogeneous cases. The statistical associating fluid theory is used to describe the polymer intrachain correlations and short-range forces due to weakly polar or association interactions. The additional long-range dispersion contributions are included by a mean-field approach. We have applied our DFT to the interfacial properties of polystyrene−CO2 and poly(methyl methacrylate)−CO2 systems. In the absence of long-range interaction, our DFT predicts oscillatory density profiles for the polymer segments at the liquid−vapor interface, which is not captured in some other versions of DFT.19 The amplitudes of these oscillations are greatly reduced upon adding the additional long-range dispersion interaction. The interfacial tension decreases with increasing pressure, and its slope shows a discontinuity at the triple point. Our calculated interfacial tension is in good agreement with experimental data for coexistence between a polymer-rich liquid and a CO2 vapor at low to moderately high pressures. We find that including the association energy of CO2 is necessary for satisfactorily reproducing the experimental data. As the parameters of the PR-SAFT EOS were optimized to reproduce bulk thermodynamic data at low to moderate pressures, our DFT based on this EOS is similarly valid in this pressure range, which is rather broad for many practical applications. At high pressures (P > 35 MPa), however, the PRSAFT EOS produces some unphysical features with regard to

Figure 8. Excess adsorption isotherms of CO2 on PS interface as a function of pressure.

isotherms of CO2 as a function of the pressure. The excess CO2 adsorption is defined as Γ=

z

∫0 G [ρ1(z) − ρ1V] dz + ∫z

L G

[ρ1(z) − ρ1L] dz

where zG is the Gibbs dividing surface of PS density profiles, obtained by solving the equation z

∫0 G [ρ2(z) − ρ2V]dz + ∫z

L G

[ρ2(z) − ρ2L] dz = 0

The results show that the excess adsorption first increases with the pressure, but decreases as the pressure is further increased. At T = 300 K and T = 310 K, the excess adsorption increases sharply with pressure at low pressures and drops sharply at the triple pressure. With a further increase in pressure, the excess adsorption decreases as the adsorbed layer is in equilibrium 3838



Article

simple fluids and polymer solutions Monte Carlo studies of the phase behavior. Phys. Chem. Chem. Phys. 2009, 11, 1923−1933. (16) Müller, M.; MacDowell, L. G.; Virnau, P.; Binder, K. Interface properties and bubble nucleation in compressible mixtures containing polymers. J. Chem. Phys. 2002, 117, 5480−5496. (17) Park, H.; Thompson, R.; Lanson, N.; Tzoganakis, C.; Park, C.; Chen, P. Effect of temperature and pressure on surface tension of polystyrene in supercritical carbon dioxide. J. Phys. Chem. B 2007, 111, 3859−3868. (18) Gubbins, K. E.; Liu, Y. C.; Moore, J. D.; Palmer, J. C. The role of molecular modeling in confined systems: Impact and prospects. Phys. Chem. Chem. Phys. 2011, 13, 58−85. (19) Talreja, M.; Kusaka, I.; Tomasko, D. L. Density functional approach for modeling CO2 pressurized polymer thin films in equilibrium. J. Chem. Phys. 2009, 130, 084902. (20) Müller, M.; MacDowell, L.; Yethiraj, A. Short chains at surfaces and interfaces: A quantitative comparison between density functional theories and Monte Carlo simulations. J. Chem. Phys. 2003, 118, 2929−2940. (21) Li, Z. D.; Firoozabadi, A. Interfacial tension of nonassociating pure substances and binary mixtures by density functional theory combined with Peng−Robinson equation of state. J. Chem. Phys. 2009, 130, 154108. (22) Shin, H. Y.; Wu, J. Z. Equation of state for the phase behavior of carbon dioxide polymer systems. Ind. Eng. Chem. Res. 2010, 49, 7678− 7684. (23) Mognetti, B. M.; Yelash, L.; Paul, W.; K., B.; M., M.; Macdowell, L. G. Efficient prediction of thermodynamic properties of quadrupolar fluids from simulation of a coarse-grained model: The case of carbon dioxide. J. Chem. Phys. 2008, 128, 104501. (24) Rosenfeld, Y. Free-energy model for the inhomogeneous hardsphere fluid mixture and density-functional theory of freezing. Phys. Rev. Lett. 1989, 63, 980−983. (25) Rosenfeld, Y. Free-energy model for the inhomogeneous hardbody fluidApplication of the Gauss−Bonnet theorem. Mol. Phys. 1995, 86, 637−647. (26) Tarazona, P.; Cuesta, J.; Martinez-Raton, Y. Density functional theory of hard particle systems. Lect. Notes Phys. 2008, 753, 247−341. (27) Ebner, C.; Saam, W.; Stroud, D. Density-functional theory of simple classical fluids. 1. Surfaces. Phys. Rev. A 1976, 14, 2264−2273. (28) Tang, Y. P.; Wu, J. Modeling inhomogeneous van der Waals fluids using an analytical direct correlation function. Phys. Rev. E 2004, 70, 011201. (29) Yu, Y. X.; Wu, J. Z. A fundamental-measure theory for inhomogeneous associating fluids. J. Chem. Phys. 2002, 116, 7094− 7103. (30) Tripathi, S.; Chapman, W. G. Microstructure of inhomogeneous polyatomic mixtures from a density functional formalism for atomic mixtures. J. Chem. Phys. 2005, 122, 094506. (31) Yu, Y. X.; Wu, J. Z. Density functional theory for inhomogeneous mixtures of polymeric fluids. J. Chem. Phys. 2002, 117, 2368−2376. (32) Fischer, J.; Methfessel, M. Born−Green−Yvon approach to the local densities of a fluid at interfaces. Phys. Rev. A 1980, 22, 2836− 2843. (33) Rathke, B. Bubble formation in gas-saturated fluid: The detection of homogeneous nucleation. Prog. Rep. Assoc. German Eng.,. Ser. 3 2003, 770. (34) Otake, K.; Kobayashi, M.; Ozaki, Y.; Yoda, S.; Takebayashi, Y.; Sugeta, T.; Nakazawa, N.; Sakai, H.; Abe, M. Surface activity of myristic acid in the poly(methyl methacrylate)/supercritical carbon dioxide system. Langmuir 2004, 20, 6182−6186. (35) Jaeger, P.; Eggers, R.; Baumgartl, H. Interfacial properties of high viscous liquids in a supercritical carbon dioxide atmosphere. J. Supercrit. Fluids 2002, 24, 203−217. (36) Chacon, E.; Fernandez, E.; Duque, D.; Delgado-Buscalioni, R.; Tarazona, P. Comparative study of the surface layer density of liquid surfaces. Phys. Rev. B 2009, 80, 195403.

the relative location of the coexistence curves and the spinodal. To study the phase behavior and interfacial properties at high pressures in polymer−CO2 mixtures, a DFT based on more sophisticated EOS-like perturbed-chain SAFT,41 which better describes the local packing at high densities, should be developed. Such an effort is under way.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The Dow Chemical Company is acknowledged for funding and for permission to publish the results. REFERENCES

(1) Han, J. H.; Han, C. D. Bubble nucleation in polymeric liquids. 1. Bubble nucleation in concentrated polymer-solutions. J. Polym. Sci., Part B: Polym. Phys. 1990, 28, 711−741. (2) Goel, S. K.; Beckman, E. J. Generation of microcellular polymeric foams using supercritical carbon-dioxide. 1. Effect of pressure and temperature on nucleation. Polym. Eng. Sci. 1994, 34, 1137−1147. (3) Stafford, C. M.; Russell, T. P.; McCarthy, T. J. Expansion of polystyrene using supercritical carbon dioxide: Effects of molecular weight, polydispersity, and low molecular weight components. Macromolecules 1999, 32, 7610−7616. (4) Binder, K.; Müller, M.; Virnau, P.; MacDowell, L. G. Advanced Computer Simulation Approaches for Soft Matter Sciences I; Advances in Polymer Science; Springer-Verlag: Berlin, 2005; Vol. 173, pp 1−104. (5) Nawaby, A.; Handa, Y. P.; Liao, X.; Yoshitaka, Y.; Tomohiro, M. Polymer-CO2 systems exhibiting retrograde behavior and formation of nanofoams. Polym. Int. 2007, 56, 67−73. (6) Yang, Y.; Liu, D. H.; Xie, Y. B.; Lee, L. J.; Tomasko, D. L. Lowtemperature fusion of polymeric nanostructures using carbon dioxide. Adv. Mater. 2007, 19, 251−254. (7) Tsutsumi, C.; Fukukawa, N.; Sakafuji, J.; Oro, K.; Hata, K.; Nakayama, Y.; Shiono, T. Impregnation of poly(L-lactideran-cyclic carbonate) copolymers with useful compounds with supercritical carbon dioxide. J. Appl. Polym. Sci. 2011, 121, 1431−1441. (8) Kasturirangan, A.; Koh, C. A.; Teja, A. S. Glass-transition temperatures in CO2 and polymer systems: Modeling and experiment. Ind. Eng. Chem. Res. 2011, 50, 158−162. (9) Duarte, A. R. C.; Martins, C.; Coimbra, P.; Gil, M. H. M.; de Sousa, H. C.; Duarte, C. M. M. Sorption and diffusion of dense carbon dioxide in a biocompatible polymer. J. Supercrit. Fluids 2006, 38, 392− 398. (10) Binder, K.; Mognetti, B. M.; Macdowell, L. G.; Oettel, M.; Paul, W.; Virnau, P.; Yelash, L. Towards the quantitative prediction of the phase behavior of polymer solutions by computer simulation. Macromol. Symp. 2009, 278, 1−9. (11) van Konynenburg, P. H.; Scott, R. L. Critical lines and phase equilibria in binary van der Waals mixtures. Philos. Trans. R. Soc. London, Ser. A 1980, 298, 495−540. (12) Virnau, P.; Müller, M.; MacDowell, L. G.; Binder, K. Phase separation kinetics in compressible polymer solutions: Computer simulation of the early stages. New J. Phys. 2004, 6, 7. (13) MacDowell, L. G.; Virnau, P.; Müller, M.; Binder, K. Critical lines and phase coexistence of polymer solutions: A quantitative comparison between Wertheim's thermodynamic perturbation theory and computer simulations. J. Chem. Phys. 2002, 117, 6360−6371. (14) Virnau, P.; Müller, M.; MacDowell, L. G.; Binder, K. Phase behavior of n-alkanes in supercritical solution: A Monte Carlo study. J. Chem. Phys. 2004, 121, 2169−2179. (15) Mognetti, B. M.; Virnau, P.; Yelash, L.; Paul, W.; Binder, K.; Müller, M.; MacDowell, L. G. Coarse-graining dipolar interactions in 3839



Article

(37) Katsov, K.; Weeks, J. D. Incorporating molecular scale structure into the van der Waals theory of the liquid−vapor interface. J. Phys. Chem. B 2002, 106, 8429−8436. (38) Müller, E. A.; Gubbins, K. E. Molecular-based equations of state for associating fluids. A review of SAFT and related approaches. Ind. Eng. Chem. Res. 2001, 40, 2193−2211. (39) Gao, W. H.; Butler, D.; Tomasko, D. L. High-pressure adsorption of CO2 on NaY zeolite and model prediction of adsorption isotherms. Langmuir 2004, 20, 8083−8089. (40) Humayun, R.; Tomasko, D. L. High-resolution adsorption isotherms of supercritical carbon dioxide on activated carbon. AIChE J. 2000, 46, 2065−2075. (41) Gross, J.; Sadowski, G. Perturbed-chain SAFT: An equation of state based on a perturbation theory for chain molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260.

3840


Density-Functional Theory for PolymerâCarbon Dioxide Mixtures

Recommend Documents

Density-Functional Theory for PolymerâCarbon Dioxide Mixtures