EXPERIMENTAL TECHNIQUE
DIFFUSION COEFFICIENTS IN HYDROCARBON SYSTEMS Homogeneous Phases at Elevated Pressures G . R . G A V A L A S , H . H . R E A M E R , A N D B. H . S A G E Chemical Engineering Laboratory, California Institute of Technology, Pasadena, Calif. 97709
The phenomenon of change in pressure accompanying gaseous diffusion in a closed isochoric, isothermal system has been studied experimentally and theoretically. Special equipment has been developed for measuring pressure differences with an accuracy better than 0.01 p.s.i. a t total pressure levels of hundreds o f pounds per square inch. The theoretical analysis shows that as the system approaches equilibrium, the -br)lzt pressure as a function of time assumes the simple exponential form p - pes= -ae , where pes i s the equilibrium pressure, D l 2 is the binary diffusion coefficient, a i s a thermodynamic quantity, and b i s a constant depending on the size and shape of the diffusion cell. This expression can be used for the determination of Dlzfrom pressure-time measurements. Measurements upon the system CH4-C3Hs a t J = 160' F. and pes = 570 p.s.i. have shown the applicability of this method, but more experimental work i s required for ascertaining its accuracy.
IFFUSION
coefficients in homogeneous systems involving
D both gases and liquids have been measured for many years (Jeffries and Drickamer, 1954; Maxwell, 1890). I n many cases these methods entailed measurement of the composition at a given point in a homogeneous system as a function of time. These techniques have yielded satisfactory results, but, except when radioactive isotopes are employed (Timmerhaus and Drickamer, 1951), they necessitate the inconvenient withdrawal of samples at several times during diffusion. The sample withdrawal usually involves some disturbance of the system. A method has been developed to use a measured change in pressure with respect to time during diffusion in a closed isochoric system under isothermal conditions. The change in pressure experienced is related to the deviation of the system from an ideal solution (Lewis, 1960). I n principle, the system is brought to a given temperature and pressure and a quantity of one of the components of the binary system is introduced a t one end of the uniform system. The pressure is measured as a function of time as equilibrium is approached. Theoretical
Real gases exhibit a pressure change upon mixing at constant weight, total volume, and temperature. At low pressures, this pressure change is small, but can be appreciable a t higher pressures because of pronounced deviations from ideal solutions (Lewis, 1960). A typical case is the diffusion of a binary mixture 1-2 in a diffusion cell of length z, as shown schematically in Figures 1 and 2. Initially, the composition described by the mole fraction X I of component l is a step function, a condition that experimentally can be realized only approximately. At any later time the mole fraction profile is smooth and tends asymptotically with time to the equilibrium value ~ 1 , ~The ~ . initial part of the process is accompanied by a rapid pressure change. The rate of pressure change rapidly 306
l&EC FUNDAMENTALS
- 0.5
8
0.4
0
+ U 0.3
4
a
LL
0.2 J
0
= 0.1 0.2
0.4 0.6 REDUCED POSITION
0.8
z/z,
Figure 1. Schematic representation of composition profiles during a binary diffusion experiment
decreases and the equilibrium value, peq,is approached asymptotically with time. The rate of pressure change, as a measure of the speed of the diffusion process, can be used for determining the diffusion coefficient of the components of the binary mixture. The evolution of the one-dimensional composition profile and the pressure with time is described by the following equations of change in molar units: ac
at
+ a-az (cv*)
c -a x 1 + c v * - az =
at
az
a2 -
=
0
( 2) Dl*C
(4)
equilibrium. Thus Equation 9 should be used simultaneously with Equations 1 through 4. T h e conservation of component 1 provides an additional equation:
TIME
where ~ 1 is the , ~ known ~ equilibrium concentration of component 1. Equation 10 is redundant and can be used conveniently as a check of the numerical calculations. Two methods of solution are presented. T h e first is numerical and as a practical matter assumes a compositionindependent diffusion coefficient; otherwise it is general. T h e second method applies only for large times when the system is near equilibrium. I t provides an analytical expression for the pressure-time relationship and can be used effectively for the determination of the diffusion coefficient from the experimental data of the pressure as a function of time. Since the system is near equilibrium, it gives values of the Chapman-Cowling diffusion coefficient for the equilibrium pressure and composition of the system.
SEC
Figure 2. Schematic representation of change during a binary diffusion experiment
pressure Numerical Solution for All Times (Constant
with boundary and initial conditions:
(5) XI(Z, 0) = f ( z ) (6) T h e various symbols are defined in the nomenclature, which is that of Bird et al. (1960). T h e molar average velocity, u*, negligible a t low pressures, can be appreciable a t high pressures because of the concentration dependence on composition. However, u* is smaller than the momentum or mass average velocity and this is the reason for writing the continuity equations in molar quantities. Equation 3 is the equation of state of the mixture. Equation 4 is used in place of the rriomentum equation because equalization of pressure takes place much more rapidly than the diffusion process. T h e pressure, as a function of time, still remains a n unknown variable. Assuming for the moment that xl(z,t) is known, the pressure a t a given time can be calculated from the equation:
c
(XI(Z,
)
i'),p(t), T dz
= ceqzo
The numerical calculation is essentially a procedure whereby a diffusion coefficient is assumed and the pressure as a function of time is calculated from Equations 1 through 6. The calculation can be repeated with various values of the diffusion coefficient until the best possible fit is obtained between calculated and experimental pressure-time curves. Although this procedure can be used to determine the diffusion coefficient, its application is not convenient. One difficulty is the formulation of a mathematical initial condition, such as Equation 6, which describes well the experimental initial condition, and another is the composition dependence of the diffusion coefficient, which complicates the numerical solution of Equation 2. The calculation requires the use of partial volumetric data which often are not available. For the foregoing reasons the numerical solution has been used mainly for checking the range of validity of the analytical asymptotic solution to be presented. Since it also has some interest from the numerical point of view, it is briefly presented below, while the details are given in an appendix. First, the following transformation is introduced:
(7)
where ceq, the equilibrium concentration, is a known quantity. Equation 7 is simply an equation in the unknownp(t). I n most situations, the total pressure change is small and we can linearize Equation 3 to:
XI =
sb;'
c
( x ' , p ( t ) , T ) dx'
which for fixed p establishes a 1 to 1 correspondence between the variables x 1 and XI.T h e derivatives of XI and XIare related as follows:
1'
( x ' , p ( t ) , T ) dx'
c - = - -
at
at c-
where po is a given fixed reference pressure-for example, the initial or the equilibrium pressure. This expression is introduced in Equation 7 to obtain:
D12)
ax1 = a x1 az bz
Thus Equations 2, 5, and 6 can be rewritten as:
z=o,z=zo
a x1
az
=o
Equation 9 is the basis for all subsequent calculations. However, x1 and p cannot be calculated separately except near VOL 7
NO. 2
MAY 1960
307
where Table 1.
The advantage of the above transformation lies in the fact that Equation 1 4 is the standard diffusion equation with a nonlinear source term, G, which is generally small. The relative magnitude of G for the systems investigated is initially about 20% of the other terms and decreases with time as equilibrium is approached. The details of the numerical solution of Equation 14 are presented in the appendix. Figure 3 compares calculated and measured pressure-time curves for the experimental conditions of Table I. I n the experiment, component 1 was injected as a liquid into a cell occupied by pure gaseous component 2 at a pressure above the critical pressure of the binary system for the temperature in question. A short time after the injection there is only a single phase, whose composition profile can be estimated only crudely for the purpose of starting the numerical solution. The calculations utilized the Benedict-Webb-Rubin equation of state with the coefficients given by Opfell et al. (1959). A diffusion coefficient, Dlz = 0.596 X 10-5 sq. ft. sec.-l, was used in the calculations. This value has been obtained from the experimental curve by the analytical asymptotic technique described below. Considering the uncertainty of the initial conditions and the composition dependence of the diffusion coefficient, the agreement between experimental and calculated curves in Figure 3 is very good. Asymptotic Analytical Solution
For large times the state of the system is near equilibrium and the following simplifications are permissible :
1. T h e convective term, c v * ( b x l / d z ) , may be neglected as being the product of two small terms. 2. The diffusion coefficient, D12,may be considered equal to its value at the equilibrium composition. Even when the over-all composition changes by 10 to 207, from z = 0 to z = z o , it is possible to use, as a space-average value, the diffusion coefficient at the equilibrium composition. 3. The effective concentration may be taken as the spaceaverage or equilibrium value.
Conditions of Binary Diffusion Experiment
Component 1. Propane Component 2. Methane Weight of component 1. 0.0023898 Ib. Weight of component 2. 0.0044650 Ib. Volume of cell. 0.00357898 cu. ft. Length of cell, to. 5.986 inches Temperature, T . 160" F. Equilibrium composition, XI, e q . 0.1630 Pressure at end of run. 569.969 p.s.i.
Assumptions 1 to 3 are equivalent to a linearization of Equation 2 around equilibrium. The linearized equation is:
T h e general solution of Equation 18 with the boundary conditions given in Equation 5 is:
The coefficients a, depend on the initial conditions and need not be specified a t this point. Near equilibrium the pressure change is an effect much smaller in magnitude than the composition change. For example, the numerical solution discussed above shows that a change in composition, [xl(O, t l ) xl(0, ~ z ) ] / x I , = ~ ~ 0.135 corresponds to a change in pressure, [p(tz) - p(tl)]/p,, = 0.00042. This effect is also implied by conditions 1 to 3 for linearization and permits calculation of the composition change without regard to the pressure change. Subsequently, the pressure change can be calculated from the composition change as follows: Near equilibrium we may expand c in Taylor series and neglect terms of order 3 or higher:
Integration of this equation between 0, zo gives:
(21) The first and third terms in the numerator can be evaluated by integrating Equation 19 : PZ.
I
I
I
I07
IO8
I
I
I09 TIME
I10
I&EC FUNDAMENTALS
1
SEC
Figure 3. Experimental and calculated curves for the methane-propane system 308
I I I I x 103
pressure-time
r----i
Figure 4. equipment
Schematic
r---- 1
representation
of
a t pes - p = 6 p.s.i. and becomes smaller as equilibrium is approached. By using Equations 22, 2,3, and 24, Equation 21 is reduced to:
/ 3%1
This can also be written a$:
-
where tl is any reference time such that p p ( t J is small. By taking logarithms in Equation 25 there is obtained:
By arranging the experimental data in accordance with Equation 27 we can obtain the diffusion coefficient, &, from the slope of the resulting line, graphically or numerically. This calculation does not require information on the initial condition or on partial volumetric data because the constant al and the
(
(
$)eq
I
I
experimental
At large times, the second and higher terms in the series are negligible compared to the first. I n fact, the series converges much more rapidly than the series in Equation 19 because each exponential is squared. T h e second term in the numerator of Equation 21 has been found from the numerical calculations to be much smaller than the third term. For example:
thermodynamic quantities
I
determine only the
intercept of the line of Equation 27. Experimental
T h e binary mixture of methane and propane was confined in vessel A of Figure 4 and the heavier component was introduced through valve B at the bottom of this vessel to prevent convection. Valve C connecting vessels A and D was closed and diaphragm E was used as a linear transducer to indicate when the pressure in A was identical to that in D . T h e details of this equipment have been described by Reamer and Sage (1 9%).
T o adjust the pressure in D to that in A , a second vessel partially filled with mercury was employed. Valve F remained open and mercury was introduced into vessel G from injector H, which was actuated through nut J by a motor K. Nitrogen a t approximately the same pressure as the binary mixture in A filled vessels D, L, and G. From the indications of a counter attached to motor K , the change in volume of this system to maintain the pressure in D equal to that in A was established. Such volumetric measurement permitted a precise measure of the change in pressure in A to be realized. T h e basic pressure level of the system was determined by means of a balance, M , provided with a piston-cylinder combination (Reamer and Sage, 1958), which as shown in Figure 4 was connected to the system through the oil-mercury interface in the auxiliary vessel, X. A manual injector, P,was provided to adjust the elevation of the oil-mercury interface in vessel N . Valve Q was closed after the pressure level had been established, in order to avoid any changes in the quantity of mercury in the system except from the movement of injector H. Vessel AD and associated valves were located in an agitated silicone bath, the temperature of which was controlled from the indication of the strain-free platinum resistance thermometer through a suitable modulating circuit (Reamer and Sage, 1958). The actual temperature was measured by means of a second strain-free platinum resistance thermometer (Reamer and Sage, 1958), which had been compared with the indication of a similar device calibrated at the National Bureau of Standards. I t is believed that the temperature was known within O.0lo F. of the International Platinum Scale and that variations in temperature were less than 0.001' F. over a 24hour period. Vessel GL was located in a similar agitated bath. The entire system was calibrated at the pressure level to be employed in the experiment by determining from the pistoncylinder combination balance at A4 the change in pressure corresponding to the full travel of injector H. For the work carried out to date a total travel from 3000 revolutions of motor K corresponding to 30,000 counter readings resulted in a change in pressure of approximately 20 p.s.i. in system DGL. The pressure in A was adjusted to correspond to that in D by the admission or withdrawal of gas. The change in pressure was established from a knowledge of the counter reading and the volumetric behavior of the nitrogen at the temperature and pressures in question (Sage and Lacey, 1950). T h e change in pressure corresponding to a change in counter readings could be easily established for each pressure level. T h e sensitivity of the slack diaphragm was approximately 0.001 p.s.i. I t was found possible to detect and quantitatively evaluate changes in pressure of the order of 0.002 p.s.i. The accuracy of measuring differences in pressure is of the order of 0.5% or 0.002 p.s.i., whichever is the larger measure of uncertainty. A sample set of measurements involving the methane-propane system a t 570 p.s.i.a. a t a temperature of 160' F. is shown in Figures 3 and 5. T h e conditions associated with the experiVOL. 7
NO. 2
MAY 1960
309
I I 0 CALCULATED
\-
-0
1
I
:
0.8
0.6 e
0.4 e 0,
2 0.2
u
s
0
IO00
2000 TIME
3000 SEC
Figure 8. Comparison of Equation solution of Equations 1 through 6
TIME
Figure 6.
SEC
Residual drawing of data of Figure 3
TIME
SEC
Figure 7. Comparison of Equation 27 with data ane-propane system
for meth-
Straight lines obtained by least squares
ment are given in Table I. The equipment employed is relatively stable with respect to time and permits measurements of changes in pressure to be established within 0.03 p.s.i. over a period of 12 hours. The larger uncertainty over the longer periods of time is associated with the stability of the electronic circuit utilized in connection with the linear transducer associated with measurements of the position of the slack diaphragm, E, of Figure 4. Discussion of Results
Because of the linearization of Equation 2 it is known that Equation 27 is asymptotically valid near equilibrium, but its exact range of validity can be determined only by numerical 310
l&EC FUNDAMENTALS
27
with numerical
calculations and experiment. If Equation 27 is valid in the range peq - Ap to peq and if Sp is the standard error in relative pressure measurements, Equation 27 is useful for interpreting experimental data only when Sp/Ap is small-e.g., 0.05 or less. Another experimental question regards the accuracy of measuring pes, on a relative basis, and the time needed for this measurement. The experiment of Table I and the numerical calculations will now be discussed relative to the experimental requirements stated above. Figures 5 and 6 show the experimental pressuretime curve and Figure 7 shows the same curve plotted in accordance with Equation 27. Equation 27 appears to agree with the experimental data a t pressures from 2 to 0.5 p.s.i. of equilibrium. Section ab of the curve in Figure 5, corresponding to times smaller than 108 X lo3 seconds, is removed from equilibrium by more than 2 p.s.i. and is not represented well by Equation 27. I n section cd, corresponding to times larger than 111 X lo3seconds, slope dp/dt remains very nearlyconstant. Thus, section cd is not represented well by Equation 27. This behavior may be due to some disturbance, such as a drift in the electronic circuit or a vibration, and does not permit a very accurate determination of peq, The middle section, bc, is represented well by Equation 27. The standard error in this section has been computed as Sp = 0.0108 p.s.i., so that the criterion Sp/Ap < 0.05 is easily satisfied. Figure 8 shows the plot of a complete numerical solution of Equations 1 through 4 using the method outlined. The calculation was performed with constant diffusion coefficient and for the experimental conditions of Table I. Excellent agreement with Equation 27 can be observed up to 3.3 p s i . from equilibrium. This corresponds to a composition profile with x,(O,t) = 0.413, xl(z,,,t) = 0.003. The experimental curve fits Equation 27 only up to 2 p.s.i. from equilibrium. This difference is probably due to the fact that the computational model, ignoring the composition dependence of the diffusion coefficient, is more nearly linear than the physical system. Since the equilibrium pressure, p,,, was not obtained with sufficient accuracy in the present experiment, the data in Figure 7 are plotted for several values of peq. A calculation was performed to determine simultaneously the values of Diz and peq which give the best straight-line fit in the least squares sense. This calculation has given p e s = 569.88 p.s.i., Dl2 = 0.596 X 10-5 sq. ft. set.-' The sensitivity of the computed value of D12 to the selected value ofp,, is such that a change in peq of 0.1 p.s.i. results in a 10% change in D,z. This sensitivity is the main limitation in the accuracy of Dl2 as determined from
Equation 27, since the present equipment does not allow a sufficiently accurate measurement of pes. An alternative procedure for the determination of D12 without the need of estimating or measurin pes can be based on the relations
The vector and matrix quantities of this equation are defined as follows: I
A =
(35)
r
2
0
0
.
.
l
r
1
0
.
.
L o o 0 0
.
. .
0
0
0
+ -
The time interval 6 should be chosen such thatp(t 6) p(t) is much larger than the measurement error. By performing calculations with different values of 6 one can obtain a measure of the accuracy of the method. A linear regression analysis, using least squares, has been performed on the data arranged in accordance with Equation 29 with the following results: 6 = 200 seconds; D1z = 0.697 X lO-5sq. ft. sec.-I; 6 = 1000 seconds; 2 3 1 2 = 0.715 X 10-5 sq. ft. sec.-l No other experimental data are available for comparison. The value of DIZ estimated by the Chapman-Enskog formula and corrected for elevated pressure by using pseudocritical properties and the generarlized chart of Bird et al. (1960) is Dl2 = 0.395 X 10-5 sq. f i . set.-' The method based on Equation 29 seems preferable, since it eliminates the need of estimating p e s and the concomitant sensitivity problem. Future work will concentrate on improving the accuracy of the method by extending the useful interval, bc, of Figure 5 to pressures closer to equilibrium and by using a more systematic statistical analysis. Acknowledgment
T h e help of Virginia Berry with the statistical calculations is appreciated. The authors are indebted to one of the reviewers for suggesting the use of the alternative procedure, known in chemical kinetics as Guggenheim’s method (Frost and Pearson 1965). Appendix
The numerical solution of Equations 14, 15, and 16 is based essentially on the Crank-Nicholson implicit finite difference scheme (Lapidus, 1962). Thus letting:
t
=
lAt;
t
= mAz;
NAz = zo
(1 = 0 , 1, . . . ; m = 0, 1,
X(z,t) = Xm,i, x?(z,t) =
Xm-1, Z + I
X I , ml,
etc.
+ .E,+I,z - 2 X m , + X m - I ,
11
(32)
Equation 11 in finite difference form is written as:
AXl+i
=
-BXz
-
bGz
1 = 1,2, , . .
.
.
0
2
~
1
where
q = - Z ( l + - i T 1 ) . 1 = - 2 ( 1 - ~ ) , 6 = ~2 Az2 912At
DizAt
Dl2
(3 7) I n order to solve Equation 33 for X L + lin terms of X z we should first calculate Gi from the data at previous times. This calculation is performed as follows: 1. A preliminary value for p at time I& is chosen-e.g., the same value as at t = ( I - 1)At. Then Equation 11 can be used to obtain a table for interpolating betneen x 1 and X . From this table Lve obtain by interpolation X I , ~ J , . , x ~ , . ~ ~ . T h e concentration profile, c o t , . . , c.y~ is calculated from Equation 3. 2. A new improved value for the pressure is calculated by using Equation 9 and calculations 1 are repeated. 3. Steps 1 and 2 are repeated until the pressure changes negligibly between two iterations. The final pressure obtained in this fashion corresponds to time lAt. Usually, only two to three iterations are necessary. 4. From data at times (At, ( I - l)At, and ( I - 2)At we may compute the derivative bc, bt and calculate GI. 5. With G J known, Equation 33 is solved by any one of the methods of solving linezr algebraic equations. Nomenclature
matrix defined by Equation 35 first coefficient in the series, Equation 20 matrix defined by Equation 36 quantity defined by Equation 37 total molar concentration molar concentration of component 1 in binary mixture 1-2 equilibrium value of total molar concentration equilibrium value of concentration of component 1 initial mole fraction profile of component 1 quantity defined by Equation 17 molar flux of component 1 molar flux of component 2 pressure fixed reference pressure equilibrium pressure quantity defined by Equation 37 quantity defined by Equation 37 temperature time molar average velocity defined by u* = 1 (.Vi
+ Y*)
VOL 7
NO. 2
(33) MAY 1 9 6 8
311
x l ( z , t ) = mole fraction of component 1
X(z,t) 2
zo 912
= = = =
variable related to x, by Equation 11 position along diffusion cell length of diffusion cell Chapman-Cowling binary diffusion coefficient of components 1 and 2, defined by Bird et al. (1960) and Longwell and Sage (1965) :
SUBSCRIPTS 1 = component 1 2 = component 2 eq = equilibrium
Jeffries, Q . R., Drickamer, H. G., J . Chem. Phys. 22, 436 (1954). Lapidus, L., “Digital Computation for Chemical Engineers,” McGraw-Hill, New York, 1962. Lewis, G. N., J . A m . Chem. Soc. 30,668 (1960). Longwell, P. A., Sage, B. H., A.Z.CI1.E. J . 11,46 (1965). Maxwell, J. C., “Scientific Papers,” p. 625, University Press, Cambridge, England, 1890. Opfell, H. B., Pings, C. J., Sage, B. H., “Equations of State for Hydrocarbons,” American Petroleum Institute, New York, 1959. Reamer, H. H., Sage, B. H., Rev. Sci. Instr. 29, 710 (1958). Sage, B. H., Lacey, W. N., “Thermodynamic Properties of the Lighter Paraffin Hydrocarbons and Nitrogen,” American Petroleum Institute, New York, 1950. Timmerhaus, K. D., Drickamer, H. G., J . Chem. Phys. 19, 1242 (1951). RECEIVED for review October 12, 1967 ACCEPTED March 14, 1968
Literature Cited
Bird, R. B., Stewart, W. E., Lightfoot, E. N., “Transport Phenomena,” Wiley, New York, 1960. Frost, A. A., Pearson, R. G., “Kinetics and Mechanism,” Wiley, New York, 1965.
Work supported in part by American Petroleum Institute Project 37 a t the California Institute of Technology.
MEASUREMENT OF TRANSIENT CHANGES IN LIQUID HOLDUP IN PACKED BEDS NICHOLAS S T A N D I S H
Department of Metallurgy, Unicersity College, Wollongong, N . S. W., Australia
An apparatus is described which records the transient behavior of liquid irrigation in packed beds with and without a countercurrent gas flow. The apparatus measures weight changes of the bed by means of strain gages mounted on a cantilever beam. The signal from the gage bridge is amplified and recorded on a chart, Examples of recorder traces show a number of interesting features of holdup and loading phenomena. The apparatus, for which an accuracy of better than 1 % is claimed, could also be used in conjunction with transient mass transfer studies.
KNOWLEDGE
of liquid holdup in packed beds is necessary
A from the standpoint of packed tower design and is also important from the viewpoint of mechanism of mass and heat transfer. Many investigators have studied the various aspects of liquid holdup and tower irrigation but none appears t o have considered the transient behavior of these phenomena. A number of investigators, notably Shulman et al. (1963), have indicated the importance of holdup for providing qualitative and quantitative treatment of the transient behavior of packed columns. These authors, however, did not study such transient holdup in columns. Numerous correlations of liquid holdup in packed beds have been reported. Unfortunately, these correlations were obtained using mean or time averaged results and cannot, therefore, indicate directly the unsteady state behavior of liquid irrigation. Similarly, the more recent investigations of transient behavior of packed columns using pulse, step, or frequency response techniques, while yielding valuable data on liquid dispersion, do not give a complete picture of the transient flow of the liquid irrigant. Therefore, an apparatus has been developed which will continuously record the transient nature of liquid irrigation of packings. 312
l&EC FUNDAMENTALS
Apparatus
The basic arrangement of the apparatus is shown in Figure 1. I t consists of a metal beam, M , one end of which is firmly secured to a rigid support, B , while the free end carries a sus-
d r h
Figure 1 . paratus
Diagram of ap-
pended cage, C, containing the packed column, P, and manometer, U. A set of strain gages, S, is mounted on the beam a short distance away from B. Two gages on the tension side and two on the compression side of the beam are wired to form a resistance bridge. T h e output leads of the bridge are connected to a
