Ind. Eng. Chem. Res. 1990,29,977-988 etry. Znd. Eng. Chem. Res. 1988,27,1936-1941. VDI W h e a t l a a . Bereehnungsbhtter ftir den WheBbergang, 4th ed.; VDI-Verlag: Diisseldorf, 1984. Vortmeyer, D.; Schuster, J. Evaluation of Steady Flow Profiles in Rectangular and Circular Packed Beds by a Variational Method.
977
Chem. Eng. Sei. 1983,38,1691-1699. Received for review September 13, 1989 Revised manuscript received February 12, 1990 Accepted February 22, 1990
Structure of Dilute Supercritical Solutions: Clustering of Solvent and Solute Molecules and the Thermodynamic Effects RongSong Wu,? Lloyd L. Lee,**+and Henry D. Cochrant School of Chemical Engineering and Materials Science, University of Oklahoma, Norman, Oklahoma 73019, and Chemical Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831
The objective of this paper is to exhibit, on the molecular level, the relationships between the microscopic structure and thermodynamic properties of dilute supercritical solutions by application of the integral equation theories for molecular distribution functions. T o solve the integral equations, we use Baxter’s Wiener-Hopf factorization of the Ornstein-Zernike equations and then apply this method to binary Lennard-Jones mixtures. A number of closure relations have been used: such as the Percus-Yevick (PY), the reference hypernetted chain (RHNC), the hybrid mean spherical approximation (HMSA), and the reference interaction-site (RISM) methods. We examine the microstructures of several important classes of supercritical mixtures, including the usual “attractive”-type and the less known “repulsive”-type solutions. The clustering of solvent molecules for solvent-solute structures in the attractive mixtures and, correspondingly, the solvent cavitation in the repulsive mixtures are clearly demonstrated. These are shown to be responsible for the large negative growth of the solute partial molar volumes in the attractive case and the positive growth in the repulsive case. Integral equations also afford us a unique opportunity to study the microstructures of solute-solute interactions. There is strong evidence of solute-solute aggregation in extremely dilute supercritical mixtures, a picture consistent with experimental fluorescence spectroscopy evidence. 1. Introduction It has been known for over a century (Hannay and Hogarth, 1879, 1880) that gases at slightly supercritical temperatures can exhibit, for a solute of low volatility, solubility that increases upon compression from that predicted for an ideal gas mixture with solvent power orders of magnitude higher than the ideal gas prediction. Other properties of supercritical solutions (Eckert et al., 1983, 1986; Johnston et al., 1988; Subramanian and McHugh, 1986; Collins et al., 1988) can also vary dramatically with changes of temperature or pressure near the critical point (CP) of the pure solvent. Recent results (Kim and Johnston, 1987a,b;Debenedetti, 1987; Cochran et al., l987,1988a,b; Pfund et al., 1988) have interpreted these observations in terms of solvent structure about a dissolved solute molecule: properties change dramatically near the CP as the correlation in density fluctuations becomes long-ranged and the size of the cluster of solvent molecules about each solute increases. The purpose of the present paper is to use the integral equation calculations to explore the relationship between structure and properties of dilute supercritical solutions. We briefly review the pertinent experimental observations and the previous theoretical results relating to the sol u b o l v e n t cluster in typical supercritical solutions. Then we extend this examination in two areas: first, we consider the relationship between structure and properties when the characteristics of the solvent and the solute are reversed (i.e., when solvent T,> solute T,and solvent V, > solute VJ. Second, we propose to consider the distribution University of Oklahoma. *Oak Ridge National Laboratory.
of solute molecules about a central solute molecule to understand the nature and effects of solute-solute aggregation for states near the solvent CP. To begin, we review the experimental observations that have suggested clustering of solvent molecules around a solute molecule in dilute supercritical solutions (solute partial molar volume (PMV) measurements, solubility measurements, solute absorption and fluorescence spectra, and others) as well as some observations that have suggested a high degree of solute-solute aggregation in dilute supercritical solutions. We shall also review briefly the prior theoretical interpretations of these observations. The experimental observations of the large negative growth of the solute PMV by Eckert and co-workers (Eckert et al., 1983, 1986) of dilute solutions of large organic molecules dissolved in supercritical C02 and CzH4 were seminal in the developing understanding of the relationship between properties and structure of supercritical solutions. The partial molar volume was determined from density measurements as a function of solute concentration for dilute solutions; the density was measured by the vibrating tube technique. These workers observed that the partial molar volume of solutes (naphthalene, tetrabromomethane, camphor) became negative and very large in magnitude (ca. 100 times the solvent’s bulk molar volume) near the solvent CP. They concluded “...that the solution process represents a disappearance of more than 100 mol of solvent/mole of solute. This suggests some sort of clustering process, perhaps akin to electrostriction about an ion in a protic solvent.” Kim and Johnston (1987a,b) and Debenedetti (1987) both presented molecular interpretations of the partial molar volume results to suggest that each solute is surrounded by a cluster of solvent molecules which grows to 0 1990 American Chemical Society
978 Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990
about 100 excess solvent molecules near the CP. The cluster size is obtained from the experimental partial molar volume data and the isothermal compressibility of the pure solvent by means of the solution theory of Kirkwood and Buff (1951) as follows: r = pkTKT - PO; (1) where 0; is the partial molar volume of solute 2 at infinite dilution. KT = -V-l dV/dQT is the isothermal compressibility. The cluster size (of the solvent molecules surrounding a solute molecule) is defined as
r
= plG 3 = p1 Jmdr 4nr2 [gt ( r ) - 11
(2)
where g12(r)is the pair correlation function (pcf) of a pair of molecules of species 1 (=solvent) and 2 (=solute). G12 is called the solute-solvent fluctuation integral. The isothermal compressibility could be obtained from the 11fluctuation integral, GI, lim plkTKT = 1 + plGll = P2+
Similarly, the PMV at infinite dilution is (4) P ~ U ; = 1 - ~ i ( G 1 2- GiJ Thus, from the KB fluctuation integrals Gij, we could calculate the thermodynamic quantities such as the compressibility and the PMV, and from the pair correlation functions gii(r),we could calculate the fluctuation integrals. This connection opens the door for a molecular interpretation of the bulk behavior, since we could use integral equations to obtain the pcfs. The large increase of the (negative) PMV was observed to parallel closely the increase in the slope of the solubility curve and the increase in the isothermal compressibility near the CP of the solvent; these phenomena are related by well-known thermodynamic expressions; for example, see Gitterman and Procaccia (1983). Cochran et al. (1987, 1988a) and Pfund et al. (1988) have also interpreted the dilute supercritical solubility results within the framework of Kirkwood-Buff theory as follows:
Again the solute-solvent fluctuation integral, G12, is the key quantity in the theoretical interpretation. Studies of chemical reactions in supercritical solvents (Collins et al., 1988) have also been suggestive of solvent clustering around a dissolved solute molecule near the solvent CP. Cochran and Lee (1989)and Cochran et al. (1988b, 1989) have earlier calculated pair correlation functions and fluctuation integrals for model supercritical solutions using integral equation theories. Consistent with experimental observations and earlier theoretical interpretations, they found growth of the solute-solvent fluctuation integral near the solvent CP. The growth results, in large part, from the increasing correlation length of density fluctuations near the CP. Other experiments (Yonker et al., 1986; Yonker and Smith, 1988; Kim and Johnston, 1987a,b Schneider et al., 1987; Kajimoto et al., 1988; Johnston et al., 1989) have explored the changes in the solute spectrum as the solvent is compressed. These works have shown a shift in the t absorbance or solvent-induced wavelength of maximum E changes in the UV fluorescence of a solute as the CP of the solvent is approached. The solvent-induced effects
were interpreted as reflecting an increased density of solvent molecules surrounding a solute molecule. Phenomenological interpretations of these results suggest that near the CP there is an increase in the local (short-range) solvent density around the solute. Such an increase was observed in the results of integral equation calculations (Cochran and Lee, 1989; Cochran et al., 1988b, 1989), as well as in a recent molecular dynamics simulation study (Petsche and Debenedetti, 1989). Other spectroscopicresults (Brennecke and Eckert, 1988, 1989) indicate that solute excimers (excited state dimers) are formed in very dilute solutions, well below saturation, as the CP is approached. The broad fluorescence band centered at ca. 480 nm indicates that about half of the pyrene dissolved in supercritical ethylene exists as excimers. These results have been interpreted to suggest that a high degree of solute-solute aggregation exists in dilute supercritical solutions near the solvent CP. The integral equation results of Cochran and Lee (1989) and Cochran et al. (1988b, 1989) are consistent with this interpretation. Like the solute-solvent correlation function, the solutesolute correlation function indicates both increased longrange correlation and increased short-range local density near the CP of the solvent. In the present paper, we shall extend the earlier theoretical interpretations of structure and properties of supercritical solutions using principally further results from integral equation calculations. Section 2 introduces the use of integral equation calculations for supercritical solutions. In Section 3 we describe the calculation method we have employed. Baxter's Wiener-Hopf factorization of the Ornstein-Zernike (OZ) equation has been used together with the Percus-Yevick (PY) closure. Section 5 presents the calculations for the "attractive" and "repulsive" mixtures by various closure relations (see, e.g., Lee, 1988). Considerable details are given for the latter case, due to its novelty and interesting features. Section 6 gives a full analysis of the solute-solute aggregation, a behavior sufficiently characterized by its microstructures. 2. Integral Equation Calculations Recently, Cochran and Lee (1989) and Cochran et al. (1988b, 1989) showed that pair correlation functions and fluctuation integrals for systems with interaction potentials representative of supercritical solutions can be calculated by solution of OZ equations:
where yij(r) = hij(r)- cij(r). Throughout this work except where noted otherwise, we use PY closure to study the supercritical solutions:
cij(r) = [l + y&)]
( [exp
(7)
When we applied the highly efficient Labik-Gillan numerical technique (Gillan, 1979; Labik et al., 1985),which combines Newton-Raphson (NR)and direct iteration (DI) methods, for solution of the OZ equation, the maximum value of the isothermal compressibility that could be approached was pkTKT = 15; beyond that the calculation algorithm diverges, owing to instability caused by truncation of hi@) in the Fourier transform space. The Labik-Gillan method requires the discretization of yij(r) and ci;(r)into a set of N values representing hi,(r;) and cij(ri), where ri = iAr and i = 1, ..., N . The LabikGillan method also requires the assumption hij(r)= 0 for r L R = NAr. A similar truncation is not always necessary
Ind. Eng. Chem. Res., Vol. 29, No. 6 , 1990 979 for cij(r),i.e., when using PY or MSA closures with finite-ranged potentials. From the comparison by Cummings and Monson (1988) between numerical and analytical solutions of OZ equations, it is clear that any numerical scheme that requires that hij.(r)= 0 for r > R is unable to approach the CP closely. Thls suggests that a method not dependent on r-space truncation of hij(r)should be employed in the vicinity of the CP. On the basis of factorization of the OZ equation for a single component, Baxter (1968) proposed a numerical algorithm for the calculation of h(r) under the sole assumption that ~ ( ris) finite in range. He showed that if c(r) = 0 for r > R (i.e., the direct correlation function is finite in range, as it is in the PY and MSA approximations for finite potentials) then the OZ equation can be factored into two equations:
rh(r) = -q'(r)
+ 2 a p i R d t q ( t ) ( r- t ) h(lr - ti)
rc(r) = -q'(r)
+ 2apJRdt
q'(t) q ( t
- r)
(8) (9)
where q ( r ) = 0 for r < 0 and r I R (10) and q'(r) denotes dqldr. Baxter suggested that the two equations (8) and ( 9 ) could be used, in conjunction with a closure relation, to iterate on the function q ( r ) , thus yielding a numerical method for solving the OZ equation where the solution of h(r)at r 2 R is not required. An original version of such an algorithm for the singlecomponent case was derived by Baxter (1967),employing a single integral equation which coupled eqs 8 and 9. It was solved by Watts (1968) with the PY closure using the NR method with grid 0 . 0 5 ~and with the potential cutoff varied from 3 . 5 to ~ 6 . 0 ~ To . approach the CP, a series of runs with various densities around different isotherms was required. Recently a variant of this method was introduced by Cummings and Monson (1990) which turns eqs 8 and 9 into a DI scheme for function q(r). They indicated that their DI method had the following difficulties: small grids, Ar = 0.005, were required, so a large number of points, N = 300, had to be used. When kT/t < 1.4, the DI scheme encountered convergenceproblems. They also coupled eqs 8 and 9 and solved with the NR technique. When they applied the N R algorithm, it was computationally intensive when N > 0 ( 1 0 2 ) ,and a good initial guess was necessary to ensure convergence. When we tried to extend this method to typical, binary supercritical solutions, very dissimilar molecules were encountered, and the above problems became worse. When the solute molecule has a large size (az >> all), we needed a large number of points, N , to cover the solute potential. Direct application of the NR method became harder because of the need to handle 4 times as many variables and a larger matrix to solve the simultaneous equations. For typical supercritical solutions, the solute molecules have large ea values compared to the solvent ell. When we reach the solvent CP, kT/tll is approximately 1.3, but we encounter very low reduced temperatures with respect to solute molecules (kT/t2, < 1 ) . These difficulties led us to seek a more powerful algorithm to solve the OZ equations for binary solutions near the solvent CP.
3. Method of Calculation Baxter reformulated the OZ equations of mixtures into two sets of equations given by
for Rij 2 r 2 Sij and RU
rhij(lrl) = -q$j(r)
2a&k&, k
dt qik(t)(r - t)kj(lr - tl)
(12) for r 1 Sit In these equations, Rij = (Ri + R j ) / 2 , Sij = (Ri - R j ) / 2 , and R , = min (Rki,Rkj - r]. The Ri are range parameters chosen such that cij(r)= 0 for r > Rij and qij(r) = 0 for r I Ri.. Our method for the numerical solution of the integral equation approximations in supercritical solutions is as follows. Step 1. By using the Labik-Gillan method to solve the OZ equations of mixtures, we can approach the CP to pkTK, = 15. The solutions are gij(r),ci.(r),and The direct correlation functions are used as input to solve functions qi.(r)in the next step. Step 2. 'he use the integral form of eq 11 Cij(r)=
1 "dr rcij(r) R..
r
=
The analytical solution of q lij(r)for hard spheres is used as an initial guess q$j(r)= air + bi for Sij < r < R , (14) ai = (1 - E3
+ 3Ri&J/(1 - &J2
bi = -3R:[2/2(1 -
(15) (16)
and 1 bi(r - Ri) qij(r)= -ai(r2 - R:) (19) 2 We discretize the functions qij(r),q:j(r),cij(r), and hij,(r) into qij(ru),q$j(ru),cij(ru),and kij(ru),where r, = uAr, m t h Ar = 0 . 0 4 ~ ~To~ simplify . the calculation, we use R1 = 2 . 5 6 and ~ ~ R2 ~ = 2 . 5 6 ~ ~(an - all). Using the trapezoidal rule, we write eq 13 in discrete form
+
+
Rm
Cij(rJ = Qij(r,) - 2aCPkAt C qki(tp)qkj(ru + t p ) { p k
tp=S,
u = 0, 1, ..., R i j / A r (20)
with {, = 1/2, when t, is at the upper or lower limit; otherwise {, = 1. Then we applied the NR method to solve eq 20. The Jacobian of eq 20 is dcij(rJ/aqmn(rJ = SimSjnSuu 2aAtCPk[tu;skmsinqkj(ru+ ru) + tuakmajnqki(rv - ru)] (21) k
The terms that survive inside the square brackets are as follows. The first term is nonzero when (22) Rkj - r, 2 T u 2 Ski, if Tu > Rkj - Rki or
Rki 1 ru 2 Ski, if ru
Rkj - Rki
(23)
980 Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990
and the second term is nonzero when Rk; ? Tu 2 ski + Tu, if Tu > Rk, .- R k i
We checked the quantity q2, defined by (24)
or Rki
+ P,
+ Tu,
2 ru 2 s k i
if
Tu
5 Rk; - Rki
(25)
ru
and = 1/2 if ru is at the upper or lower limit; otherwise, 1; = 1. To save time, we used Ar = 0 . 0 4 ~and ~ ~a short truncated potential to get an approximate solution of qij(r). As a measure of convergence of the NR cycles, we define the variable v1 as
The NR iteration is stopped when v1 is less than loa. Step 3. We interpolate q-(r,) from step 2 to a smaller ~ ~ extena ~ it to the desired cut-off pogrid Ar = 0 . 0 2 and tential -ycut = 5.12~11. Since qij(ru)is already a good approximation, direct iteration of eq 20 converges quickly to the final solution. We differentiate functions qij(r)to get q ’ij(r)It takes a long time to reach this point, but these three steps guarantee accurate solutions of hi,(r),cij(r),and q{j(r). This set of solutions is then used as a starting point to approach the CP. Step 4. Using hij(ru)and cij(ru)from step 1and q’&J from step 3, we have a set of solutions satisfying eqs 8 and 9. We proceed with direct iteration of eqs 7, 11, and 12. Suppose that c$(ru),h$(r,) and q’?(r,) are the nth iterates for the functions ci,(r),hij(r),and 9:j(r). The PY closure, eq 5, yields
r,hgUT(rU)= fi;(ru)[ruh$r)+ Tu - rUc$(rJ]- ru
(27)
where = exp[-ui;(r,)/kTl
fij(rJ
(28)
The ( n + 1)st approximations for hNEW(r,)is given by damping the functions h$(r,) and hi;$GT ( r J :
hflEW(r,)= ahgUT(r)+ (1 - a)hg(r,)
(29)
where (Y is a mixing parameter designed to ensure convergence. The damping is only required for the following region: ai; 5 r I1.6uij (30) Using eq 12 to solve for h?+l(ru)gives
ruh$+l(lrul)= -q’$(ru) + q$(t)(ru - t)hgEw((r,- t l ) (31)
2aEpkSsf?t k
The functions q’??l(r,) are u dated by adding the difference of rUh$+l(ru)and ruhij REW (Tu) with the damping factor 6: q’$+l(r,)
+ 9’3ru)
= @[rh;+l(r,)- rhflEW(r,)]
(32)
Then we can integrate 9’5+l(r)to get 95++’(r).The c$+’(rU) are solved by eq 11 r,c;+l(r,)
=
-9’$+’(rU) + 2*xPkS k
Rm
shi
dt 9ffi +‘(t) 9’tT1(ru+ t) (33)
The DI cycle was stopped when v2 was less than lo*. If eq 34 is not satisfied, we return to eq 27. The damping of functions q I-,(r)is a general requirement to guarantee convergence in the DI algorithm. The damping of functions h y w ( r )is made necessary, owing to the low reduced temperature with respect to solute molecules eZ2, where fij >> l. All integrations must be performed by Simpson’s rule to minimize the integration inaccuracies. This prevents the accumulation of errors that otherwise would cause the iteration algorithm used here to diverge when we approach the CP. (Note that we have ~ ~ of Ar = 0 . 0 0 5 ~ ~The ~. used grid Ar = 0 . 0 2 ~instead smaller grid was used by Cummings and Monson, who performed integrations with the trapezoidal rule.) Since the total correlation functions are long ranged near the CP, the integrals Gij in eqs 2 and 3 cannot be integrated directly. We calculate Gi; indirectly as follows. Define hijand t i .as Fourier transforms of functions hij and cij, respectively. Then (35) and
-s -
4T ti,(k) = k o dr r sin (kr)cij(r)
A t the limit k = 0, we have hi;(0) = 4 x S0m d rr2 hij(r)
(37)
zij(O) = 4 a X m d rr2 cij(r)
(38)
and
Since cij(r)is finite-ranged,eq 38 can be integrated directly; then using OZ equations in Fourier transform space, we have exact value of hij(0),which equals Gij. We have used the Baxter method to approach the critical point of Lennard-Jones (LJ) fluids and were able to determined their critical constants based on the PY approximation (Brey and Santos, 1985). The values are Tc* = 1.291 and pc* = 0.270. These were determined by plotting the inverse compressibility versus T * and p* (see Figure 1). (All units are reduced by the size parameter u and energy parameter e.) These numbers can be contrasted with the exact simulation values for LJ fluid: T* = 1.35, p* = 0.35. The differences are due to the approximate nature of the P Y approximation. The PY critical values determined here can be compared with other known estimates (Rowlinson and Swinton, 1982) for PY Tc*= 1.315 f 0.002; pc* = 0.280 f 0.002 (39) 4. Solute-Solvent Interactions in Attractive and Repulsive Mixtures Supercritical mixtures have been classified into attructiue, weakly attractive, and repulsive types (Debenedetti and Mohamed, 1989). Virtually all of the supercritical solutions in technological applications are cases of attractive clustering mixtures-e.g., decaffeination of coffee and other examples of supercritical fluid extraction, su-
Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990 981
O
'
0.ool
1 .za
O
V
8
HMSA
I
1.32
1.30
1
m
1.34
REDUCED TEMPERATURE T'
Figure 2. Solutesolvent pair correlation function at T* = 1.4and p* = 0.35. Xe in Ne (attractive mixture). HMSA closure, solute fraction = lo+. SJ
O.Ij
0.0 0.0
,
0.1
,
1
,\/&70 0.2
0.3
0.4
HMSA
0.5
REDUCED DENSITY p'
Figure 1. (a, top) Inverse of isothermal compressibility plotted against reduced temperature, for isochores p* = 0.27 (+), p* = 0.275 (b, bottom) Inverse of square root of iso(X), and p* = 0.28 (0). thermal compressibility plotted against reduced density, for isotherms T * = 1.291 (X) and T* = 1.308 (+).
percritical fluid chromatography, sorbent regeneration by supercritical fluids, condensation polymerization in supercritical fluids, etc. These mixtures are characterized by large negative solute PMV near the solvent critical point. By contrast, the "repulsive" mixtures exhibit large positive PMV near the solvent critical point. The question naturally arises whether repulsive mixtures exist in nature and whether they may be of importance. A thermodynamic analysis by Debenedetti and Mohamed (1989) showed that if dilute A dissolved in near-critical B forms an attractive system, then a dilute solution of B in supercritical A necessarily will form a repulsive mixture. For example, a dilute solution of xenon in supercritical neon is a typical attractive system;then a dilute solution of neon in supercritical xenon must be repulsive. There were also recent experiments (Biggerstaff and Wood, 1988a,b) on repulsive mixtures. The "existence" of repulsive mixtures is thus resolved. Since the Ne-Xe system can be well modeled as an LJ mixture, theoretical investigation of exemplary repulsive systems is also within reach. (The particular case of mixtures where the condition 0 < pD; < pkTKT is satisfied is called "weakly" attractive mixtures.) Attractive Mixtures. Petsche and Debenedetti (1989) carried out computer simulation for both attractive and
I/ I
1
,
I
r/all
Figure 3. Solute-solvent correlation function at T* = 1.4 and p* = 0.8. Xe in Ne (attractive mixture). HMSA closure, solute fraction = 104.
repulsive solutions, with parameters simulating the Ne-Xe system (Table I). We have solved the HMSA (Zerah and Hansen, 1986) closure for the attractive Ne-Xe case at the following conditions: T * = 1.4, p* = 0.35; T * = 1.4, p* = 0.8. Figures 2 and 3 show the solvent-solute correlation functions g12(r). One observes that, at T * = 1.4 and p*
982 Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990 Table I. Parameters for Neon and Xenon Used in the Simulationa 1
neon 1 I
size ratio, .j/oNe energy ratio, ci/bNs
xenon 1.435 7.04
"The MD used 864 particles with 863 solvent molecules and 1 solute molecule (Petache and Debenedetti, 1989).
0
5
10 L 011
15
20
Figure 5. Number of excess solvent molecules within a sphere of radius L around a solute, for the pure Lennard-Jones solvent (P), an attractive mixture (A; Xe in Ne), and a repulsive mixture (R; Ne in Xe), at T * = 1.4 and p* = 0.35. PY closures. Solute mole fraction = 10-9.
b
I
0.0
j
j
c
!
~
2
2
3
3
i
6
1
c
5
6
i
u
b
u
i
5
b
5
i
6
b
6
5
7
b
1
i
8
b
8
i
b
9
i
I 9
1
0 d
r I GII Figure 4. Pair correlation functions for solvent-solvent (gll), sol u h o l v e n t (gu),and solute-solute &) interactions at T * = 1.37 and p* = 0.4248. Pyrene-carbon dioxide attractive mixture, RHNC closure, solute mole fraction = 3 x io-'.
= 0.35 (near the solvent CP), g12exhibits a higher first peak and a longranged tail (with persistent tail values > 1). This long-rangedness disappears at a higher density (p* = 0.8). This behavior has been alluded to earlier as manifestation of clustering (Cochran et al., 1988b). We have also made calculations using the RHNC closure (Lado 1973) for the carbon dioxide-pyrene mixture studied by Brennecke and Eckert (1988). Solute pyrene is very dilute in C02 (mole fraction 3 X lo-'). The results for the three correlation functions, solvent-solvent, solvent-solute, and solutesolute, are shown in Figure 4. Again we clearly see the long-rangedness of the solvent-solute correlations. This enrichment translates into a large value for the fluctuation integral Glz. We note that G12is closely related to the excess coordination number, IVY2:
This number gives the excess of solvent molecules 1 radially surrounding the central solute molecule 2 up to radial distance of L. For random distributions, N",(L=m) = 0, meaning that there is no excess of solvent molecules around the center. The density of solvent molecules is given simply by the uniform density ( p l ) , and its radial number is given by (s/6)plL3, L being the coordination distance from the center. When N12or G12is greater than zero, there is positive correlation of 1 around 2; less than zero, there is negative correlation of 1 around 2. Figure 5 shows the Gij calculated from the PY closure for three
types of fluids: A = attractive mixture GI2,P = pure solvent Gll, and R = repulsive mixture G12 at T* = 1.35 and p* = 0.4 (all variables reduced according to solvent parameters). The curves represent the excess solvent neighbors at increasing coordination distances ( L )from the central solute molecule. Clearly, for case A, G12 becomes very large (.: llo), indicating substantial clustering of solvents around the solute molecule. This value could be compared with the pure solvent-solvent Cll E 8 and repulsive solvent-solute G12 E -6. We note that the major contribution to excess number comes from long-range contributions. Namely, large clusters did not derive from first neighbors alone. Upon reflection on geometry, this stands to reason since there is only so much room in the first neighbor shell for packing molecules. For clusters of size approaching 100, they must come from second, third, and higher neighborhoods. Repulsive Mixtures. Similar calculations are carried out for the repulsive mixtures using PY. This time, infinitely dilute neon is dissolved in supercritical xenon. The following conditions are used: T* = 1.4, p* = 0.35; T* = 2.0, p* = 0.35; T * = 1.4, p* = 0.80; T * = 1.34, p* = 0.27. The results are plotted in Figure 6. At high density, p* = 0.8, the pcf g,, is oscillatory. Near the CP of xenon, gI2(r) becomes suppressed in magnitude after the first peak. The . a second peak barely makes it above unity ( ~ 1 . 0 0 7 )At lower temperature (T* = 1.34),the second peak is less than 1 (see Figure 6b). Thenceforth, the correlations do not rise above 1, representing undercorrelation between the neon solute and the xenon solvent. This behavior is exactly the opposite of the attractive case discussed above. The solvent molecules instead of clustering around the solute molecule stay away (cavitate) from neon. This behavior is clearly seen in Glz. Figure 7 gives the GI2 for the four repulsive mixtures. The deficit in solvent molecules as the critical point of xenon (in PY approximation, T,* = 1.29, pc* = 0.27) is approached reaches a -12 for the lowest temperature studied. Again, the deficit is made out of both short-range and long-range contributions, with the major part coming from long-range solvent depletion. We have compared our PY calculation with the simulation data of Petsche and Debenedetti (1989) (Figure 8). While there is statistical scatter in the simulation data, the two results reinforce each other in the height of the first peak and
Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990 983
‘ 1
5
0 1.11
20
15
10
L 1611
,
Figure 7. Number of excess solvent molecules within a sphere of radius L around a solute, at T * = 1.34,p* = 0.27 (- - -); T * = 1.4, p* = 0.8 (--); T* = 1.4,p* = 0.35 (-); and T* = 2.0,p* = 0.35 (---) for dilute Ne in Xe (repulsive mixture). PY closure. Solute mole fraction = IO4.
J r’=l,
Figure 6. (a, top) Solute-solvent pair correlation functions gI2(r)at T* = 2.0,p* = 0.35 (---); T* = 1.4,p* = 0.35 (--); and T* = 1.4, p* = 0.8 (---) for dilute Ne in Xe (repulsive mixture). PY closure. Solute mole fraction = lo4. (b, bottom) Solute-solvent pair correlation functions gI2(r)at T * = 2.0,p* = 0.35 (-- -); T * = 1.4,p* = 0.35 (--); T* = 1.34,p* = 0.27 (-); and T* = 1.4,p* = 0.8 (---) for dilute Ne in Xe (repulsive mixture). PY closure. Solute mole fraction = lo*.
locations of the peaks and valleys. In repulsive mixtures, we expect the wavelengths of maximum absorbency or the solvent effects on fluorescence to shift in directions opposite to those in the attractive mixtures. It will be of interest for future spectroscopic experiments to confirm the repulsive supercritical behavior. In addition, one expects to observe the concomitant rapid positive increase of the PMV near the solvent CP and sharp drop of solubility with increasing pressure, characteristic signatures of repulsive mixtures (Biggerstaff and Wood, 1988a,b). By using integral equations, solved for model supercritical mixtures, we are able to show on the molecular level the clustering and cavitation of solvent molecules around the solute molecules depending on whether the mixture is attractive or repulsive. The story does not end here. Further interesting structural manifestations are found for the solute-solute correlations and are discussed next. 5. Solute-Solute Aggregation Not only do solvent molecules tend to cluster about solute molecules in attractive mixtures, we also find a high
0 : 0
I
I
I
2
1
4
I 6
rloll
Figure 8. Solute-solvent pair correlation function for dilute Ne in Xe (repulsive mixture) at T* = 1.4, p* = 0.35. (-) Molecular dynamic simulation result (Petsche and Debenedetti, 1989);(- - -) integral equation result.
.0.27
Figure 9. Solute-solute pair correlation functions gn(r) at T* = 1.415 = 1.096Tc*,and different densities. Ppene-carbon dioxide mixture. PY closure. Solute mole fraction = IO+.
degree of solute-solute aggregation near the CP of the solvent. We shall discuss the attractive case first and then
984 Ind. Eng. Chem. Res., Vol. 29, No. 6. 1990 Table 11. Lennard-Jones Parameters elk, K co*-co2 225.3 pyrene-pyrene 662.8 C02-pyrene 386.4
Solute Solute
A
6,
3.794 7.14 5.467
the repulsive case. Both cases show interesting solutesolute correlations. Attractive Mixtures. Figure 9 shows the correlation function g2,(r) for the LJ mixture simulating the carbon dioxide-naphthalene system, which is inaccessible by computer simulation and has been obtained only through integral equations (e.g., by PY). For the attractive mixture, the solutesolute pair distributions exhibit increased height of the first peak near the CP, indicating increased shortrange solute concentration about a solute molecule, consistent with interpretations of the excimer fluorescence spectra by Brennecke and Eckert (1988, 1989). This short-range solute-solute pair structure is similar to but larger in magnitude than the short-range solute-solvent structure presented earlier. Likewise, the long-range solutesolute pair structure increases near the CP, suggesting a solute-solute cluster commingled with the solvent-olute cluster. Because the solute concentration is very small, the solute-solute cluster is only a statistical association. Despite the low mole fraction, this statistical solutesolute cluster could be expected to exhibit macroscopic effects, for example, in solute-solute dimerization reactions, near the CP. In cases when the correlation length becomes large (Le., near the CP), the notion of dilute solutions being approximated as infinitely dilute must be employed with greater caution since despite the rarefaction of solutes, detectable aggregates do occur. In addition to showing increased height of the first maximum near the CP, the short-range solute-solute pair distributions for attractive mixtures (Figure 9) show a much shallower minimum between the first and second maxima; this is similar to the behavior of the solutesolvent pair distributions near the CP but much exaggerated. As mentioned earlier, the fluorescence spectroscopic experiments of Brennecke and Eckert (1988, 1989) indicated significant excimer formation for very dilute pyrene dissolved in supercritical ethylene and carbon dioxide. We shall report now some additional integral equation calculations for the carbon dioxide-pyrene system. The Lennard-Jones parameters are given in Table 11. The pressure used is 82.5 bar (or 1197.1 psia), concentration y yr = 3.0 X lo-', and density (pure CO,) = 0.8059 lb-mol/&. The RHNC equations (Lado, 1973) were solved for this case. The three pair correlation functions gll, g,, and g, are shown in Figure 4. The solute-solute correlation function shows a curious double peak with a shallow minimum in between. We shall infer that this feature (first minimum greater than unity) characterizes the affinity of pyrene molecules for one another. To determine the origin of this shallow minimum, we need to examine the geometry of packed molecules. For this case the radial packing geometry is shown in Figure 10. The first peak of gZ2(r)represents the direct contact between a pair of pyrene molecules (direct pyrene-pyrene correlation). The second peak represents the geometry of a pair of pyrene molecules sandwiching a CO, molecule in between (pyrene-pyrene correlation mediated via a solvent molecule: i.e., pyrene-CO,-pyrene). There is a perceptible shoulder after the second peak. It is the normal location for triplet pyrene correlations (concatenations such as pyrene-pyrene-pyrene). Far from the CP, the first minimum in correlation functions usually falls below 1.0 because of the volume excluded by an intervening solvent
+' Solute
Solute
Solute Solute Solute
'
Shouldemrd peak) L.*2a,,I
Figure 10. Radial packing geometry of g&)
peaks.
Table 111. Lennard-Jones Parameters
elk. K monomer-monomer, mm dimer site-dimer site, ss monomer-dimer site, ms
225.3 225.3 225.3
A
6,
3.794 5.691 4.743
molecule. This does not happen here. In this case at r = 2.7g11, where normally a CO, molecule would intervene between the pyrene molecules, the usual gap with g,, < 1.0 fails to occur. Certainly, this is not due to scarcity of solvent molecules (as pyrene is infinitely dilute in CO,). This could be explained by the aggregation of pyrene molecules at this state near the CP, in such a way that causes exclusion of CO, and simultaneously relative enrichment of pyrene molecules. This behavior is consistent with and strengthens the picture of excimer formation. To determine how dilute is dilute, we have also calculated a case at a smaller concentration of pyrene (using RHNC at the same temperature and pressure), yp = 2.35 x lo-". Very similar gy(r) are obtained. In &ct, the differences in the numerical values of gij(r)for the two concentrations are in the fourth decimal place. We also examined the temperature effect on radial structure. Figures 4 and 11 show the three radial distribution functions gll, g12,and g , at three temperatures-T* = 1.37, 1.50, and 2.00. Besides the ordering of peak heights, the noticeable feature is the deepening first minimum of g, as the temperature is increased. Starting at g p = 2.0 (enrichment), the value at the minimum finally drops to g,, = 1.0 at the highest temperature (T* = 2). Thus, the aggregation (detected by lack of exclusion of CO, molecules) of pyrene molecules is reduced at higher temperatures when conditions are farther removed from the CP. We find that these solute-solute correlations show intriguing similarities to the correlation functions observed for diatomics (while we want to make the distinction here that the solute-solute aggregation observed with integral equation theories does not represent bonding). To contrast the structural idiosyncrasies with the diatomics, we carried out RISM calculations for dimer-monomer mixtures. For the monomers (CO,) and dimers, the LJ parameters are given in Table 111. The bond length between the sites of the dimer is taken to be 1 = 0.66~- = 2.5 A. Figure 12 shows the site-site correlation functions g-, g,,, and g , (site s refers to an atom in the dimer). The first peak of g,, occurs at the collision diameter between
Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990 985 2
or------
- _ _ ~
3.5.
3.0.
Figure 12. Site-site pair distribution functions of monomer-mo(,), and dimer-dimer (gu)for dinomer (g-), monomer-dimer g mer-monomer mixtures. RISM theory. PY closures.
b
I
1 g12
1
j j j
0.0 0
1
I
Z
2
3
3
U
V
5
S
6
6
1
1
8
0
9
9
i
b
i
b
i
b
i
b
i
b
i
b
i
b
i
b
i
b
i
1
?
r 1011 Figure 11. (a, top) Pair correlation functions for solvent-solvent and solute-solute (gZ2)interactions at T* (sll),eolute-solvent = 1.50 and p* = 0.4248. Pyrenecarbon dioxide attractive mixture, RHNC closure. Solute mole fraction = 3 X lo-'. (b, bottom) Pair solute-solvent (gll), correlation functions for solvent-solvent (gI1), and solute-solute (gZ2)interactions at T* = 2.00 and p* = 0.4248. Pyrene-carbon dioxide attractive mixture, RHNC closure. Solute mole fraction = 3 x io-'.
the sites in the two dimers (=1.86u-). The second peak is due to the correlation of one site of the first dimer with the second (farther) site of the other dimer, r = 2.34~". We shall designate one site in the first dimer as the "host", the nearby site in the second dimer as the "guest", and the farther site in the second dimer as the "companion" (to the guest). In this language, the second peak is due to the correlation between the host (in molecule 1) and the companion (in molecule 2). The increment in distance is, of
course, the bond length 0 . 6 6 ~ ~ In.fact, the second peak occurs at a shorter distance became all possible nonaligned angular orientations of the second dimer tend to reduce the distance from 2.34~- to about 2.20- The third peak comes from a pair of sites mediated by a monomer molecule and occurs at =2.8um. An interesting point is the similarity in the qualitative features of g, with those of the low-temperature pyrene pyrene g., In both cases, the first minimum is shallow and much greater than unity. This phenomenon has been explained geometrically earlier. If we focus our attention on the correlation between the site called the host in the first dimer and the companion in the second dimer, then it is clear that in between these two there is the guest site of the second dimer that intervenes and excludes a solvent molecule from entering. This explains the probabilistic enrichment that forms the shallow first minimum; i.e., probabilistically, it is highly likely to have a neighborfrom a site in the second dimer-localized here, with a probability density above and beyond the average density of dimers. The same was observed for the low-temperature LJ mixtures modeling COz-pyrene. If pyrene and pyrene tend to aggregate, then there is a increased likelihood of having a third pyrene molecule sandwiched between a host pyrene molecule and a companion pyrene molecule. This third pyrene molecule is tied to the companion pyrene molecule because of aggregation. There are, of course, differences between the structures of the monomer-dimer case and the COz-pyrene case. The three peaks in the former are due to host-guest, host-companion, and hostsolvent-guest correlations, whereas the peaks for COzpyrene are due to correlations of host-guest, host-solvent-guest, and host-guest-second guest. Repulsive Mixtures. It is of interest to investigate the solute-solute interactions in the repulsive mixtures discussed above. We have seen the cavitation of the solvent molecules away from the solute molecule in the Ne-Xe system. How do the solute molecules (in this case neon) behave toward each other? We have calculated the gzz using PY for the four state conditions listed in section 4. The behavior is shown in Figure 13. As the density of solvent (xenon) is lowered from 0.8 to 0.35, the oscillations in Ne-Ne correlations are damped. If we calculate the excess number (eq 40; see Figure 141, we note that there is enrichment of excess numbers when the critical condition of xenon is approached. A t high temperatures, T *
986 Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990
0 0
IJ, 0
1
I ,
,
2
,
,
,
41
3
4
5
6
r'Q11
-
2
-
10
15
20
L 1611
Figure 14. Number of excess solute molecules within a sphere of radius L around a solute, at T* = 2.0, p* = 0.35 (---); T* = 1.4, p* = 0.35 (--); T* = 1.34, p* = 0.27 (-); and T* 1.4,p* 0.8 (---) for dilute Ne in Xe (repulsive mixture). PY closure. Solute mole fraction = 10".
, . ,
I
: i
09 0
2
4
r I QII
Figure 13. (a, top) Solute-solute pair correlation functions gaz(r)at T* = 2.0, p* 0.35 (---); T* = 1.4, p* = 0.35 (--); and T* 1.4, p* = 0.8 (- -) for dilute Ne in Xe (repulsive mixture). PY closure. Solute mole fraction = IO!. (b,bottom) Solute-eolute pair correlation functions g&) at T*= 2.0, p* = 0.35 (---), T* = 1.4, p* = 0.35 (-4; T* = 1.34, p* = 0.27(-1; and T* = 1.4, p* = 0.8 (---) for dilute Ne in Xe (repulsive mixture). PY closure. Solute mole fraction = 10-9.
-
= 2.0, there is deficit in neon-neon correlations (G22< 0). At high densities, Gz2rises above zero. Namely, neon molecules tend to "separate" out of the mixture. Near the
critical point of xenon, the separation increases rapidly. Since neon molecules are rejected by the xenon solvent, evidenced by the solvent cavitation discussed earlier, Ne tends to form aggregates too. By analogy with nonpolar molecules dispersed in an aqueous medium, this aggregation is caused by xenon-phobic attraction (cf. sic hydrophobic attraction). If we compare the solute-solute aggregation in repulsive mixtures with that in the attractive case, our calculations show that the former is several orders of magnitude smaller than the latter. The excess GZzfor the attractive case is on the order of lo3 near CP. 6. Concluding Remarks
The algorithm we have used for solving the OZ equations of mixtures near the CP has proved to be robust and efficient compared with other methods we have tried. The results by the Wienel-Hopf factorization have been shown to agree closely with results from the Labik-Gillan method,
but the new method can approach the CP much more closely. The CP for dilute LJ mixtures can be located quite precisely. This algorithm is applicable generally to mixtures, the components of which may be modeled as interacting with a pair potential finite in range. Using integral equation theories (PY, HMSA, RHNC, and RISM),we have demonstrated here the detailed microscopic behavior of the interesting repulsive supercritical mixtures proposed elsewhile by Debenedetti and Mohamed (1989), as well as the usual attractive mixtures. The results of the correlation functions have proved helpful in gaining understanding of the distribution of molecules surrounding a central solute molecule in dilute supercritical solutions. These results are consistent with observed macroscopic properties which are sensitive to the short-range structure-particularly spectroscopic results. The longrange solute-solvent structure for attractive mixtures near the CP from integral equation calculations is consistent with the thermodynamic properties which depend on the solute-solvent fluctuation integral G12-e.g., the solute partial molar volume and solubility. These long-range structures are as for now inaccessible by computer simulation techniques. The solute-solute distributions for attractive mixtures from integral equation calculations are also consistent with the observation of existence of a high proportion of solute excimers near the CP in very dilute supercritical solutions. To date there are no solute-solute distributions near CP from molecular simulation calculations (which must be performed at infinite dilution), nor are there results with which to compare the long-range solute-solute pair distributions. Judging by the integral equation results, one might expect dramatic effects on other applications such as chemical reactions in supercritical media. For repulsive mixtures near the CP, the integral equation results for solute-solvent pair distributions are in agreement with the simulation results at short range. We have clearly shown the solvent cavitation around the solute molecules. There is recent experimental work (Biggerstaff and Wood, 1988a,b) on repulsive mixtures. It would be interesting for future work to deduce the quantitative information from the molecular approach. The postulated solvent shift in spectroscopic studies of repulsive mixtures near the CP would be interesting to seek experimentally,
Ind. Eng. Chem. Res., Vol. 29, No. 6, 1990 987
as would possible effects on solute-solute bimolecular reactions. Acknowledgment This work was supported in part by the Chemical Sciences Division of the Office of Basic Energy Sciences of the U S . Department of Energy under Contract DEAC05-84R21400 with Martin Marietta Energy Systems, Inc., and subcontract with the University of Oklahoma. We thank the courtesy of Pablo G. Debenedetti for providing the preprint on repulsive mixture simulations. We thank David M. Pfund for making the HMSA calculations. We also acknowledge the use of computer time from a grant by the San Diego Supercomputer Center. Nomenclature a = hard sphere constant in eq 15 b = hard sphere constant in eq 16 c = direct correlation function f = Meyer function g = pair distribution function h = total correlation function k = Boltzmann constant k = wave number in Fourier space p u = vapor pressure q = Baxter’s function in eq 11 and 1 2 r, s, t = spatial separation u = potential energy u = discretization index u = volume y = mole fraction C = direct correlation function integral E = enhancement factor G = Kirkwood fluctuation integral KT = isothermal compressibility N = number of discrete points P = pressure R, S = range parameter T = absolute temperature V = volume 2 = compressibility factor Greek Symbols a = damping parameter in eq 29 p = damping parameter is eq 32 y = indirect correlation function y = cut-off distance 6 = Kronecker delta c = molecular energy parameter { = parameter in eq 21 q1 = tolerance in eq 26 q2 = tolerance in eq 34 p = number density u = molecular size parameter t2 = parameter in eq 17 t3 = parameter in eq 18 r = cluster size Subscripts c = critical point property cut = cutoff i = discretization index i, j , k , 1, m = component index m = minimum m = monomer site max = maximum p = index in eq 20 s = dimer site u = discrete u = index in eq 21 1 = solvent index
2 = solute index
Superscripts n = iteration counter n + 1 = iteration counter NEW = iteration counter OUT = iteration counter s = solid phase * = reduced quantity m = at infinite dilution ’ = first derivative = Fourier transform - = partial molar quantity A
Literature Cited Baxter, R. J. Method of Solution of the Percus-Yevick, Hypernetted-Chain, or Similar Equations. Phys. Reu. 1967,154, 170. Baxter, R. J. Ornstein-Zernike Relation for a Disordered Fluid. Aust. J . Phys. 1968,21,563. Baxter, R. J. Ornstein-Zernike Relation and Percus-Yevick Approximation for Fluid Mixtures. J. Chem. Phys. 1970,52,4559. Biggerstaff, D. R.; Wood, R. H. Apparent molar volumes of aqueous argon, ethylene, and xenon from 300 to 716 K. J. Phys. Chem. 1988a,92, 1988. Biggerstaff, D. R.; Wood, R. H. Apparent molar heat capacities of aqueous argon, ethylene, and xenon at temperatures up to 720 K and pressures to 33 MPa. J. Phys. Chem. 198813,92,1994. Brennecke, J. F.; Eckert, C. A. Molecular Interactions from Fluorescence Spectroscopy. Proc. Znt. Symp. Supercritical Fluids 1988,263. Brennecke, J. F.; Eckert, C. A. Fluorescence Spectroscopy Studies of Intermolecular Interactions in Supercritical Fluids. In Supercritical Fluid Science and Technology;Johnston, K. p., Penninger, J. M. L., Eds.; ACS Symposium Series 406;American Chemical Society: Washington, DC, 1989;Chapter 2. Brey, J. J.; Santos,A. On the Critical Behavior of the Percus-Yevick Equation for Nontruncated Potentials. J. Chem. Phys. 1985,82, 4312. Cochran, H. D.; Lee, L. L. Solvation Structure in Supercritical Fluid Mixtures based on Molecular Distribution Functions. In Supercritical Fluid Science and Technology; Johnston, K. P., Penninger, J. M. L., Eds.; ACS Symposium Series 406;American Chemical Society: Washington, DC, 1989;Chapter 3,pp 27-38. Cochran, H. D.; Lee, L. L.; Pfund, D. M. Application of the Kirkwood-Buff Theory of Solutions to Dilute Supercritical Mixtures. Fluid Phase Epuilib. 1987,34,219. Cochran, H. D.; Fund, D. M.; Lee, L. L. Theoretical Models of Thermodynamic Properties of Supercritical Solutions. Sep. Sci. Technol. 1988a,23,2031. Cochran, H.D.; Lee, L. L.; Pfund, D. M. Study of Fluctuations in Supercritical Solutions by an Integral Equation Method. Proc. Int. Symp. Supercritical Fluids 1988b,245. Cochran, H. D.; Lee, L. L.; Pfund, D. M. Structure and Properties of Supercritical Fluid Mixtures from Kirkwood-Buff Fluctuation Theory and Integral Equation Methods. Fluctuation Theory of Mixtures; Matteoli, E., Ed.; Advances in Thermodynamics Series; Taylor and Francis: New York, 1989,submitted. Collins, N. A.; Debenedetti, P. G.; Sundaresan, S. Disproportionation of Toluene over SZM-5 under Near-Critical Conditions. AZChE J. 1988,B34, 1121. Cummings, P. T.;Monson, P. A. A New Method for the Numerical Solution of Integral equation approximations. Int. J. Thermophys. 1990,11, 97. Debenedetti, P. G. Clustering in Dilute, Binary Supercritical Mixtures: A Fluctuation Analysis. Chem. Eng. Sci. 1987,42,2203. Debenedetti, P. G.;Mohamed, R. S. Attractive, Weakly Attractive, and Repulsive Near-Critical Systems. J. Chem. Phys. 1989,90, 4528.
Eckert, C. A.; Ziger, D. H.; Johnston, K. P.; Ellison, T. K. The Use of Partial Molar Volume Data to Evaluate Equations of State for Supercritical Fluid Mixtures. Fluid Phase Equilib. 1983,14,167. Eckert, C. A.; Ziger, D. H.; Johnston, K. P.; Kim, S. Solute Partial Molar Volumes in Supercritical Fluids. J. Phys. Chem. 1986,90, 2798. Gillan, M. J. A New Method of Solving the Liquid Structure Integral Equations. Mol. Phys. 1979,38,1781. Gitterman, M.; Procaccia, I. Quantitative Theory of Solubility in Supercritical Fluids. J . Chem. Phys. 1983,78, 2648.
I n d . Eng. Chem. Res. 1990,29, 988-994
988
Hannay, J. B.; Hcgarth, J. On the Solubility of Solids in Gases. Proc. R. SOC. London 1879,A29, 324. Hannay, J. B.; Hogarth, J. On the Solubility of Solids in Gases. Proc. R. SOC. London 1880,A30,484. Johnston, K. P.; Flarsheim, W.; Hmjez, B.; Mehta, F.; Fox, M.; Bard, A . Solvent Effect on Chemical Reactions at Supercritical Fluid Conditions. Proc. Znt. Symp. Supercritical Fluids 1988, 907. Johnston, K. P.; Kim, S.; Combs, J. SpectroscopicDetermination of Solvent Strength and Structure in Supercritical Fluid Mixtures, k Review. In Supercritical Fluid Science and Technology; Johnston, K. P., Penninger, M. L., Eds.; ACS Symposium Series 406; American Chemical Society: Washington, DC, 1989; Chapter 5. Kajimoto, 0.; Futakami, M.; Kobayashi, T.; Yamasaki, K. Chargetransfer-State Formation in Supercritical Fluid: (N,N-Dimethy1amino)benzonitrile in CF,H. J . Phys. Chem. 1988,92, 1347-1352. Kim, S.; Johnston, K. P. Molecular Interactions in Dilute Supercritical Fluid Solutions. Znd. Eng. Chem. Res. 1987a, 26, 1206-1213. Kim, S.; Johnston, K. P. Clustering in Supercritical Fluid Mixtures. AZChE J. 1987b,33, 1603-1611. Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Solutions. I. J. Chem. Phys. 1951,19,774. Labik, S.; Malijevsky, A.; Vonka, P. A Rapidly Convergent Method of Solving the OZ Equation. Mol. Phys. 1985,56,709. Lado, F. Perturbation Correction for the Free Energy and Structure of Simple Fluids. Phys. Reu. A. 1973,87,2548. Lee, L. L. Molecular Thermodynamics for Nonideal Fluids; Butterworths: Boston, 1988. Nicolas, J. J.; Gubbins, K. E.; Streett, W. B.; Tildesley, D. J. Equation of State for the Lennard-Jones Fluid. Mol. Phys. 1979,37, 1429.
Petsche, I. B.; Debenedetti, P. G. Solute-Solvent Interactions in Infinitely Dilute Supercritical Mixtures: A Molecular Dynamics Investigation. J. Chem. Phys. 1989,91,7075. Pfund, D. M.; Lee, L. L.; Cochran, H. D. Application of the Kirkwood-Buff Theory of Solutions to Dilute Supercritical Mixtures, 11. The Excluded Volume and Local CompositionModels. Fluid Phase Equilib. 1988,39,161. Rowlinson, J. S.; Swinton, F. L. Liquids and Liquid Mixtures; Butterworths: Boston, 1982; p 254. Schneider, G. M.; Ellert, J.; Haarhaus, U.; Holscher, I. F.; Katzenski-Ohling,G.; Koppner, A.; Kulka, J.; Nickel, D.; Rubesamen, J.; Wilsch, A. Thermodynamic,Spectroscopic and Phase Equilibrium Investigations on Polar Fluid Mixtures at High Pressures. Pure Appl. Chem. 1987,59,1115-1126. Subramanian, B.; McHugh, M. A. Reactions in Supercritical Fluids-A Review. Znd. Eng. Chem. Process Des. Deu. 1986,255, 1.
Watts, R. 0. Percus-Yevick Equation Applied to a Lennard-Jones Fluid. J. Chem. Phys. 1968,48, 50. Yonker, C. R.; Smith, R. D. SolvatochromicBehavior of Binary Supercritical Fluids: The Carbon Dioxide/2-Propanol System. J. Phys. Chem., 1988,92,2374-2378. Yonker, C. R.; Frye, S. L.; Kalkwarf, D. R.; Smith, R. D. Characterization of Supercritical Fluid Solvents Using Solvatochromic Shifts. J. Phys. Chem. 1986,90,3022-3026. Zerah, G.; Hansen, J.-P. Self-consistent integral equations for fluid pair distribution functions: another attempt. J. Chem. Phys. 1986,84,2336. Received for review September 12, 1989 Revised manuscript receiued December 5, 1989 Accepted January 5, 1990
Reynolds Shear Stress for Modeling of Bubble Column Reactors Thomas Menze1,t Thomas in der Weide, Oliver Staudacher, Ondra Wein,$and Ulfert Onken* Fachbereich Chemietechnik, Lehrstuhl Technische Chemie B, University of Dortmund, Postfach 500500, 0-4600Dortmund 50, FRG
A general model for the design and scale-up of bubble columns and airlift loop reactors requires the exact prediction of their flow structure. The deterministic and stochastic components of the liquid flow in these gas-liquid contactors have therefore been measured. With the aid of a newly developed measuring technique, the hot-film anemometry with triple split probes, the local mean velocities, the turbulence intensities and, for the first time in bubble-column reactors, Reynolds shear stresses were determined. On the basis of a large amount of experimental data, a hydrodynamic model for the radial profile of the mean axial liquid velocity and an empirical correlation for the prediction of the turbulence intensities have been developed. The hydrodynamic model and the correlation are valid for various reactor geometries, gas and liquid flow rates, and various properties of the liquid. 1. Introduction
Because of their simple construction and ease of maintenance, bubble columns and airlift loop reactors are widely used as gas-liquid contactors. Their application covers typical chemical processes such as hydrogenation and oxidation of hydrocarbons, as well as aerobic fermentations of various types in the field of biotechnology (Deckwer, 1985). The widespread application leads to various constructions and to reactor sizes that extend from a few liters, as in the growth of plant or animal cells (Katinger et al., 1979), to thousands of cubic meters, as in wastewater treatment (Zlokarnik, 1982). 'Present address: Bayer AG, ZF-TVT 4, Bayerwerk, D-5090 Leverkusen, FRG. *Presentaddreas: Institute of Chemical Process Fundamentals of the CSAV, Suchdol, CS-16OOO Prague, Czechoslovakia. 0888-5885/90/2629-0988$02.50/0
Although bubble columns and airlift loop reactors have been the topic of a large number of publications in the last decade, no reliable general model for their design and scale-up has been developed until now. The main reason for this lies in the lack of data about their complicated flow structure. The description of the flow in these reactors has often been based on visual observations or on measurements with techniques that do not meet the special requirements of two-phase flow. Therefore, older models disregard the flow structure or use very simplified assumptions (Joshi and Shah, 1981). In the past years, progress in the field of computer design and data processing induced the development of a number of advanced measuring techniques for the investigation of two-phase flow (Schiigerl et al., 1985;Buchholz, 1987). While a visual observer of bubble columns must draw the conclusion that the flow is of a totally stochastic nature, a closer examination using one of the novel 0 1990 American Chemical Society
