Ind. Eng. Chem. Res. 1987,26, 855-857 Registry No. NO, 10102-43-9; NO,, 11104-93-1.
Literature Cited Chen, S. L.; Heap, M. P.; Pershing, D. W.; Martin, G. B. Fuel 1982, 61, 1218. DeSoete, G. G. Private Communication, Newport Beach, CA, Jan 1980. Giammar, R. D.; Hopper, D. R.; Webb, P. R.; Radhakrishnan, E. Proceedings of the Third Stationary Source Combustion Symposium; 1979; EPA: Washington, DC, 1979 Report EPA-600/779-050a (NTIS PB292539). Glass, J. W.; Wendt, J. 0. L. Nineteenth Symposium (International) on Combustion;The Combustion Institute: Pittsburgh, PA, 1983. Langsjeon, P. L.; Burlingame, J. 0.; Gabrielson, J. E. Report EPA600/7-81-111a (NTIS PB82-115312),1981; EPA, Washington, DC. Levy, J. M.; Chan, L. K.; Sarofim, A. F.; Beer, J. M. Eighteenth Symposium (International) on Combustion; The Combustion Institute: Pittsburgh, PA, 1981. Maloney, K. L. Proceedings of the 1982 EPAIEPRI Joint Symposium on Stationary Combustion NO Control;E P A Washington, DC, 1985; Report EPA-600/9-85-022b (NTIS PB85-235612). Munro, J. M. Ph.D. Dissertation, University of Utah, Salt Lake City, 1982. Quartucy, G. C.; Buening, H. J.; Yang, R. J. Proceedings of the 1985 Joint Symposium on Stationary Combustion NO, Control;1986; E P A Washington, DC, 1986; Report EPA-600/9-86-021b (NTIS PB86-225059).
855
Rees, D. P.; Smoot, L. D.; Hedman, P. 0. Eighteenth Symposium (International) on Combustion; The Combustion Institute; Pittsburgh, PA, 1981. Starley, G. P.; Bradshaw, F. W.; Carrel, C. S.; Pershing, D. W.; Martin, G. B. Combust. Flame 1985a 59(2), 197. Starley, G. P.; Bradshaw, F. W.; Carrel, C. S.; Pershing, D. W.; Martin, G. B. Combust. Sci. Technol. 1985b,43, 21. Starley, G. P.; Slaughter, D. M.; Munro, J. M.; Pershing, D. W.; Martin, G. B. Nineteenth Symposium (International) on Combustion; The Combustion Institute: Pittsburgh, PA, 1983. *Author to whom all correspondence should be addressed.
Gregory P. Starley, D. Robert Brodbeck David W. Pershing* Department of Chemical Engineering University of Utah Salt Lake City, Utah 84112 G. Blair Martin Air and Energy Engineering Research Laboratory Environmental Protection Agency Research Triangle Park, North Carolina 27711 Received for review September 27, 1985 Accepted January 12, 1987
Nonlinear Diffusion in the Solid Phase A Fourier series expansion method has been used to develop simple and relatively accurate solutions for nonlinear Fickian diffusion in the solid phase. For the fundamental harmonic, the first term in the Fourier expansion is qualitatively similar but quantitatively quite different from the conventional exponential function. In contrast, the Fourier coefficients of the higher order harmonics are combinations of fast-decay and slow-decay components. The simple analytical results agree well with numerical solutions for one sample calculation. Analytical solutions for transient-state Fickian diffusion and for heat conduction in the solid phase are well-known when the diffusion coefficient, y, is a constant. These solutions have been reviewed in detail by Carslaw and Jaeger (1959), Crank (1975), and others. However, the situation can be extremely complicated in the case of variable properties (for example, concentration-dependent diffusion coefficients are common for many sorbentpolymer systems). A number of solutions have been reviewed by Crank (1975), and they can be quite complex because of the nonlinearity of these partial differential equations. For many applications, it is convenient to have analytical solutions that are both simple and relatively accurate. It has been pointed out previously (Tsang, 1961) that the Fourier series expansions are still valid for these nonlinear problems. However, the Fourier coefficients are no longer purely exponential functions of time. These coefficients cannot be evaluated exactly. In the present work, approximate methods have been developed to calculate these Fourier coefficients. It may be seen that for the fundamental harmonic, the Fourier coefficient is qualitatively similar but quantitatively quite different from the pure exponential functions. In contrast, the Fourier coefficients of the higher order harmonics are linear combinations of fast-decay and slow-decay exponents. Such results are compared with numerical solutions for a one-dimensional system in the following. The agreement appears to be quite good. The method is also applicable to nonlinear diffusion and heat conduction, differential equations in general. Consider the diffusion in a slab of thickness L , and large area, with the initial uniform sorbent concentration, &, 0888-5885/87/2626-0855$01.50/0
and in contact with an environment of sorbent concentration $o. The concentration, 4, inside the slab is then given by the Fick’s diffusion equation
where z is the distance, T is the time, and the diffusion coefficient, y, may depend on the concentration, 4. De) and introducing dimensionless variables, fining ~ ( 4 =~yo, x = m / L , t = K 2 y o T / L 2 , C = (4 - 40)/(dl- 40), and D = y/yo, then the diffusion equation (1)may be written as
ac = at
”(.g) ax
with the initial condition, C = 1at t = 0, and the boundary conditions, C = 0 at x = 0 and at x = a. The specific example, D = 1 + XC, will be used where X is a constant. The set of functions, (2/a)lI2 sin nx (where n = 1, 2, etc., is an integer), is orthonormal and complete for 0 < x < a and goes to zero at the boundaries x = 0 and x = a. Any arbitrary function, F ( x ) , which satisfies the boundary conditions ( F = 0 at x = 0 and x = K ) can be expanded by the Fourier series F(x) = E,$, ( 2 / ~ ) sin ~ /nx ~ where the Fourier coefficients, f,, are constants. The dimensionless concentration, C ( t ,x ) , is a function of two variables, t and x . The variable, t , may be considered at first as a fixed parameter (Lowdin, 1959). Because of the boundary conditions, C(t, 0) = 0 and C(t, a) = 0, Fourier coefficients may be calculated at any arbitrary time, t. Thus, the Fourier coefficients,f,, are functions of the dimensionless time, t , 0 1987 American Chemical Society
856 Ind. Eng. Chem. Res., Vol. 26, No. 4, 1987
c ( t , x ) = C f , ( t ) ( 2 / ~ ) ' / sin ~ nx
(3)
n
Because of the symmetry about the middle of the slab, the summation is over the odd integers, n = 1, 3, 5, etc. Substitution of (3) into the diffusion equation (2) gives
+ n2fnDsin nx -
C(2/?r)l12[(dfn/dt)sin nx n
close to its initial value, (8/?r)1/2.If fl is replaced with ( 8 / ~ ) ' /then ~ , (13) may be written as dgk/dk = -k2[1
dfk/dt = C[akn(t) + Pkn(t)]fn(t) n
(5)
where k and n are odd integers and ffkn
= (-2n2/r)$TD sin kx sin nx dx 0
= (2n/?r)ST(dD/dx) sin kx
gk(t) = C eXp(-k2[1+ ( 8 X / ~ ~ ) ] t ]
D=1
-k2[1
(7)
+ XC = 1 + Xf1(2/x)'/2 sin x + ...
(8)
COS
and the initial conditions are f k = (8/ak2)'/' at t = 0. For constant y or = 0, then P k n = 0, cikk = -k2, and (Ykn = 0 for k # n. Thus, the solution of (5) has the very simple )'/~ for the linear exponential form fk(t) = ( 8 / ~ k ~ exp(-k2t) diffusion equation, where k is an odd integer. Although (5) is an infinite matrix equation for X # 0, it can be simplified by considering only the two dominant terms, the diagonal term ( n = k) and the term coupled to f l ( n = 1). The result is particularly simple for f l where the sum in (5) is over n = 1 only: dfi/dt = (ai1 + 8ii)fi = -fi[l
+ (4X/3~)(2/~)"'fiI (9)
With the initial condition, f l = ( 8 / x ) ' / ' at t = 0, the solution is
fl(t) = (8/x)1/2(3?r2e-t)/(3?r2+ 16X - 16Xe-t) (10) This function differs considerably from the usual exponential form that occurs in the case of constant y or X = 0. The higher order Fourier coefficients, such as f3, f5, etc., may now be calculated. Even larger deviations from exponential functions have been found. Including the two dominant terms, n = 1 and n = k, on the right side of (5) and then using (6)-(8) to calculate the matrix elements, cy and P, we get dfk/dt =
(akk
+ Pkk)fk +
+ Pkl)fl
dfk/dt = -k2[1 + (8/?r3)1/2kfl]fk + (32/~~)'/'[k/(k' - 4)]Xfl2 (11) in which a very small term, -k2(8/?r3)1/2(4k2- 1)-'/'Xf j k , has been neglected. This is a linear differential equation for f k where f l is given by (10). The general solution is the sum of a complementary function (CF), gk, and a particular solution (PS),hk: fk(t)
= gk(t)
hk(t)
(12)
For gk(t), the second term on the right side of (11) is omitted dgk/dt = - k 2 [ l
+ (8/?r3)'/2Xf1]gk
(13)
Since k2 >> 1, the function, gk(t), decays more quickly with time than fl does. For the short-time region where &(t) differs significantly from zero, the function f l is relatively
(15)
where c is the integration constant. The function, hk(t),is approximately proportional to f 1 2 and behaves approximately as e-2t;thus, we have the approximate solution, dhk/dt = -2hk, which is quite small compared to the first term on the right side of (11). If the dhk/dt term is neglected, then (11) may be written as
(6)
nx dx
0
(14)
The general solution is
nfn(dD/dx) cos nx] = 0 (4)
After multiplication by ( 2 / ~ ) ' /sin ~ kx and integration from x = 0 to x , the orthonormal property gives
(8X/~')]gk
+ (8/?r3)'/'Xfl]hk + (32/x3)'/'[k/(k2 - 4)]Xf12 = 0 (16)
The solution is hk(t) = ( 3 2 / ~ ~ ) ' / ~ X f l ~ / (-k4)[1 ( k ~+ ( 8 / ~ ~ ) ' / ~ X f(17) l]) The general solution for fk is given by the sum of (15) and (17): fk(t)
+ hk(t) = C exp(-k2[1 + ( 8 X / x 2 ) ] t )+ (32/a3)'/'~f1'k-'(k2 - 4)-'[1 + ( 8 / ~ ~ ) ' / ~ X f l (18) ]-'
= gk(t)
The integration constant, c, is determined by the initial condition that f l = ( 8 / x ) ' f 2 and f k = ( 8 / ~ k ~ )at' /t ~= 0. The result is c = (8/?rk2)'iz(2/x)'/2(32/~2)k-'(k'- 4)-'[1
+ (8X/x2)]-'
(19)
For k = 3, the function, f 3 ( t ) ,is f 3 ( t )= (0.53 - [1.55X/(9
+ 7.2X)])e-(9+7.2A)t + 0.61Xf12/(9 + 4.6Xf1) (20)
The dimensionless concentration, C, may be obtained from (3) by using f l from (10) and f3, f5, etc., from (18). The average concentration, CA (averaged over the distance, x ) is then
where the summation is over odd integers. These analytical results may be compared with numerical solutions. Equation 2 has been solved numerically for X = 10. For a given time, t , the numerically computed concentrations, C, have been Fourier-analyzed. The first and third Fourier coefficients, fl(t)and f3(t),are shown by the crosses and circles in Figure 1 on semilogarithmic scales. The average concentration, CA, is also shown by the solid triangles. The open triangles (with spikes) are the CA values that were calculated numerically by Crank (1975). The curves represent the analytical results from (lo), (201, and (21). The general agreement is quite good. From these semilogarithmic plots, we note the gradual change in the slope of the f l curve. In contrast, there are two distinct regions of rapid and slow decays for f3 due to the interaction of the CF and PS terms. It may be concluded that simple analytical solutions may be obtained for the nonlinear Fick's diffusion equation by an eigenfunction expansion method. These solutions appear to be reasonably accurate. The choice of the eigen-
I n d . Eng. Chem. Res. 1987, 26, 857-859
857
D = dimensionless diffusion coefficient f = Fourier coefficient (function of time) F = arbitrary function g = complementary function h = particular solution k = odd integer L = thickness n = odd integer t = dimensionless time x = dimensionless distance z = distance from surface Greek Symbols a = matrix element /3 = matrix element y = diffusion coefficient X = dimensionless constant
'
0.01 0
0.2
0.4
0.6
t Figure 1. Fourier coefficients,fi(t) and f3(t),and average concentration, CA, for D = 1 + 10C or X = 10. Solid curves: analytical solutions from (lo), (20), and (21). Crosses, circles, and solid triangles: numerical solutions of Fick's diffusion equation. Open triangles with spikes: numerical solutions for CAgiven by Crank (1975).
function, ( 2 / 1 r ) ' /sin ~ nx, was dictated by the slab geometry and the boundary conditions. For geometries such as cylinders and spheres, different types of eigenfunctions are necessary (Crank, 1975).
Acknowledgment The partial financial support by National Aeronautical and Space Administration Grant NAG-5-156 is gratefully acknowledged.
Nomenclature c = integration constant C = dimensionless concentration CA = average dimensionless concentration
T = time 4 = concentration of sorbent &, = sorbent concentration in the environment = initial sorbent concentration
Literature Cited Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solids; Clarendon: Oxford, 1959; Chapters 3-5. Crank, J. The Mathematics of Diffusion, 2nd ed.; Clarendon: Oxford, 1975; Chapters 4-9. Lowdin, P. 0. Adu. Chem. Phys. 1959,2, 260. Tsang, T. J. Appl. Phys. 1961, 32, 1518.
Tung Tsang* Department of Physics and Astronomy Howard University Washington, D.C. 20059 Carol A. Hammarstrom Code 6120, Naval Research Laboratory Washington, D.C. 20375 Received for review November 21, 1985 Accepted December 12, 1986
A Consistent Rule for Selecting Roots in Cubic Equations of State In several papers by Robinson (1979), Peng and Robinson (1976), and Soave (1972), it has been taken for granted when selecting roots of cubic equations of state (EOS) that the largest root is for the compressibility factor of the vapor, while the smallest positive root corresponds to that of the liquid. Some recent studies by Lawal (1985) reveal that this assumption is only valid for van der Waals (VDW) EOS. Therefore, the purpose of this communication is to illustrate the reason why this assumption does not always hold for all cubic EOS of van der Waals type. Both the analytical and graphical techniques are presented to illustrate a consistent rule to be applied in choosing the smallest positive root of EOS. Examples of applying the rule in Peng-Robinson (PR) EOS for the high-density, single-phase fluid are illustrated. T h e error in predicting the fluid density with the selection of an incorrect compressibility factor resulted in an average absolute percent error of 17 142%. Let us consider a cubic EOS of the form a(T) p = - RT u - b v 2 abu pb2
+
+
Equation 1 can be rewritten as 2 3 - [I + (1 - 413122 + [ A - a~ + (p - 4~212 [AB + p ( ~ + 2 ~ 3 1 1=
o
(2)
where
B = -bP
RT
(4)
0888-5885/87/2626-0857$01.50/0
By assigning certain integral values to a and 0, eq 1 and 2 take the same form as the familiar cubic EOS. Examples of these are VDW (a = 0, p = 0, a = a ( T ) ) , RKS (a = 1, p = 0), PR ( a = 2 , p = -l),and Harmens (a = 3, /3 = -2). Equation 2 yields one or three roots depending upon the temperature of the pure component. In the two-phase region as dictated by the pressure being less than the vapor pressure when T < T,,the largest root is for the compressibility factor of the vapor. The smallest positive root, which must be greater than or equal to B of eq 4, corresponds to that of the liquid. This is the criterion for testing the roots of eq 2, but when this criterion is ignored, the use of the absolute value becomes as necessary in the fugacity equation as it is in the GPA program (Bergman, 1976). Therefore, 'complex roots and roots with a value 0 1987 American Chemical Society
