Ind. Eng. Chem. Res. 1987,26, 456-462
456
for di = 2.58 mm, u = 0.466 m/s, h!f = kg of c, t b L = 103 s for Beypazari lignite, tbCh = 120 s for Beypazari char (400 "c),and tbC = 124 s for Beypazari coke (900 "c).As seen from these data, the burn-out time for coke particles is about 20% greater than that for coal particles. It is believed that this result is due to the rapid devolatilization of the coal particles in situ; thus, a more porous structure is left behind, and oxygen can access easily through the pores. Another reason may be the higher heating rate of the char particles due to combustion of volatiles, compared with coal particles which do not initially contain volatiles. If the volatiles are heated to temperatures in excess of about 900 "C as soon as they are released, under conditions such that there is an oxygen deficiency, part of their carbon content is converted to soot (Field et al., 1967). In this study, observation of the luminous flame during the combustion of volatiles may be attributed to soot formation because there is an oxygen deficiency in the bed, lasting for a short time after the introduction of coal particles. Therefore, the combustion of volatiles in a FBC with soot formation needs further investigation, since it is important in heat transfer by radiation.
Conclusions The burn-out time of volatiles increases with the masses charged, particle size, and volatile matter content and decreases with the air velocity. It was found that the surface flux of volatiles increases with particle size. The radial velocity of the gases coming out of the particles during the devolatilization process is near or greater than U,. This may be the reason why coal particles float on the bed surface during combustion of their volatile materials. The subsequent effect of the volatiles on the remaining char particles is to increase the rate of combustion of these char particles. Acknowledgment I am grateful to Prof. J. F. Davidson for giving me an opportunity to carry out this work at the Department of Chemical Engineering, University of Cambridge, and for his valuable comments. I also thank Dr. R. D. LaNauze of the National Coal Board for providing the samples of Arigna coal and M. Ghadiri and K. V. Thambimuthu for their valuable comments. The financial support of a WHO
fellowship from the World Health Organization is acknowledged. The views expressed are mine and not necessarily those of WHO.
Nomenclature d, = initial particle diameter ( d = 2R) M = mass of the coal charged N = number of coal particles in a batch n = number of moles of volatiles Q = volumetric flow rate of volatiles r = radial distance from the center of a particle T = temperature of the bed t b = total burn-out time of coal particles t , = burn-out time of volatiles U = superficial velocity at T Umf= minimum fluidization velocity at T V = volume of volatiles emitted from coal particles ur, u8, u m = velocity of volatiles in r, 0, and C$ directions Greek Symbols p , p c = apparent and carbon density of coal particles Literature Cited Avedesian, M. M.; Davidson, J. F. Trans. Znst. Chem. Eng. 1973,51, 121. Basu, P. Fuel 1977, 56(10),390. Basu, P.; Broughton, J.; Elliott, D. E. Fluidized Combustion; Symposium Series 1; The Institute of Fuel: London, 1975; paper A3. Bird, R. B.; Steward, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 1960; p 83. Campbell, E. K.; Davidson, J. F. Fluidized Combustion; Symposium Series 1; The Institute of Fuel: London, 1975; Paper A2. Chakraborty, R. K.; Howard, J. R. J. Znst. Fuel 1978 (Dec), 220. Chakraborty, R. K.; Howard, J. R. Chen. Eng. Commun. 1980, 4 , 705. Field, M. A.; Gill, D. W.; Morgan, B. B.; Hawskley, P. G. W. Combustion of Pulverized Coal;The British Coal Utilization Research Association: Leatherhead, 1967; Chapter 5, p 175. Howard, J. B.; Essenhigh, R. H. Eleventh Symposium (Znternational) on Combustion; The Combustion Institute: Pittsburgh, 1967; p 399. LaNauze, R. D. Fuel 1982, 61, 771. Pillai, K. K. J . Znst. Energy 1981, 54, 142. Ross, I. B.; Davidson, J. F. Trans. Znst. Chem. Eng. 1981, 59, 108. Ross, I. B.; Patel, M. S.; Davidson, J. F. Trans. Znst. Chem. Eng. 1981, 59, 83. Yates, J. G.; Walker, P. R. In Fluidization; Davidson, J. F., Keairns, D. L., Eds.; Cambridge University Press: New York, 1978; p 241.
Received far review May 30, 1985 Accepted August 19, 1986
Turbulent Diffusion of a Traveling Wave of Charged Aerosol Richard M. Ehrlich and James R. Melcher* Laboratory for Electromagnetic and Electronic Systems, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
The effective turbulent diffusion coefficient of a duct flow is deduced from measurement of the longitudinal mixing of a stream of aerosol particles, charged by a sinusoidally varying corona source. The air-stream turbulence conditions in an electrostatic precipitator (ESP) can play an important role in determining its overall collection efficiency (Feldman et al., 1976; Marietta and Swan, 1976). In order to predict the behavior of such devices, one must adequately model the effects of turbulent mixing on particle distributions. It is 0888-5885/87/2626-0456$01.50/0
common to account for these effects by the use of an empirically determined equivalent diffusivity (Pyle et al., 1980; Leonard et al., 1980; Friedlander, 1977). This paper presents a specialized technique for measurement of the effective turbulent diffusivity of a duct flow, based upon the longitudinal mixing of a stream of 0 1987 American Chemical Society
Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987 457 particles charged by a sinusoidally varying corona source. The technique and an experiment are described in the next section; the data are in the following section. Subsequent sections present a continuum model for the turbulent mixing of a passive scalar contaminant (in this case, the particle charge density, p ) , some of the more important approximations, and values for the equivalent diffusivity, D,, found for the various flow conditions. The concept of accounting for flow turbulence via an equivalent diffusivity is often attributed to Taylor (1921). Early models for electrostatic precipitation (Deutsch, 1922) pictured the effective flow diffusivity, D,, as so large as to eliminate any transverse particulate concentration gradients in the precipitator. In recent years, efforts have been directed toward modeling the behavior of precipitators in the presence of moderate diffusion so that turbulent mixing is significant but not the dominant particle transport mechanism. Leonard (1982) provided a review of the theoretical and experimental work that led from the Deutsch model for precipitation to his own, which he termed a “convective diffusion model”. Assuming a flow with uniform mean velocity, effective turbulent diffusivity, and particle migration velocity, he solved analytically for the particle concentration as a function of position in the idealized precipitator. To verify these predictions, he injected precharged particles into a collection region with a uniform imposed electric field and measured the concentration profiles in situ by using an optical particle counter. The flow turbulence was varied by the use of grids and baffles placed upstream of the precipitation region. With flow parameters similar to those used in this research ( H = 5 cm, U N 3 m/s), he obtained good agreement between theory and experiment. The effective transverse diffusivities fitted to the data ranged from 6 X lo4 to 3 X m2/s, depending upon the imposed flow conditions. Here, the emphasis is placed upon longitudinal mixing of particles in a flow that has no imposed transverse particle migration velocity. Otherwise, the model is consistent with Leonard’s.
Charging and Detection Configuration The experimental setup and procedure used for the diffusivity measurements are similar to those used previously by us to observe the statistics of ac particle charging (Ehrlich and Melcher, 1984a,b). The major difference between the two lies in the data reduction techniques used. In the previous experiments, only time-independent statistical quantities, such as the mean and root mean square (rms) spread of the migration velocity, were interpreted. No velocity-vs.-time plots or velocity-frequency spectra were generated, as turbulent mixing downstream of the charger partially destroyed the coherence of the ac charging wave form. It is precisely this decay in the coherence of the particle charge wave form that is to be measured here and used to deduce the turbulent diffusion coefficient. A monodisperse aerosol of dioctyl phthalate (l-pm diameter) is mixed with air before it enters a flowstraightening section at the entrance to a small tunnel (Figure 1). The channel is composed of metal plates of width 14 cm on two sides and glass walls of height 3 cm on the other two. The ac charging potential is applied to a linear array of steel pins which are suspended on a brass rod extending across the width of the channel, halfway between the metal plates. The geometry and applied potentials are contrived so as to permit corona ion generation only near the tips of the pins. Ions that carry the corona currents from pins to plates charge the incoming
particles, a small fraction of which are precipitated immediately. It is assumed, for the purposes of the diffusivity computation, that the particles leave the charging region with a mobility that varies (from particle to particle) sinusoidally with time, at the frequency of the ac corona excitation, wo. A small fraction of ions is entrained in the gas and will neutralize the aerosol as it is carried downstream. To prevent this, the pins face upstream so that the field originating downstream of the tips serves to precipitate the ions. (They are, of course, far more mobile than the aerosol.) The dc ion trap also serves to ensure that these entrained ions are removed. As the flow proceeds downstream, turbulent mixing introduces progressively greater amounts of disorder into the sequence of particle charges. It is assumed that the charge and identity of each particle remains unaltered. Only their ordering changes. Figure 2 illustrates idealized typical progressions of particle charge vs. time (for fixed observers) at locations immediately after the charger (Figure 2a), and for points farther downstream (Figure 2c,e). The Fourier transforms correspondingto these charge-vs.-time functions are parts b, d, and f, respectively, of Figure 2. Note that the same particle charge statistics are valid at all three locations. Only the ordering is changed by the turbulent mixing. The Fourier component at the base frequency shrinks,however, as the coherence of the particle ordering decays. The charge plots of Figure 2 were generated by beginning with a sinusoidal function, sampled at regular intervals, and then adding a random component to the time of each data point. In this way, the charge statistics were preserved. The Fourier transforms were produced by the same computer programs that were used to analyze the data to be presented here. The charged particles in the flow eventually pass through a region in which a thin aluminum plate is suspended halfway between the grounded walls. With this plate, used to impose a dc electric field (E N 5 X lo6V/m), a laser-Doppler velocimeter (LDV) is then used to determine the particle mobility by measuring the transverse particle velocity induced by the imposed electric field. Because the aerosol is monodisperse, this is tantamount to a measurement of the charge. The details of the LDV data processing system are described elsewhere (Ehrlich, 1984). For present purposes, it is sufficient to say that the recorded information from a single run of the experiment consists of roughly 3000 ordered pairs of measured migration velocities vs. time of the measurement (ui,ti). The time span of a run was typically about 2 s. Once a set of measured particle velocity-time pairs is obtained, the Fourier transform, F ( o ) ,is approximated as 0 L
N-1
-
F(o)= -~ - u i e j w q t i + l - ti] t N - tl i=1
where N is the total number of pairs (ui,ti) obtained. The factor in front of the summation in eq 1 ensures that the transform of a sinusoidal velocity of frequency wo will be very nearly equal to the magnitude of that sinusoid when w = wg. Due to the discrete nature of this formula for F(w) and the finite measurement time span, it is useful only for frequencies, w , which are low enough so that the mean particle interarrival time, ( t N - t l ) / N ,is much smaller than 2 ~ / wbut high enough that 2 r / w is much smaller than (tN - tl).
If the particle charge were truly a sinusoidal function of time, then the Fourier component of measured migration velocity at the charging frequency, ufeXc = F(w,J,would
458
Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987
A P i n - bar
Flow st ra iphtener
I
\
I
/
4
I
Exhaust f a n
M iqra t i o n plate
/
\ I\ I
\
m t e r control l e d relay
I
I
-
\
I I
\
Ion t r o p
HeNe loser Hiqh-voltaqe DC supply
- -
Transmitting optics
T I I
I
LDV observation
/ velum
U U
Series blockiq c a w c i tor
H qh- v d t q e A t supply Rece I v I nq OPtKS
L D V /' slqnal P ocesior and computer
Figure 1. Side and top views of the wind tunnel. Aerosol, entrained at the left, sequentially passes the pin bar (where it is sinusoidally charged) and the ion trap (where entrained ions are removed before they can neutralize the aerosol) through a mixing region and finally through a dc analyzer (the migration plate) where its mobility is measured by means of the laser Doppler velocimeter.
be exactly equal to 2lI2 times the rms spread in velocity, u'. Thus, the ratio y = uf.exc/21/2u'istaken as a measure of the decay in coherence due to turbulence. In the section on the model, a relationship is derived between the decay factor, y, the effective diffusivity, D,, and other physical parameters of the flow.
0bservations Figures 3 and 4 illustrate Fourier transforms of migration velocity data. For these, the horizontal air-flow velocity was 2.5 m/s and the electrical excitation frequency was w0/27r = 20 Hz. For the first one, the pin tips were about 100 cm from the observation volume, while for the second the pins were "near" the LDV observation volume (about 25-cm separation). Although the rms spread in the measured migration velocity was nearly the same for both [15.7cm/s for Figure 3,16.1 cm/s for Figure 41, the Fourier component of velocity at the charging frequency was obviously much smaller for the former. Figures 5 and 6 show typical observations of the variation of the decay coefficient, y, with excitation frequency, f. Respectively, these are with the charging pins in the "near" and "far" positions. Figure 7 illustrates the frequency dependence of the ratio between the decay coefficients with the pins in the two positions. The points on that figure were computed using "smoothed" versions of the data from Figures 5 and 6. The
solid curve corresponds to the variation of this relative decay that is predicted by the model to be considered next.
Model The situation is as pictured in Figure 8. A turbulent air flow of uniform mean velocity, U , and effective diffusivity, D,, carries neutrally buoyant particles of number density, n. As particles pass through the plane at z = 0, they are charged by a time-varying (corona) source so that the charge, q , on a particle that is at z = 0 when the time is t is given by where Qo is the amplitude of the charging resulting from the charging process and wo is its angular frequency. The particles are assumed to act as inert markers in the flow. Effects of self-precipitation (due to repulsion between particles of like charge) and recombination (due to attraction between particles of opposite sign) are neglected, with good reason. It is desired to predict the time dependence of the particle charge at some fixed location downstream from the charger. Specifically, the following analysis predicts the behavior of the particle charge density, p ( z , t ) , averaged over a volume, AV, of dimensions much larger than the mean interparticle spacing, n-'I3, but much smaller than any macroscopic dimensions of the
Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987 459
I
h / s d
I
OD
0.0
E 8.0
(Hzl Figure 4. Fourier transform of measured transverse particle velocity with an excitation frequency of 20 Hz and mean and rms flow velocities of 2.5 m/s and 16.1 cm/s, respectively, and with the corona pins ''nearn the migration plate. f
w-
W
Figure 2. Theoretical illustration of progression of particle charge with time (a) immediately downstream of the charger and (c-f) for points further downstream. I
I
I
go0
I
00
0.
I
Q Q
0 0
t
Y
Q Q
O.
Q
Q
Q
0
0
0.0
IO0
50 f (Hz)-
101 f
(Hz)
Figure 3. Fourier transform of measured transverse particle velocity with an excitation frequency of 20 Hz and mean and rms flow velocities of 2.5 m/s and 15.7 cm/s, respectively, and with the corona pins "far" from the migration plate.
system [such as the duct height, 2H,or the excitation wavelength, X = 2 ~ U / w ] .
l M p(r,t) E -Eqi AV;=1
(3)
In eq 3, q i is the charge on the ith particle, with a total of M charged particles contained in the volume, AV, which is centered downstream at F. The equation of charge conservation, in the face of convection and turbulent diffusion, is (4) where
is the unit vector in the z direction. It is assumed
Figure 5. Decay factor as a function of charging frequency for corona pins "near" the migration plate and mean velocity of 7.5 m/s.
that no significant gradients of flow or charging parameters exist in either transverse direction. Equation 4 becomes
aP(z,t) ap(z,t) at + U - = D az, - - g -
a2p(z,t)
(5)
In the steady state, the response of the system will be harmonic at the excitation frequency, wo,
p(z,t) = PO Re @ (z)ejwot1
(6)
where z(g) satisfies the (normalized ordinary) differential equation
The complex density, z(g), is dimensionless, having been normalized to p o = Qon,the complex particle charge density existing at z = 0. The distance is normalized to the length, L , of the mixing region. Characteristic diffusion and
460 Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987
I
I
I
0
1
1
I
I
I
I
I
0
0
+I
y
0.2
0 0
0
I : / 01 0 0
0
0
-
0
0
0
8
0.05
0.00 . 0
20
40
0 00.0
60
3 0
(Hz)-
f
G-
Figure 6. Decay factor as a function of charging frequency for corona pins “far” from the migration plate and mean velocity of 7.5 m/s.
Figure 9. Plot of eq 11.
convection times, based upon that length, are rD= L 2 / D , and r C5 L / U. The boundary conditions imposed on the normalized complex charge density are p ( ~ ) l = ~ =1 ~and
P(z)Iz-m
m.
The solution to the elliptic ordinary differential equation, eq 7, subject to the given boundary conditions, is =
A
(8)
where the normalized wavenumber,
\
0.6C
$2~
determines the rate of spatial decay of the particle charge density. The decay coefficient, y, is given by the ratio of charge density a t z = L to that a t z = 0 y = Ip (1)l = e”l!l (10) -
or f (Hz)--
Figure 7. Ratio of decay factors of “far” pins to that with “nearn pins as a function of charging frequency. The flow velocity was 7.5 m/s, and the solid line is based on the theoretical model with a m2/s. turbulent diffusion coefficient of D,= 1.7 X Charqer
//
Particles
I \
I
I
Figure 8. Schematic depicting longitudinal turbulent mixing where particles charged in the plane z = 0 suffer turbulent mixing as they are convected the length L.
where G 4w0rC2/rDis a dimensionless constant. The form of eq 11is convenient for relating experimentally measured values of y to values for the equivalent flow diffusivity. Figure 9 shows a plot of the quantity on the right side of eq 11 as a function of G. Given values for y, wo, and T for a particular experiment, one can use the figure to find the corresponding value of G, which then allows determination of D,. For small values of D , G N 4 In (l/y)/(worc) and D, If3 In ( l / y ) / ( w ; L ) . Since the particle charge density decays exponentially with the distance downstream from the charger, one can, in principle, compute a decay Coefficient by comparing p ( z ) between any two points separated by the distance, L. It is not necessary to measure p ( z ) a t the charger itself. Discussion of Model Both the approximation of the particle cloud as a continuum of noninteracting elements and the assumption of a sinusoidal charging process are addressed briefly in this section. In the experiments reported here, the typical particle densities used were n = 2.5 X lo9 mW3.Thus, the mean
Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987 461 Table I. Summary of Turbulent Diffusion Coefficients
U, m/s
Yr
2.5 4.3 7.5
0.432 0.291
D,,m2/s 1.09 x 10-3 8.04 x 10-3 17 x 10-3
interparticle spacing, n-lI3 N 0.7 mm, was much smaller than the duct height, 2H N 30 mm. With a 2 m/s flow rate and 20-Hz charging, X 300 mm. Thus, it is reasonable to treat the cloud of charged particles as a continuum. It is argued by Ehrlich (1984) that the charge repulsion and recombination processes occur on time scales that are of the order 11 7 e-18 tTev where = 4 / 3 ~ R 3n, is the volume fraction of the flow occupied by the particles, and reve q/EE? is the electroviscous time for the charged flow (E,, is the charging electric field, q is the viscosity of the air, and e is the permittivity). With q = 2 X lod, Eo= l@, e = 8.85 X W2, R = 0.5 X lo4 [MKS units], and n as given previously, r p 80 s, which is much longer than the convection time defined earlier, TC 0.5 s. Thus, it is reasonable to ignore effects of particle repulsion and attraction. The corona charging of particles is a rate-dependent process (Melcher, 1981). In the case of a practical ac corona charger, driven by a sinusoidal voltage, the particle charge observed at the charger outlet is likely to contain significant harmonics of the fundamental charging frequency, cog. Thus, to represent the harmonics, eq 2 might be replaced by a Fourier series. Because of the assumption of negligible particle interactions, each charging harmonic could be dealt with separately, using the approach of the previous section. For the sake of simplicity, the effects of such charging harmonics have been ignored here. It can be shown that, for small values of the dimensionless parameter, G, the characteristic decay length of the particle-charge mode falls in inverse proportion to the square of the charging frequen cy.
Application of Longitudinal Mixing Model to Experimental Data Table I contains results from the application of the model to the experimental data. For flow velocities of 2.5 and 4.3 m/s, the ratio of the values for the decay coefficients obtained with the charging pins far from the migration plate to those obtained with the pins close to it was measured at only the one frequency of 20 Hz. These values resulted in the diffusion coefficients shown for the fiist two velocities in Table I. For the highest velocity, when wo ranged from 2 ~ 2 to 0 2 ~ 4 0the , plot of Figure 9 was used, along with the definitions of TC and TD and L = 75 cm. For this velocity, measurements of the relative decay coefficient were made for frequencies ranging from 20 to 40 Hz. The value of D, = 17 X m2/s was obtained by fitting a curve of the predicted decay coefficient to the actual data. The solid curve on Figure 7 shows the variation of the relative decay coefficient with frequency as predicted by the model, using this value of D,. The effective diffusivities obtained range from 1 to 17 x 10-~ m2/s. Conclusions The flow diffusivity measurement technique introduced here produced reasonably self-consistant results, with
diffusion coefficients, D,, in the range 1-17 X lo3 mz/s for mean flow velocities of 2.5-7.5 m/s. This is in rough agreement with values reported by Leonard et al. (1983) and Robinson (1975). These diffusivities can thus be used in quasi-one-dimensional models which predict the behavior of devices such as electrostatic precipitators, in which turbulent mixing is superimposed upon orderly particle migrations.
Acknowledgment This work, which formed part of a thesis submitted by R.M.E. to the Department of Electrical Engineering and Computer Science at M.I.T. in partial fulfillment of the requirements for the ScD, was sponsored by Combustion Engineering Inc., with Gerald Driggers both administering the project and functioning informally as a thesis committee member. Formal members of the committee were Profs. J. Haritonidis, L. Smullin, and M. Zahn. Support also came in the form of a 1-year Cabot Solar Energy fellowship.
Nomenclature D, = turbulent diffusivity, mz/s E = electric field intensity, V/m Eo = charging electric field, V/m f = frequency, Hz F ( w ) = Fourier transform G = 4wO7c/?D = half-channel height, m i, = unit vector in z direction
j = (-1)lP k = wavenumber, m-l
L = length of mixing region, m qi = charge of ith particle, C q article = particle charge, C = amplitude of charging source, C F = position, m t , t N , ti = time, s U = mean gas velocity, m/s u, ui = particle velocity, m/s ~ f = velocity . ~ ~ at excitation ~ frequency, m/s u' = rms spread in velocity, m/s A V = observation volume, m3
do
Greek Symbols y = decay coefficient, eq 10 q = gas viscosity, kg/(ms)
F/m charge density, C/m3 rc = convection time L / U , s rD = turbulent diffusion time L2/D,, s rev = electroviscous time q/eEo2,s w = angular frequency, s-l wo = excitation angular frequency, s-l 6 = particle volume fraction, m3 t
= 8.85 X
p =
Literature Cited Deutsch, W. Ann. Phys. 1922,68, 335. Ehrlich, R. M. Sc.D Thesis, M.I.T. Cambridge, MA, 1984. Ehrlich, R. M.; Melcher, J. R. "AC Electrostatic Precipitation", 1984a, "AC Corona Charging of Particles", 1984b, Presented at the Annual Institute of Electrical and Electronic Engineers Industrial Applications Society Conference, Chicago, IL. Feldman, P. L.; Kumar, K. S.; Cooperman, G. D. AIChE Symp. Ser. 1976, 73(165), 120. Friedlander, S. K. Smoke, Dust and Haze:. Wilev: " New York. 1977: Chapter 3. Leonard. G. L. High TemDerature Gasdvnamics Laboratorv ReDort 196, NSF Grant No. 'CPE-7926290: and EPRI Contract 'No. EPRI-RP-533-1, 1982; Mechanical Engineering Department, Stanford University, Palo Alto, CA. Leonard, G . L.; Mitchner, M.; Self, S. A. Atmos. Enuiron. 1980, 14, 1289.
462
I n d . Eng. Chem. Res. 1987,26, 462-469
Leonard, G. L.; Mitchner, M.; Self, S. A. J.Fluid Mech. 1983,127, 123. Marietta, M. G.; Swan,G. W. Chemical Engineering Science; Pergamon: Elmsford, NY, 1976; Vol. 31, p 795. Melcher, J. R. Continuum Electromechanics; M.I.T. Press: Cambridge, MA, 1981; pp 5.2-5.16. Pyle, B. E.; Pontius, D. H.; Snyder, T. R.; Sparks, L. E. Presented
at the 73rd Annual Meeting of the Air Pollution Control Association, Montreal, Quebec, 1980. Robinson, M., PhD Thesis, Cooper Union, NY, 1975. Taylor, G. I. Proc. London Math Soc., 1921, A20, 196.
Receiued for review June 3, 1985 Accepted July 17, 1986
Solid-Vapor Azeotropes in Hydrate-Forming Systems Jashwantsinh L. Thakore and Gerald D. Holder* Chemical and Petroleum Engineering Department, University of Pittsburgh, Pittsburgh, Pennsylvania 15261
Three-phase vapor (V), water-rich liquid (LJ, and hydrate (H) equilibrium conditions and solid hydrate azeotropic compositions were determined for various ternary mixtures in the temperature range from 274.15 t o 281.15 K. The experimental measurements demonstrate that the systems, methane-propane-water and krypton-propane-water, form hydrate azeotropes in the temperature range 274.15-281.15 K. The systems, methane-cyclopropane-water and methane-isobutane-water, do not form hydrate azeotropes in this temperature range. Azeotropes occur as a result of the hydrate crystal structure and are not caused by hydrate-phase nonidealities. Gas hydrates are crystalline compounds which form from mixtures of light, nonpolar gases and water. In particular, hydrate crystals can form from mixtures of natural gases (methane, ethane, and propane) and water. Hydrates can form in pipelines (Hammerschmidt, 1934) or in association with an underground hydrocarbon reservoir, where they can impede production by blockage of rservoir pores (Katz, 1971, 1972; Holder et al., 1976). Because hydrates can form at temperatures well above the freezing point of ice, they have been considered as a method of removing salt from sea water (Barduhn et al., 1962; Barduhn, 1969) and as a potential vehicle for the storage of natural gases (Parent, 1948). In any of these applications, the conditions determining the stability of a hydrate are pressure (P),temperature (T) and equilibrium gasphase composition (y). Any thorough evaluation of the conditions of hydrate stability requires knowledge of the extreme limits of such stability. This limit will be the maximum temperature at a given pressure or the minimum pressure at a given temperature for which hydrates are stable. Most thermodynamic studies of gas hydrates have focused on the conditions where a vapor (V), a water-rich liquid (L), and a solid hydrate (H) are all in equilibrium. The pressure and temperature conditions describing such three-phase equilibria are univariant-a line on the P-T coordinates as in Figure 8-for any binary gas-water mixture or for any multicomponent mixtures of fixed gas-phase composition. For a specified gas mixture, such equilibrium conditions determine the limits of stability (maximum temperature or minimum pressure) for the hydrate phase. The limiting composition conditions for hydrate formation have not been thoroughly studied however. Thermodynamic studies have shown, in general, that larger hydrate-forming molecules, such as propane, form hydrates with relatively low equilibrium dissociation pressures, and smaller molecules, such as methane, form hydrates with relatively high dissociation pressures. Hydrates from mixtures of such larger and smaller molecules will have dissociation pressures nearer to the dissociation pressure of the larger molecules than might be expected. For example, a 99% methane-1% propane gas mixture, when contacted with water, will form hydrates with a dissociation
* Address
t o whom correspondence should be addressed.
0SSS-5SS5/S7/2626-0~62$01.50/0
pressure which is less than half the dissociation pressure of pure methane (Deaton and Frost, 1949). A few studies (Snell et al., 1961; Verma, 1974; van der Waals and Platteeuw, 1959; Holder and Grigoriou, 1980) have shown that some ternary (gas + gas + water) mixtures can, in fact, have dissociation pressures which are lower than those of either gas-water binary pair which comprise the ternary mixture. For example, at 278 K a methanepropane mixture containing 25 mol% propane in the gas phase forms hydrates a t a pressure which is about 10% lower than the hydrate-formation pressure of pure propane (Verma, 1974). For this ternary system, there will be some mixture which will have a lower dissociation pressure than any other mixture; this pressure will be lower than the dissociation pressure of both propane and methane hydrate. Such a mixture will be an azeotropic mixture, but the azeotrope’s existence does not correspond to nonidealities in the classical equilibrium model but rather to the nature and structure of crystalline hydrates in general. The greater stability (as evidenced by lower equilibrium pressures) of the hydrates formed from certain mixtures is due to competing effects. For example, in propane hydrates, propane molecules can enter only the large cavities of structure I1 hydrate, and as in all hydrate phases, the hydrate stability is derived from its cavities being occuppied. Methane, a relatively small molecule, normally forms structure I but can enter the small cavities of structure I1 and contribute to the stabilization of that structure in a way that propane cannot, since propane does not enter the small cavities and therefore can never stabilize them. The net effect of decreasing the mole fraction of propane in a propane-methane mixture is to increase the stability of the hydrate provided by methane and to decrease the stability provided by the propane. Under certain conditions, including the specific point mentioned above, the increase in stability by methane is greater than the decrease in stability by propane, and a lower pressure is required to form the hydrates. If propane could enter the small cavities in a crystal structure, the increase in methane mole fraction would result in a higher, not a lower, equilibrium pressure.
Model The stability of a crystalline hydrate phase depends on the fraction of its cavities occupied by gas molecules. If an insufficient fraction is occupied for a given set of 0 1987 American Chemical Society
