Langmuir 1994,10, 1837-1840
1837
Monte Carlo Simulations of the Laplace Pressure Dependence on the Curvature of the Convex Meniscus in Thin Unwetted Capillaries Elena N. Brodskaya and Elena M. Piotrovskaya' Department of Chemistry, St. Petersburg University, Petrodvoretz, St. Petersburg 198904, Russia Received March 5, 1993. In Final Form: February 21, 1994" The conditions of the two-phase liquid-vapor equilibrium of a Lennard-Jones fluid in thin unwetted slitlike and cylindrical capillaries of different widths were investigated at the temperature 90 K by the grand canonical Monte Carlo method. One of the main aims of this work was the estimation of the surface tension y for the curved meniscus (cylindricaland spherical). It is shown that the form of the dividing surface does not influence the dependence of y for a simple liquid on the curvature of this surface. Properties of the systems confined to small volumes differ significantly from those of the corresponding macroscopic phases. As a rule such systems are nonuniform everywhere. Noticeable changes appear also in the phase equilibria, and it is necessary to take this effect into account when discussing disperse systems in colloid chemistry. Phase transitions in thin pores are investigated both theoretically' and by computer simulations.lz Analytical equations were first derived for the case of rather wide pores. At the same time the investigations of fluid behavior in very narrow pores (their widths can be compared to the molecular sizes) are of great practical and theoretical interest. Reviews on computer simulations of such systems have previously appeared.l14y5 A principal goal of the present work is to discuss vaporliquid phase coexistence in pores of different forms and sizes with absolutely hard walls (unwetted pores). The convex meniscus is formed in the unwetted pores, and the pressure difference AP of the coexisting liquid and vapor has to satisfy at equilibrium the Laplace equation
AP = yc,
(1)
where y is the surface tension and c8 is the average curvature of the dividing surface, that is, the surface of tension. Strictly speaking, eq 1is correct only for spherical surfaces. If the surface has the values of the main curvature radii differing substantially, then the surface tension even in equilibrium is a tensor with two components. The generalized form of the Laplace equation was proposed a Abstract published in Advance ACS Abstracts, April 1, 1994.
(1) (a) Evans, R. J. Phys.: Condem. Matter 1990,2, 8989. (b) Ball, P. C.; Evans, R. Mol. Phys. 1988,63,159. (2) Peterson, B. K.; Gubbins, K. E. Mol. Phys. 1987,62, 215. (3) (a) Heffelfinger,G.S.;van Swol,F.; Gubbins, K. E. J. Chem. Phys. 1988,89, 5202. (b) Heffelfinger,G. S.; Tan,Z.; Gubbins, K. E.; Marini Bettolo Marconi, U.; van Swol, F. Znt. J. Thermophys. 1988,9,1051. (c) Heffelfinger,G.5.;van Swol, F.;Gubbins,K. E. Mol. Phys. 1987,61,1381. (d) Peterson, B. K.; Gubbins, K. E.; Heffelfinger, G.S.; Marini Bettolo Marconi, U.; van Swol, F. J. Chem. Phys. 1988,88,6487. (e) Schoen, M.; Diestler,D. J.; Cushman,J. H. J. Chem. Phys. 1987,87,5464. (0Schoen, M.; Cushman, J. H.; Diestler, D. J.; Rhykerd, C. L., Jr. J. Chem. Phys. 1988, 88, 1394. (9) Schoen, M.; Rhykerd, C. L., Jr.; Cushman, J. H.; Diestler,D. J. Mol. Phys. 1989,66,1171. (h) Walton, J. P. R. B.; Quirke, N. Mol. Simul. 1989,2, 361. (4) Brodskaya, E. N.; Piotrovskaya, E. M. Chemistry and Thermodynamics of Solutiom (in Russian);LGU Publishing House: Leningrad, 1986; Vol. 6, p 54. (5) Nicholson, D.; Parsonage, N. G.Computer Simulation and the Statistical Mechanics ofddsorption; Academic Press: New York, 1982; p 398.
e a r l i e ~ ~But J certainly for the estimation of surface tension we used approximation 1. In very narrow pores, the meniscus curvature becomes rather large and the surface tension depends on the curvature. It is necessary to stress that our systems are strongly nonuniform, and so a bulk liquid phase does not appear anywhere inside the system. Nevertheless, the Laplace equation is still valid according to Gibbs' method of the description of nonuniform systems. In this case instead of a real two-phase system with a nonuniform zone between phases, one considers the model system which consists of two uniform phases and a membrane with zero width having a certain tension. If a system is large enough, the uniform phases really exist at a long distance from the interfacial zone. But in the case of small systems with the overlapping surface layers the real uniform bulk phases may not exist. Then according to Gibbs' approach the model system contains hypothetical bulk phases with the values of temperature and chemical potential of the system under investigation. Namely, this approach was employed for the microdroplets of a Lennard-Jones fluid.gJ0 Here we intend also to follow the same method. The problem of the dependence of the surface tension on the surface curvature radius, first formulated by Gibbs, is still discussed in the l i t e r a t ~ r e . ~The t ~ dependence of y on the radius of a droplet was obtained by computer simulations of clusters of a simple liquid.lOJ1 One more way of solving this problem by investigating two-phase equilibrium in narrow pores is given in this paper. It also accounts for the influence of the form of the surface on the surface tension. That is why we considered both parallel-wide slits and cylinders, so that the form of the meniscus was either cylindrical or spherical.
Simulations The present work is devoted to Monte Carlo (MC) simulation of a liquid in thin slit and cylindrical capillaries with absolutely hard unwetted walls. Pair molecular (6) Krotov, V. V.; Rusanov, A. I.; Blinovsky, A. Colloid J . (in Russian) 1982,44,471. (7) Stuke, B. Chem.-Zng.-Tech.1961, 33, 173. (8) Gibbs,J. W. The collected works in two uolumes;Longmans,Green and Co.: New York, 1928. (9) Rueanov, A. I. Phasegleichgewichte und Grenzfldchenerscheinungen; Akademie-Verlag: Berlin, 1978; 4655. (10) Rusanov, A. I.; Brodskaya, E. N. J. Colloid Interface Sci. 1977, 62, 542; Colloid J. (in Russian) 1977, 39, 642. (11) Thompson, S. M.; Gubbins, K. E.; Walton, J. P. R. B.; Chantry, R. A. R.; Rowlinson, J. S. J. Chem. Phys. 1984,81, 530.
0743-7463/94/2410-1837$04.50/00 1994 American Chemical Society
Brodskaya and Piotrouskaya
1838 Langmuir, Vol. 10, No. 6, 1994
interactions are described by a truncated Lennard-Jones potential, @(r)(the cutoff distance is rc, at r > rc @(r)= 0): +(r) = 4c[(a/r)'* - (a/r)'l
(2)
where E and u are the parameters of the potential and r is the intermolecular distance. All calculations were held in the reduced units for which the parameters t and u were taken equal to unity. The slit pore is formed by two planar surfaces parallel to the xy plane. The walls were considered to be absolutely hard. The sizes of the basic MC cell in x and y directions were constant and equal to 1, = I, = 4.4or 8.8. Periodic boundary conditions (PBC) were used for the system in these directions. The slit spacings (pore widths) are H = 3.0, 6.0, and 12.0. For cylindrical pores the cylinder axis coincided with the x axis and PBC were used in this direction with the cell size of 8.8. We considered the pores of the following radii: R = 1.5, 3.0, and 6.0. Computer simulation was performed in the grand canonical ensemble (at the given volume V, temperature T , and chemical potential 1).The length of the MC chains for averaging was about lo6 configurations.
Two-Phase Equilibria The grand canonical ensemble describes open systems. That means that at the given value of chemical potential one can observe during the simulation only a one-phase state. In our case it is either a liquid-like or a vapor-like phase. The real coexistence of two phases can be obtained only for closed systems, which corresponds to a canonical NVT ensemble. But we use the grand canonical ensemble, and the condition of the two-phase coexistence (namely, the chemical potential pex) was gained by calculation of the grand thermodynamic potential Q separately for both phases. The main aim of the present work was to investigate the dependence of the two-phase liquid-vapor coexistence in thin unwetted pores on the pore size. The calculations were carried out for slit pores with H = 3.0, 6.0, and 12.0 and cylindrical pores with R = 1.5, 3.0, and 6.0 at the temperature 0.75. We have used the procedure proposed in ref 2 for the computer simulation of capillary condensation in adsorbing (wetted) pores. This method is connected with the determination of the dependence of the grand potential Q on the chemical potential p for gas and fluid branches of the adsorption isotherm. The intersection of two branches of Q is the point of phase equilibrium. The dependence Q ( p ) can be estimated by integrating Gibbs fundamental equation for Q. It is well known that the grand thermodynamic potential Q is a function of the volume V, temperature T ,chemical potentialp, and surface area A.18299As all our calculations were held at constant volume V and area A we can write the fundamental equation for R in the form dQ = -S d T - N dp
(3)
where S is the entropy and N is the number of particles. As an initial point for integration it is possible to choose a point on the supercritical adsorption isotherm (T = 1.2) at very low pressure, so that the gas may be considered to be ideal. Then Qid(p) = -kTN%). The supercritical adsorption isotherm is continuous and reversible, so one can integrate along this isotherm up to the pressure of the dense fluid according to the equation (aniap)T,V,A= -N
(4)
t
/
lm
-15
-14
__
- ---.- --@-iI/
-r
r'
-@
-12
-13
-10
-1:
'l
-9
-8
-?/LC
Figure 1. Adsorption isotherms for slit capillaries of different widths: H = 12.0 (I),H = 6.0 (II),H = 3.0 (III),at supercritical temperature T = 1.2 (la, IIa, IIIa) and the temperature of investigation T = 0.75 (Ib, IIb, IIIb).
N
500 r
-15
-14
-13
-12
-11
-10
-9
-8
-: -6
/"
Figure 2. Adsorption isotherms for cylindrical capillaries of different radii: R = 6.0 (I), R = 3.0 (II), R = 1.5 (III), at supercritical temperature T = 1.2 (Ia, IIa, IIIa) and the temperature of investigation T = 0.75 (Ib, IIb, IIIb).
The next step of this procedure is connected with the transition from the supercritical to a lower temperature at the given chemical potential. This dependence can be derived from eq 3 by substituting S = ( U - Q - ph9/ T and considering p constant. It yields
where U is the total energy of the system. After that according to eq 4 it is possible to calculate Q at any value of the chemical potential in the liquid phase at the temperature of investigation. The assumption of gasphase ideality was used to obtain the dependence of fl(w) for the gas. The intersection of the gas and liquid branches of Q ( p ) gives the value of the chemical potential pex corresponding to the vapor-liquid equilibrium in the pore at the given temperature. The isotherms p ( N ) at the temperatures T = 0.75 and T = 1.2 are given in Figures 1and 2 for slit pores (Figure 1)and cylindrical pores (Figure 2). The low-temperature isotherms exhibit two branches with metastable states and hysteresis. I t is necessary to point out that for our systems with unwetted walls the supercritical isotherms at T = 1.2 lie above the isotherms at T = 0.75. This result is different from that of the work2for the adsorbing pores where the supercritical isotherms lie below the isotherms at lower temperatures. We fulfilled the integration along the continuous isotherm up to the value of p corresponding to the dense fluid, and after that we lowered the temperature at the constant chemical potential, first to T = 1.0,and then to T = 0.75 (temperature of investigation).
Convex Meniscus in Thin Unwetted Capillaries
Langmuir, Vol. 10, No. 6,1994 1839 Table 1
slit capillary pel 10-2 0.34 -9.40 i 0.05 0.43 -9.26 f 0.05 0.52 -8.75 f 0.05 1.02
H
P
P*X
m
12.0 6.0 3.0
R OD
0.30 f 0.05 0.43 f 0.08 0.86 f 0.15
6.0 3.0 1.5
cylindrical capillary PBl10-2 0.34 -9.26 f 0.05 0.52 -8.70 f 0.05 1.10 -7.73 f 0.05 4.00 k
P
X
0.43 i 0.08 0.90 f 0.15 1.71 f 0.30
Table 2 slit capillary
cylindrical capillary
H
AP
c8
R,
Y
YKelv
12.0 6.0 3.0
0.30f 0.05 0.43 f 0.08 0.85f 0.15
0.2 0.5 2.0
5.0f 0.5 2.0* 0.5 0.5f 0.5
1.47 0.86 0.43
0.74 0.53 0.24
4
R
1
-40
-80
5
\
R 6.0 3.0 1.5
AP
Y
YKdv
CB
R,
0.43 f 0.08 0.89f 0.15 1.67 f 0.30
0.86 0.41 0.27
0.53 0.34 0.24
0.5 2.2 6.3
4.0 0.9 0.32
first of all, that the menisci formed in these pores are of the same average curvature and, besides, that the meniscus form does not influence the phase coexistence conditions. It means that either the meniscus form does not influence the dependence of surface tension on the curvature or the surface tension does not depend on the curvature in the investigated interval of the latter. As we shall see further on, the first statement is true.
Laplace Pressure The obtained data allow the Laplace pressure, that is, the pressure difference in coexisting phases separated by the curved meniscus (in our case either cylindrical or spherical) to be estimated. The vapor pressure P is obtained from the values of paxand considering the vapor to be ideal. The pressure P on the hypothetical liquid phase is calculated by integrating the equation of state for a homogeneous liquid phase dP = p dp, where p is the liquid density (not taking into account the compressibility of a liquid). So, it yields P = PI+ P(pex - PI). A special computer simulation was carried out for the bulk liquid phase a t T = 0.75 and p1, which was chosen in the range of rather high pressures to obtain reliable data for PI, according to the virial equation of state
j 8 7-
-10
-200
-9
t
\
- 600'
-40
-9
' P
-8
Figure 3. Dependence of the grand potential i2 on the chemical potential j t for two branches of the adsorption isotherms at T = 0.75: AB (liquid branch), CD (vapor branch); (a) slit pore with H = 6.0,pex= -9.26; (b) cylindrical pore withR = 6.0,pa = -9.26. According to the described procedure we have calculated the dependences Q(p) for all systems under investigation for both the liquid (Figure 3, curve AB) and gas (curve CD) branches of the adsorption isotherm a t T = 0.75. The intersections of these curves give the values of the chemical potential for liquid-vapor coexistence (Table 1). The increase of pexwith the decrease of the pore sue corresponds to the increase of the equilibrium vapor pressure P (Table l), but the values of P are not very precise because of the exponential dependence of this value on pep Earlier we made a computer simulation of the thick flat film of a Lennard-Jones liquid, surrounded by vapor, at the temperature T = 0.75.12 The value of the chemical potential of the vapor-liquid coexistence at the flat surface was pex = -9.58, and it correlates quite well with the results of the work.13 The values of per for a slit pore with H = 6.0 and a cylindrical pore with R = 6.0 (Table 1)coincide. It shows, (12)Brodskaya, E. N.; Piotrovskaya, E. M. Vestn. LCU (in Russian) 1989, 4,
43.
(13)Adams, D. J. Mol. Phys. 1976,29, 307.
where angular brackets mean ensemble averaging. We have chosen p 1 = -9.50,and the calculated value of PI was 0.228f 0.032. The average liquid density was p~ = 0.838 f 0.004. I t is necessary to stress that the accuracy of the surface tension determination depends strongly on the accuracy of the PIvalue. That is the reason why we have held extremely long computer simulation with the length of MC chains about 20 X lo6configurations. Earlier we made estimations of y on the basis of the less accurate calculations of PIand P.Thisfact explains the differences in the calculated values of this work with our previous re~u1ta.l~ The estimated values for P are given in Table 1. And again we see the same results for the capillaries of different form, but with the same average meniscus curvature: the slit with H = 6.0 and cylindrical capillary with R = 6.0. The Laplace pressure hp = P - PB (Table 2) is almost equal to the pressure in the liquid phase a because of the small values of P.
Surface Tension The data on the Laplace pressure AP allow the surface tension to be estimated according to the Laplace equation (14)Piotrovskaya, E. M.;Brodskaya, E. N. Colloid J. (in Russian)
1991,53, 673;1992,54, 51.
Brodskaya and Piotrouskaya
1840 Langmuir, Vol. 10, No. 6, 1994
Y
2
? \
1
0 0
1 I
1
2
3
4
5
cg
Figure 4. Dependence of the surface tension y on the curvature of the dividing surface cI at T = 0.75: (1)our results for slit capillaries; (2) our results for cylindrical capillaries;(3) results of ref 10 for droplets.
(1). Due to symmetry reasons we suppose that the meniscus of the slit capillary has a cylindrical form, and for the cylindrical capillary the meniscus is spherical. Then the curvature radius R , and the average curvature c, are connected by the following equations: c, = 1/R, for the slit capillary, c, = 2/R, for the cylindrical capillary. It is necessary to know the radius of the dividing surface, that is, the surface of tension, in order to use eq 1. The position of the surface of tension is not determined strictly from these simulations, but it is possible to estimate R, for a slit pore accordingto the earlier calculation datal2on computer simulations of a Lennard-Jones liquid in the slit unwetted pores. It is possible to find the position of the equimolecular surface due to the dependence of the liquid local density on the distance from the solid wall (density profile p ( z ) ) . This position happens to coincide almost with the position of the solid unwetted wall. So the radius of the equimolecular surface of the cylindrical meniscus in this pore can be estimated as H / 2 for absolutely hard walls. I t is known that the distance between the equimolecular surface and the surface of tension (value 6) is about 1.0." Then assuming 6 = 1.0, we obtain R , = (H- 2)/2 for slit pores and R , = R - 1 for cylindrical pores. Then the error in estimation of R, is f0.5. It means that for the pores with R, = 0.5 ( H = 3.0 or R = 1.5) the uncertainty in the definition of R, is up to 100%. It is also necessary to add that for narrow capillaries the meniscus form is not strictly spherical (for cylindrical capillaries) or cylindrical (for slit capillaries), but we assume that this supposition is true for the central part of the meniscus. The calculated dependences y(c,) are given in Figure 4 together with the results for droplets.1° The correlation of the results means that for the given values of the curvature the dependence of the surface tension of the simple liquid on the meniscus form (dividing surface) does not exist. In general the correlation of the dependences y(c,) for the spherical and cylindrical menisci is not obvious. That is why this result may be considered as a significant one, in spite of all the approximations used in the process. It is necessary to mention that it would be
useful to obtain the equimolecular surface curvature radius for the meniscus. The solvation of this problem needs special calculations for the closed two-phase system in a capillary. In the absence of such calculations (which are under study now) we have to use rather arbitrary assumptions. For the cylindrical capillary we decided to estimate not the surface tension but the dividing surface radius of the meniscus from the Laplace equation, using the same values of the surface tension at the given AP as for the slit capillary from the extrapolated dependence y ( h p ) . The results of the estimations are also given in Table 2 (the last column). As it is seen from the table, we again obtained the correlation of the values of the curvature radii for the menisci of different form for H = 6.0 and R = 6.0, and also for H = 3.0 and R = 3.0. It could be awaited according to all previous results. It is necessary to repeat that for the most narrow pores, R = 1.5 and H = 3.0, the data are obtained with substantial error. It is known that the vapor pressure above the curved meniscus is connected with the surface tension by the Kelvin equation k ~ p In(P/Pm@) " = yc,
(7)
where PJ is the vapor pressure at the flat surface. Rather often this equation is used in the literature for the equimolecular surface, though strictly it can be used only for the surface of tension. If the vapor pressure in the system is available, it is possible as well to use this equation for the calculation of the surface tension. The radius of the dividing surface in this case was estimated in the same way as before. The calculated results for y ~ are given in Table 2. As it is seen from the table we obtained satisfactory agreement of the values of the surface tension calculated by the two methods for both cylindrical and slit capillaries, and this agreement for thinner capillaries is even better. This fact is connected with the more precise estimation of Pflin these cases. These two methods of the estimation of y(cs) differ from each other only in the choice of a reference state with the pressure PIor P J . As it was mentioned above the error in defining P J values is rather high, so from our point of view, the first method of calculating y is more accurate. Comparing the method proposed in this work with the other one,lOJ1it is necessary to mention as the drawback of the present method the impossibility of precise determination of the surface of tension radius. The advantage of this method is in the direct calculations of the Laplace pressure without estimations of the vapor pressure, which is calculated with large errors. In general the proposed method of calculations of the surface tension of a liquid near the strongly curved surfaces may be considered to be perspective. I t is necessary to point out in conclusion that the most important result of this work is the absence of the influence of the form of the dividing surface on the dependence of the surface tension of a simple liquid on the curvature of this surface.
~
l
~
