Article pubs.acs.org/JPCC
Density Fluctuations of Carbon Dioxide in Cylindrical Nanopore Yibing Dai,†,‡ Xiaofei Xu,*,† and Yicen Liu†,‡ †
Center for Soft Condensed Matter Physics and Interdisciplinary Research, Soochow University, Suzhou 215006, China College of Physics, Optoelectronics and Energy, Soochow University, Suzhou 215006, China
‡
ABSTRACT: It is important to understand the phase behavior of carbon dioxide in nanoporous structure, which could determine the pore shape and size distribution during a foaming process. In this work, we study the density fluctuations of carbon dioxide in a cylindrical nanopore by using classical density functional theory. It is found that the density fluctuations in solvophobic nanopore could exceed its bulk value by two orders of magnitude. The thickness of the layer with large density fluctuation is more than 3 nm to the wall. With decreasing the solvophobicity, the density fluctuation in nanopore could still be large with a decreased thickness of large fluctuation layer. The large density fluctuations happen only at the liquid side nearby the predrying line. These results are useful for the foaming process to manufacture porous material with a uniform pore shape and narrow size distribution.
I. INTRODUCTION The phase behavior of carbon dioxide (CO2) in porous structure is an issue of practical interest for manufacturing polymer foams. A typical foaming process starts with the dissolution of CO2 in polymer by pressurization, and bubble nucleation by quickly dropping the pressure to the ambient value. Then, the cell grows due to the diffusion of excess CO2 in the polymer and finally stabilizes with solidification of the cell walls by cooling.1,2 The dynamics of the cell growth and stabilization is very complicated because it is controlled by many variables, such as the polymer viscosity, CO2 concentration, and confinement of the porous structure.1,2 In contrast with bubble nucleation, relatively few works have studied about bubble growth or stabilization process for their microscopic details.1,2 It is well known that saturated liquid at the hard planar wall prefers a vapor−liquid coexistence for its relative low free energy.3−6 The solvophobic wall could enhance the depletion effect, leading to the formation of a vapor−liquid-like interface. At a vapor−liquid interface, liquid displays capillary wave fluctuations with long correlation length.7−9 As a result, solvophobic wall enhances density fluctuation. Liquid CO2 in porous structure should have similar behavior, which has a great effect on the bubble growth and stabilization of foaming process. Relevant studies are available in literature, although they are not about CO2. For hard-sphere fluids in contact with curved substrates, the surface tension shows a nonanalytic dependence on the curvature.10,11 By using molecular dynamics simulations, Patel et al.12 found that density fluctuations near hydrophobic surface are significantly larger than that in bulk and are quenched as increasing surface hydrophilicity. The local compressibility nearby the wall is large so that the hydrophobic interfaces are relatively soft.13,14 The density fluctuations of © 2016 American Chemical Society
water near solute display a nonmonotonic dependence on solute size.15 Simulations also showed that density fluctuations could affect dynamics of molecules at interfaces and thermodynamic properties at equilibrium.14 All of those results mean that density fluctuation of solvent molecules gives much useful information about the phase behavior of solvent or interface solvophobicity.16,17 Although foaming is a dynamic and nonequilibrium process, it is possible to use equilibrium thermodynamics to obtain some key information about the process.18 For the nucleation process with high-energy barrier (e.g., >10 kT), the system will experience a multitude of subcritical fluctuations before a critical nucleus is formed. With quantifying these fluctuations by means of order parameters, equilibrium thermodynamics could provide a suitable framework to describe the nucleation process. Some studies in polymer and CO2 mixture have been performed following this idea.18,19 For bubble growth or stabilization process, equilibrium thermodynamics also could provide information like the ensemble-averaged density distribution and fluctuation under the given ambient condition. This information describes the CO2 behavior near the wall of the nucleated cell. In this work, we aim to provide this kind of information, which can be used for searching the optimal ambient condition to manufacture porous material with a uniform pore shape and narrow size distribution. In a nucleated porous structure, the nanoconfinement and wall’s solvophobicity are the two dominant factors to determine the behavior of CO2. By considering these two factors, we study the density fluctuations of CO2 in a modeled cylindrical nanopore. We use classical density functional theory (DFT) in Received: March 9, 2016 Revised: April 15, 2016 Published: April 18, 2016 9520
DOI: 10.1021/acs.jpcc.6b02459 J. Phys. Chem. C 2016, 120, 9520−9526
Article
The Journal of Physical Chemistry C equilibrium thermodynamics to study the problem. The DFT could capture the necessary microscopic details of CO2 and accurately predict the phase behavior.20 In the following, we first introduce the molecular model and DFT for CO2 fluids and then give some analytical expressions to calculate the density profiles and density fluctuations in a cylindrical pore numerically. Finally, we show the results and draw some discussion.
The numerical methods to calculate these functional are given in next section. We consider a cylindrical nanopre of radius R immersed in a fluid of CO2. At a given bulk chemical potential μ, the equilibrium distribution of CO2 in the nanopore is obtained by minimizing the grand potential functional W=F−μ
II. MODEL AND THEORY The CO2 molecules are described by a molecular model developed in our previous work, which gives accurate description for the phase behavior in both sub- and supercritical regions.20 Each CO2 molecule is modeled as a dimer (N = 2) of two identical spheres, with the bonding potential given by exp[−βVB(r1 , r2)] =
4πσ
ρ (r ) =
(1)
λ (r ) = μ −
n3(r ) =
∫ dr 1 ρ(r1)Θ⎝ σ2
nV 2(r ) =
⎜
⎛ ⎜
∫
ρ (r ) =
⎞ − |r − r1|⎟ ⎠
1 =− 4
∫ dr dr1 (ρ(r ) − ρ(r1)) Θ(|r − r1| − σ )u(|r − r1|)
(10)
1 exp[βλ(r )] σ
∫0
σ
dr1 [exp(βλ(r + r1))
+ exp(βλ(r − r1))]
r1 2
σ − r12
(11)
with
(4)
where ρ(r) is the number density distribution of CO2 segments and Θ is the Heaviside step function. The nonlocal term is given by ex Fnonlocal
(9)
III. NUMERICAL DETAILS A cylindrical coordinate system (r,θ,z) is used for the numerical calculation. The point of r = 0 is placed at the axis of cylinder; see Figure 1. We suppose that the cylinder is long enough and the density distributions far away from the cylinder ends change only in r direction, that is, ρ(r) = ρ(r). In this simple onedimensional case, eq 7 becomes26
⎞ − |r − r1|⎟ ⎠
⎛σ ⎞ r − r1 dr 1 ρ(r1)δ ⎜ − |r − r1|⎟ ⎝2 ⎠ |r − r1|
(8)
Equation 10 accounts for the number density fluctuation per unit of energy perturbation, and thus χ has an unit of 1/((nm)3· kT). Here we present the data by using the dimensionless quantity χ/χb, where χb is the bulk value of χ under a given condition.
The local term includes the contribution from excluded-volume effect of hard spheres, correlation due to chain connectivity and dispersion effect.22 The details of these expressions can be found in our previous publication.20 All those expressions are functions of the weight density functionals of Rosenfeld, which are defined as23 ⎛
δF ex + V (r ) δρ(r )
⎛ ∂ρ(r ) ⎞ χ (r ) = ⎜ ⎟ ⎝ ∂μ ⎠T
(3)
∫ dr 1 ρ(r1)δ⎝ σ2
(7)
where N = ∫ dr ρ̂(r) is the total number of particles and ρ̂(r) = ∑Ni=1δ(r − ri) is the particle density operator. The bracket < > denotes a grand canonical ensemble average. Equation 9 correlates the local number density with the total number of particles in the system. Results have shown that χ could well quantify the density fluctuations nearby a surface.24,25 At a fixed temperature, χ can be calculated by24,25
ϵ is the effective energy parameter of interaction. The parameter values are ϵ/kB = 171.16 K and σ = 2.795 Å as given in our previous publication.20 We use a DFT to describe the thermodynamics and phase behavior of CO2. The DFT is constructed from the perturbedchain statistical associating fluid theory equation of state.21 The Helmholtz free-energy functional is expressed as a sum of an ideal-gas term and excess terms. The excess Helmholtz free energy is decomposed as a local contribution of equation-ofstate effect and a nonlocal contribution from the long-range dispersion interaction
n 2 (r ) =
1
χ (r ) = β⟨Nρ (̂ r ) − ⟨N ⟩⟨ρ ̂(r )⟩⟩
(2)
ex ex F ex = Flocal + Fnonlocal
∫|r−r |=σ dr1 exp[βλ(r1)]
The density fluctuation of the system is described by
⎪
⎪
1 exp[βλ(r )] 2πσ 2
with
where σ is the effective diameter of the sphere and δ is the Dirac delta function. β−1 = kT stands for the temperature multiplied by the Boltzmann constant. The excluded volume of CO2 molecules is represented by hard-core interactions. Energetic interactions are described by the attractive part of the Lennard-Jones potential. Thus, the pair interaction potential between segments is ⎧ ∞, rσ dr1 [ρ(r + r1) − ρ(r)] σr 4 1
1
Note that the numerical implementation of these expressions is efficient by computer, despite their complex expressions. There were several numerical methods to calculate the r-dependent density profiles in nanopore in literature.27,28 In contrast with those methods, our formula has the virtues of high numerical efficiency and mathematical beauty of rigorousness. We assume the nanopore’s wall is a uniform continuum with constant atomic number density np. The repulsive part of external potential V(r) is the hard sphere potential, and the attractive part is a sum of CO2 molecule and substrate atom pair interactions, each with the form of ⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ u wf (r ) = 4ϵwf ⎢⎜ wf ⎟ − ⎜ wf ⎟ ⎥ ⎝ |r | ⎠ ⎦ ⎣⎝ | r | ⎠
(14)
Thus, the attractive part is Vatt(r ) =
∫wall dr1uwf (|r − r1|)
Figure 3. Density profiles of CO2 in a nanopore with radius R = 5 nm and ϵwf ̅ = 0 kT at T = 295 K.
(15) 9522
DOI: 10.1021/acs.jpcc.6b02459 J. Phys. Chem. C 2016, 120, 9520−9526
Article
The Journal of Physical Chemistry C pressure to 30 MPa, the CO2 molecules are repelled away from the wall, giving a density profile with a wide depletion layer. The depletion width increases as further decreasing the ambient pressure P. We calculate the equilibrium vapor layer thickness at various pressure values and show the data in Figure 4. The equilibrium vapor layer thickness is defined by R
Leq = 2
∫0 (ρ(r ) − ρb )r dr (ρg − ρb )R
(17)
Figure 5. Phase diagram of CO2 in a hard nanopore ( ϵwf ̅ = 0 kT). The blue solid line is the bulk binodal curve, and the green dashed line is the bulk spinodal curve. The red and purple solid lines are the predrying line of CO2 in a nanopore with radius R = 5 and 2.5 nm, respectively. The solid circles are the critical point of the predrying state.
Figure 4. Equilibrium vapor layer thickness of CO2 in a hard nanopore (ϵwf ̅ = 0 kT) with R = 5 nm as a function of pressure at T = 295 K.
where ρb is the bulk density at the given pressure. ρg is the bulk vapor density at the same chemical potential of ρb. For comparison, in Figure 4, we also show the data from the sharpkink approximation. In this approximation, the density profile for a given vapor layer thickness l is taken to be ⎧ 0, r>R ⎪ ⎪ ρ(r ; l) = ⎨ ρg , R − l < r < R ⎪ ⎪ ρ , r 7.05 MPa and is vapor-like as P < 7.05 MPa. The whole phase diagram of the drying transition in T−ρ diagram is shown in Figure 5. The transition line is named as predrying line following the definition in literature.30 All predrying lines lie in the region of compressed liquid. With decreasing the nanopore’s radius, the line moves away the bulk binodal curve. The predrying line ends at a critical point (the black solid circles),31 where CO2 completely dries the wall of nanopore (the contact angle is 180°). Below the critical point, CO2 can only partially dry the wall (the contact angle is
