C, = c,/u,, dimensionless wave velocity for the fastest growing wave f = stream function, cm2/s F = f / h m u mdimensionless , stream function g = gravitational constant, cm/$ h , = film depth for uniform flow, cm Le = entrance length, cm N R ~ = u,h,/u, Reynolds number N o = ~ ( 2 / p ~ v ~ gsurface ) ' / ~ , tension number q = volumetric flow rate per unit width, cm2/s t = time, s u r n = free surface velocity for the uniform flow, cm/s U = dimensionless velocity profile for the primary or undisturbed flow x,y = rectangular Cartesian coordinates, cm Y = y/h,, dimensionless y direction
Greek Letters a = &hm,dimensionless wave number a , = dimensionless wave number for the fastest growi
wave a, = critical wave number pr = 2 i ~ / X ,wave number, cm-1 pi = spatial growth rate factor, cm-1 7 = position of the free surface, cm X = wavelength, cm v = kinematic viscosity, cm3/s p = density, g/cm3 u = surface tension, dyn/cm o = frequency,s-l
Literature Cited Anshus, B. E., lnd. Eng. Chem., fundam., 11, 502 (1972). Anshus, 6.E., Goren, S.L., AlChEJ., 12, 1004 (1966). Benjamin, T. B., J. FluidMech., 2, 615 (1957). Cerro, R. L., Whitaker, S., Chem. Eng. Sci., 26, 785 (1971a). Cerro, R. L., Whitaker, S.,J. Colloidlnterface Sci., 37, 33 (1971b). Cerro, R. L., Whitaker, S., Chem. Eng. Sci., 29, 963 (1974). Conte, S.D., Dames, R. T., Math. Tables Aids Comp., 12, 198 (1958). Davies, J. T., Rideal, E. K., "Interfacial Phenomena", Academic Press, London, 1963. Gaster, M., Prog. Aero. Sci., 8 , 251 (1965). Graef, M., "Mitteilungen ausdem Max-Planck Inst. fur Stromungsforschung", No. 36, Gottingen, 1966. Jones, L. O., Whitaker, S.,AlChE J., 12, 525 (1966). Kapitza, P. L., Kapitza, S. P., "Collected Papers of P. L. Kapitza", Macmillan, New York, N.Y., 1964. Krantz, W. B., Goren, S.L., lnd. Eng. Chem., Fundam., I O , 91 (1971). Krantz, W. B., Owens, W. B., AlChEJ., 19, 1163 (1973). Lee, J., Chem. Eng. Sci., 24, 1309 (1969). Lighthill, M. J., Whitham, G. B.. Proc. Roy. SOC. London, Ser. A, 229, 261 (1955). Lynn, S.,AlChE J., 6, 703 (1960). Pierson, F. W., Ph.D. Thesis, University of California at Davis, Department of Chemical Engineering, 1974. Schubauer, G. E., Skramstad, H. K., NACA Rept. 909 (1949). Strainthorp, F. P., Allen. J. M., Trans. Inst. Chem. Eng., 43, T85 (1965). Strobel, W. J., Whitaker, S.,AlChEJ., 15, 527 (1969). Tailby, S.R., Portalski, S., Trans. Inst. Chem. Eng., 40, 114 (1962) Whitaker, S.,Jones, L. O., AlChE J., 12, 421 (1966). Whitaker, S., lnd. Eng. Chern., Fundam., 3, 132 (1964). Whitaker, S.,AlChEJ., 17, 997 (1971). Yih, C-S., Phys. Fluids, 6, 321 (1963).
Received for review September 7,1976 Accepted February 9,1977 This work was supported under NSF G r a n t ENG 74-13037.
Toward a Diffusion Theory of Drying Stephen Whitaker Department of Chemical Engineering, University of California, Davis, California 956 16
A recent theoretical analysis of drying in granular porous media is used as the basis for a diffusion theory of drying. Boundary conditions at the porous media-surrounding gas interface are analyzed and the previously derived transport equations for mass, momentum, and energy are combined in an appropriate manner to yield a tractable theory of drying.
1. Introduction In a previous analysis (Whitaker, 1977) of simultaneous heat, mass, and momentum transfer in porous media, the volume averaging method (Whitaker, 1967a; Slattery, 1972; Gray, 1975) was used to derive the appropriate transport equations for the drying of granular porous media. The analysis began with the following point equations
Equations 1.1through 1.3 were applied to the liquid, gas, and solid phases illustrated in Figure 1.The averaging volume, Y, illustrated in Figure 1 is used in the definition of various averages associated with the volume averaging method. The spatial average of some function $ defined everywhere in space is given by (1.4)
+ p g + KV2V DT -= kV2T + @ Dt
0 = -vp PC,
Here eq 1.1represents the ith species continuity equation in the absence of chemical reaction, eq 1.2 is the quasi-steady, creep flow form of the equations of motion for a Newtonian fluid, and eq 1.3 is the thermal energy equation when viscous dissipation and reversible work are negligible. The term @ in eq 1.3 represents the source or sink of electromagnetic energy. 408
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
where the averaging volume is the sum of the three phase volumes Y=
v, + V , ( t ) + V J t )
(1.5)
More often, in the analysis of transport processes in multiphase systems, we are interested in the average of some quantity associated solely with a single phase. For example, we may be interested in the average temperature of the liquid phase contained within the averaging volume shown in Figure 1.The point temperature in the liquid phase is designated by T g ,and the phase average of Tp is defined as
This equation represents the governing differential equation for the intrinsic phase auerage liquid temperature. The terms To and 06 are deviations defined by
Tg = ( Tg)a + To, in the (3 phase Tg = = 0, in the u and y phases vg = ( v g )
+ Qg,
in the (3 phase
vp = 06 = 0, in the 0 and y phases
(1.15a) (1.15b) (1.15~) (1.15d)
Similar equations exist for the solid and gas phases, and it is clear that a detailed accounting of energy transport during drying is a very complicated problem. However, if one is willing to accept the assumption of local thermal equilibrium the analysis is greatly simplified. This assumption can be stated mathematically as
Figure 1. Drying process in porous media.
(1.6)
( T o ) ,= ( T ) + T, ( T p ) B = ( T ) + Tg
(1.16b)
Since To is defined in the normal way in the liquid or /3 phase and is zero in all other phases, eq 1.6 reduces to
( T , ) ? = ( T ) + py
(1.16~)
aTu V ( T ) >> VT, >>x;
(1.16d)
(1.7) The phase average temperature has the drawback that if Tp is a constant, the phase average is not equal to this constant. A quantity more representative of the liquid temperature is the intrinsic phase average temperature which is given by
(1.8) Once again since TBis zero in phases other than the (3 phase, eq 1.8 reduces to 1
,
.
The volume fractions for the solid, liquid, and gas phases are defined by c, =
V,/Y;
€&)
t,(t) = V,(t)/Y
= Vg(t)/Y;
g at
(1.16a)
Here the subscript w represents u,(3, and y, and represents the deviation of the intrinsic phase average temperature from the spatial average temperature. Although rarely stated explicitly, Slattery (1972, page 404) is an exception, the assumption of local thermal equilibrium is universally imposed on multiphase heat transfer processes when conduction is the dominant mode of heat transfer. A definitive theoretical exploration of this matter is certainly needed. Upon imposing the assumption of local thermal equilibrium it can be shown (Whitaker, 1977, Sec. 111) that the governing differential equation for the spatial average temperatures takes the form TOTAL THERMAL ENERGY EQUATION:
(1.10)
Clearly the sum of these fractions is one t,
+ tJt) +
€&)
=1
(1.11)
and the phase average and intrinsic phase average are related by r p ( T 0 ) ' = (To) (1.12) In order to obtain the volume averaged form of the thermal energy equation in the liquid phase, one integrates eq 1.3over the averaging volume 'v to obtain
Sv
PC,
(s)i, dV =
koV2T6dV +
Jy@BdV
Here ( p ) is the spatial average density, C, is a mass fraction weighted heat capacity, (h)is the mass rate of evaporation per unit volume, andKT,ffis the total effective thermal conductivity tensor. Analysis of mass and momentum transport in the liquid and gas phases leads to the following set of transport equations LIQUID PHASE EQUATION OF MOTION:
(1.13)
After considerable mathematical manipulation one can show (Whitaker, 1977, Sec. 1I.C) that eq 1.13 takes the form
Ind. Eng. Chem., Fundam., Vol. 16, No. 4 , 1977
409
C"
I,'
+ ah,,,(L R1
( T ) To
1
(1.27)
and the volume constraint E,
+ E@(t)+ €,(t) = 1
(1.28)
Equations 1.17 through 1.28 represent a rigorous theoretical representation of the transport equations for heat, mass, and momentum transport during the drying of granular porous media. In the derivation of these equations no thought was given to the boundary conditions that must be imposed at the surface of the porous media. In this paper we address ourselves to the problem of boundary conditions and the diffusion theory of drying.
2. Boundary Conditions We attack the problem of boundary conditions at the porous media-surrounding gas interface in precisely the same way that the point boundary conditions were developed by Whitaker (1977). Here we will illustrate the method of attack for the mass conservation principle. We begin by considering the two-phase system shown in Figure 2. There we have designated the porous media as the $ phase, and the surrounding gas phase is designated as the a phase. The principle of conservation of mass can be stated as
In the gas phase the density, p , is identified as pa, while in the porous media it is represented by pd. We know that the continuity equation in the gas phase takes the form
-
aP, + V (p,v,)
Figure 2. Material volume containing a singular surface,
The terms pd and vd can now be identified as
+ P @ ( V @ )+ (Py)Y(V,) v,#' = 9 P P + EY(P?)Y
(2.6)
P@ = t @ P @ t r ( P - / ) Y
(2.7)
Referring to Figure 2, we integrate eq 2.2 over the volume, V,(t) to obtain
s
s
&dV+ V.(p,v,)dV= 0 (2.8) at V'X(t) Application of the general transport theorem to the first term and the divergence theorem to the second term yields
"
n
Here we have separated the surface of VJt) into the material surface, A,(t) and the area of the singular surface, Aad. Since the porous media was assumed to be rigid, we require that w nad = 0, and eq 2.9 reduces to
-
A similar analysis for the $C phase leads to =
0 (in the surrounding gas phase)
(2.2)
at where the total density, pa, and the mass average velocity are represented more specifically as
and these two equations can be added to obtain
i=N
c va = - c Pa Pa =
i=l 1 i=N
(Pi),
(2.3)
(PiVila
i=l
a
at
+ v - [P@(VP) + (/Jr)Y(vy)l =0
410
(in the porous media) (2.5)
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
-
(pdvd - pava) n,,dA @a
The total continuity equation in the porous media is obtained by adding eq 1.19 and 1.21 leading to - (5!w3+ E,(P,)Y)
+
=0
(2.12)
Here we have used the fact that Aad = Ada and na4 = -nda. We now identify the time rate of change of mass in the material volume as
and note that eq 2.1 requires the condition
Since the area, Ad,, is arbitrary this reduces to pl#lv$l-n$,= PaV,
(2.15)
n@a
This jump condition at the phase interface can be started more explicitly as i=N
C
i= 1
(pivi),.n@,
=
fplg(vp)
+ (p,)Y(v,))-n,,
(2.16)
If the mass flux in the surrounding gas phase can be represented in terms of mass transfer coefficients and density differences, eq 2.16 may be of some use to us. However, we should note in this result that we have not obtained boundary conditions for eq 1.19 and 1.21; rather we have obtained a single boundary condition for the sum of these two equations. This approach to obtaining boundary conditions at singular surfaces, i.e., phase interfaces, has been discussed in detail by Slattery (1972) and it can be applied directly to the other scalar transport equations. The analysis for the thermal energy equation is rather long, but eventually leads to (pp(Vp)Ahvap)
'n+, + (KTeff.V(T)).n+a = (kpVTp) nga (2.17) e
Here we have assumed that the temperature is continuous across the phase interface, i.e., ( T ) = T, a t the phase interface; however, eq 2.17 does allow for jump discontinuities in the temperature gradient owing to evaporation a t the surface of the porous media. In general, the flux in the surrounding gas phase would be expressed in terms of a film heat transfer coefficient leading to an expression of the type (~lg(vb)4hvap)*n,a+
=
[Pb(Vlg)
(2.19) where (ul)is the mass diffusion velocity for species no. 1. Under these circumstances the convective transport term in eq 1.22 takes the form of a diffusive term, albeit complicated by the term multiplying (u1) in eq 2.19. By following the original derivation (Whitaker, 1977) of eq 1.22 one can eventually arrive a t
a
- (eY(p1)Y) - (riz) = V (D*(')eff* V(p1)Y)
(2.20) at where the effective diffusivity tensor is now highly concentration dependent. If we now add eq 2.20 and 1.19 we obtain
( K T e f f . V ( T ) ) . n ~ , = h ( T ,- ( T ) )
where T , is the ambient gas phase temperature. Referring back to eq 1.18, we can see that the above result represents a boundary condition in both ( T ) and q, and thus couples the thermal energy equation with the liquid phase continuity equation even in the absence of internal vaporization, i.e., (ni) = 0. The boundary condition for species no. 1, the water and water vapor, can be derived by following the procedure illustrated by eq 2.1 through 2.16. The analysis leads to (PlVl), *
nature of the drying process. First we will assume that the pressure in the gas phase is a constant given by p o and we will ignore the consequences of eq 1.20. This is the type of approach that is taken in the classic Stefan diffusion tube analysis (Whitaker, 1967b) in which the laws of mechanics and the "no slip'' boundary condition are ignored, yet reasonable values are obtained for the flux of the diffusing species. Meyer and Kostin (1975) have recently reexamined the Stefan diffusion tube analysis and they find that unusual velocity profiles result when the laws of mechanics are properly applied. Nevertheless, the flux of the diffusing species is very close to that given by the classical analysis. In addition to assuming that the gas phase pressure is a constant and ignoring eq 1.20 we will follow the Stefan diffusion tube analysis and assume that the inert gas component (usually air) is stagnant. Under these circumstances the mass average velocity can be expressed as
+ (Pl)Y(V,)
+ V - ( p p ( v g ) - D*('Ieff*V(p1)y) = 0
(2.21)
Here the second term represents the divergence of the total moisture flux and it is this quantity that we can represent at the porous media-surrounding gas interface in terms of a mass transfer coefficient times a driving force. Without analysis we will simply note that eq 2.18 can be expressed as (pivi), * n b , = ( p p ( ~ l g )- D*(l)eff.V ( p i ) y ) n4,
(2.22)
Use of eq 1.18 allows us to express the total moisture flux as
- (P,)~D(')~~~'V((P~)Y .nm, / ( P ~(2.18) )')I Once again we note that eq 2.18 provides a boundary condition for neither eq 1.19 nor eq 1.22. Of special importance is the fact that we need to solve four scalar transport equations to determine ( T ) , €6, ( p l ) y , and ( p , ) Y and we have a t our disposal only three boundary conditions involving these independent variables. This situation has resulted from the fact that transport equations for water can be written for both the liquid and vapor phase; however, the flux condition a t the phase interface can only be expressed as the sum of the two transport mechanisms. To proceed we must either generate a fourth boundary condition or reduce the number of transport equations. The complexity of the problem under consideration suggests that there may be several reasonable ways of circumventing what certainly appears to be an indeterminate set of equations and boundary conditions. The following approach is suggested as one possible method of obtaining a solvable set of differential equations and boundary conditions. I t should be clear from the nature of boundary condition given by eq 2.18 that we must combine eq 1.19 and 1.22 to obtain a moisture transport equation. Before doing so we make two assumptions about the
Some new nomenclature seems appropriate at this point and the following definitions will be used K, = €lgEKpkJFb
(2.24)
K(T) = caPsFK&(T)lPlg
(2.25)
Kg = t l g ~ p t b b- py)Klg/~g
(2.26)
so that eq 2.21 and 2.23 can be combined to give
a
-((tpplg+€,(pi)~)
at
= v . ( K , * V ( e g ~ p+) K ( T ) ' V ( T )
- K,
-g
+ D*('),ff.V(pl)y)
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
(2.27) 411
convective thermal energy flux as
The fractional moisture saturation is defined as
S=
(2.28)
and we can divide eq 2.27 by the constant, p?(tp tain
+ ( P r ) ( ( C p l y ) y (vr = I(cp)p(1 - a) + ( C,),).4
[P@kp)?(V?)
fractional
+ E?),to ob-
where
K’(T) = K(T)/P?(t?+ 6 , )
(2.30)
K’g = Kg/Pp(tp + cy)
(2.31)
Here we can see the emergence of a diffusion-like equation, and arguments are given elsewhere (Whitaker, 1977) which suggest that a reasonable form of eq 2.29 is
_ a s --V * [ D - V S ]
X
)I v ( T ) (P?(q3
+ cy)
[ D - V S + K ( T )* V ( T ) - K g * g ] } . V ( T ) (2.35)
Here the scalar parameter w is defined as magnitude of the gas phase mass flux (2.36) magnitude of the total mass flux We are now in a position to summarize the proposed theory. It consists of two transport equations, one for the fractional moisture saturation and one for the temperature, for which boundary conditions are available. In addition there is a third transport equation (the gas phase diffusion equation) which will be used to calculate the mass rate of evaporation per unit volume, (riz). The governing equations are listed in Table I. There are eight equations and eight unknowns: ( T ) ,S, to, ty, (riz ) , ( P I )7, ( p 1 ) 7, and w. Here we assume that the dependence of D and D*(l)eff on composition can be represented in terms of S rather than ( p l ) y and ( p y ) y . It should be clear at this point that eq 2.37 and 2.39 are used to determine ( T )and S, but the remaining variables may be determined in several ways. w=
at
+ V . [ K ’ ( T ) * V ( T ) -] V . [ K ’ , * g ]
(2.32)
In comparing eq 2.32 and 2.29 we can conclude that the moisture diffusion tensor, D, should be dominated by capillary action during the early stages of drying and by gas phase diffusion during the latter stages of drying. This can be expressed more explicitly as
-
D * VS ---* K, * V
-
(when t y
-
Table I. Governing Equations for the Drying of a Granular Porous Media
0)
and leads us to conclude that D will undergo large changes during the drying period. A key objective in subsequent studies will be the determination of the functional dependence of D upon S. The term in eq 2.23 involving the temperature gradient does not appear in any previous discussion of drying. It results from the fact that the surface tension, uo?, is a function of temperature, and that gradients of the surface tension can give rise to fluid motion. This has been dramatically demonstrated by Trefethen; however, it remains to be seen if this is an important effect in the drying of a porous media. In terms of the representation implied by eq 2.29 and 2.32, the flux condition given by eq 2.22 is now expressed as
Although something less than rigorous, our treatment of the moisture transport problem is now complete and we can turn our attention to the matter of convective transport of thermal energy. Returning to eq 1.17 we note that the convective transport term can be expressed as
for the case where the inert component of the gas phase is taken to be stagnant. If we compare this result with eq 2.23 we see that the total moisture flux cannot be used directly to compute the convected thermal energy since the liquid and gas flow rates are weighted by their individual heat capacities. If we make the crucial but appealing assumption that the gas and liquid mass f l u x vectors point in t h e s a m e direction, we can extricate ourselves from this difficulty and express the 412
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
Definition of the Fractional Moisture Content S = t p P p + ‘r(P1)Y P?(tP
w=
+4
Ratio of the Mass Fluxes ID*(l)effV(PI)~~ p~(eB+tr)lD-VS+ K ’ ( T ) - V ( T-K‘g.gl )
(2.45)
(2.46)
Presumably at the start of any drying process the saturation and temperature would be specified. Using eq 2.42,2.43,2.44, and 2.45 one could determine the initial values of tp, e,, (PI),, and (PI),. Since initial conditions are available for S and ( T ) , we are left with only (m)and w as unspecified a t the start of the drying process. An appropriate initial condition for (& ) would be
( A )= 0
(at the start of a drying process) (2.47)
but specification of w remains something of a puzzle. The upper and lower bounds on w would be constrained by
-
w-+O w
1
(forS-1)
K’, = liquid phase gravitational mass transport tensor,
(for S - 0 )
K(T) = K’(T)PdtL?+ 4 Kg = K’gPpJtg + 4 Keff = liquid phase effective thermal conductivity tensor, kcal/s m K KTeff = liquid phase total thermal conductivity tensor, kcal/s mK KD = liquid phase thermal dispersion tensor, kcal/s m K k , = - d ( p c ) / a t o , n/m2 k ( ~ =) - d ( p , ) / a ( T ) , n/m2 K ( m ) = mass rate of evaporation per unit volume, kg/s m3 n = outwardly directed unit normal p = pressure,n/m2 p c = p , - pa, capillary pressure, n/mz p o = reference pressure, n/m2 p1O = reference vapor pressure for species no. 1, n/m2 p = p - p o p 4 , relative ressure, n/mZ q = heat flux vector, k c a l l m2 r = position vector, m r = characteristic length for a porous media which is a function of the liquid volume fraction, m Ri = gas constant for the ith species, n m/kg K = to/(t~ e.,), fractional liquid saturation = (eppp + t r ( p l ) Y ) / [ p p ( t p t,)], fractional moisture saturation T = temperature, K To = vapor pressure reference temperature, K T o = enthalpy reference temperature, K T , = ambient temperature in the surrounding gas phase, K t = time, s ui = mass diffusion velocity of the ith species, m/s U = unit tensor v = mass average velocity, m/s vi = species velocity, m/s V , = volume of the rigid solid phase contained within the averaging volume, m3 V p ( t ) = volume of a evaporating liquid phase contained within the averaging volume, m3 V,(t) = volume of the gas phase contained within the averaging volume, m3 Y = averaging volume, m3 Y , ( t ) = material volume, m3 w = velocity of a singular surface, m/s
but beyond this there appears to be no way of specifying the initial w field. From a practical point of view, it would seem reasonable to assume some value of w and then use the governing equations to calculate the w field after the first time step. This field can then be used as the initial condition and the process repeated until the assumed and calculated w fields are the same. Essentially this implies that the initial condition for the w field is
aw -. = 0
(2.48) (t = 0 ) at Assuming that an initial w field can be established, we propose the following computational scheme: (1)Use eq 2.39 to determine s. (2) Use eq 2.37 to determine ( 5 ” ) . (3) Use eq 2.44 and the previous value of €0to determine ( P I ) , . (4) Use eq 2.43 to determine ( P I ) , . (5) Use eq 2.42 and 2.45 to determine tp and e.,. (6) Use eq 2.46 to determine w . ( 7 ) Use eq 2.41 to determine ( & ) . (8) Repeat the above steps for successive time intervals. Clearly there is some extensive numerical computation required in order to determine S and ( T )as a function of time and space. However, one should keep in mind that there are only three differential equations to be solved and only one of the algebraic relations (eq 2.44) is nonlinear. The matter of experimental determination of the various transport coefficients is not at all trivial, and there is an obvious need for some theoretical basis for estimating the value of these parameters. 3. Conclusions The transport equations for heat, mass, and momentum transport in a continuous media have been used to derive the governing differential equations for drying in a porous media. While the governing equations form a complete set, there is insufficient information to form boundary conditions for all the differential equations. One method of circumventing this difficulty has been described, and the development has led to a diffusion theory of drying.
Nomenclature Lower case boldface letters designate vectors while upper case boldface letters designate second-order tensors. A = area,m2 A,(t) = closed material surface, m2 c = constant pressure heat capacity, kcal/kg K = mass fraction weighted constant pressure heat capacity, kcal/kg K D = moisture transport diffusivity tensor, m2/s D(l)eff= gas phase total effective diffusivity tensor for species no. 1, mz/s D*(l)eff = gas phase total effective diffusivity tensor for species no. 1diffusing through a stagnant inert component, m2/s D = gas phase molecular diffusivity, mz/s g = gravity vector, m2/s
8p
h = enthalpy per unit mass, kcal/kg; film heat transfer coefficient, kcal/m2 s K bo = reference enthalpy, kcal/kg hi = partial mass enthalpy for the ith species, kcal/kg Ah,,, = enthalpy of vaporization per unit mass, kcal/kg k = thermal conductivity, kcal/s m K; mass transfer coefficient, m/s KO = liquid phase permeability tensor, m2/s K, = gas phase permeability tensor, m2/s K, = liquid phase volume fraction mass transport tensor, m2/s K ’ g ) = liquid phase thermal mass transport tensor, m2/s S-1
+
2
+
+
Greek Letters t, = VJY, volume fraction of the rigid solid phase t p ( t ) = Vp(t)/‘V,volume fraction of the evapc;ra:ing liquid phase ty = V , ( t ) / V ,volume fraction of the gas phase = a dimensionless scalar function of the topology of the liquid phase = thermal dispersion vector, kcal/s m3 p = shear coefficient of viscosity, ns/m2 p = total density, kg/m3 pi = density of the i t h species, kg/m3 = rate of heat generation, kcal/s m3 4 = gravitational potential function, m2/s2 A = unit tangent vector w = ratio of the magnitude of the gas phase mass flux to the magnitude of the total mass flux up., = interfacial tension of the P-7 interface, n/m Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
413
Subscripts i = designates the ith species in the gas phase u = designates a property of the solid phase fi = designates a property of the liquid phase y = designates a property of the gas phase C$ = designates a property of the porous media phase a = designates a property of the surrounding gas phase afi = designates a property of the a-fi interface fir = designates a property of the P-y interface yC$ = designates a property of the y-a interface $a = designates a property of the 6-a interface Mathematical Symbols dldt = total time derivative D/Dt = material time derivative alat = partial time derivative ($) = spatial average of a function, $, which is defined everywhere in space ($@) = phase average of a function, $s, which represents a property of the /3 phase
= intrinsic phase average of a function, $?, which represents a property of the fi phase Note that ($), ($o), and ($-,)y are defined everywhere in space. ($.,)y
Literature Cited Gray, W. G., Chem. Eng. Sci., 30, 229 (1975). Meyer, J., Kostin, M. D., Int. J. Heat Mass Transfer, 18, 1293 (1975). Slattery, J. C., "Momentum, Energy, and Mass Transfer in Continua", McGraw-Hill, New York, N.Y., 1972. Trefethen, L., "Surface Tension in Fluid Mechanics", Educational Services, Inc. Watertown, Mass. 02172. Whitaker, S., AIChE J., 13, 420 (1967a). Whitaker, S.,Ind. Eng. Chem., Fundam., 6, 476 (1967b). Whitaker, S., "Simultaneous Heat, Mass and Momentum Transfer in Porous Media: A Theory of Drying", Advances in Heat Transfer, Vol. 13, Academic Press, New York, N.Y., 1977.
Receiued for reuiew December 13, 1976 Accepted M a y 20,1977
This work was supported by NSF Grant ENG 74-13037 A01.
Impurity Trapping during Dendritic Crystal Growth. 1. Computer Simulation Allan S. Myerson and Donald J. Kirwan' Department of Chemical Engineering, University of Virginia, Charlottesville, Virginia 2290 1
A computer simulation of dendritic crystal growth was developed and used to derive a relationship between sobtion trapping and system and process parameters during crystallization. The relation predicts that trapping will increase with an increase in crystal growth rate and with a decrease in the absolute value of either negative or positive interfacial temperature gradients.
In many instances the distribution of a species between a crystalline solid and the interfacial liquid from which it is growing cannot be described by the equilibrium distribution coefficient, ko, as obtained from a phase diagram. Many studies have observed the existence of liquid occlusions in crystals grown from solution (Belyustin and Fridman, 1968; Brooks et al., 1968; Denbigh and White, 1966; Von Batchelder and Vaughan, 1966; Pilkington and Dunning, 1973; Wilcox, 1968). Again, in melt growth, concentrations exceeding solid solubility limits during progressive freezing and zone refining have been observed Baker and Cahn, 1969; Hellawell, 1965; Sharp and Hellawell, 1970;Wilcox and Zief, 1967; Terwilliger and Dizio, 1970; Cheng and Pigford, 1971; Edie and Kirwan, 1973; Ozum and Kirwan, 1976). I t was first postulated by Wilcox (1963) and Cheng et al. (1967) that such observations might be accounted for by the macroscopic trapping of pockets of solution, that is, by the existence of a second phase within the crystallized solid. As examples of such data Figure 1 shows our data (Myerson, 1977) for the mass fraction of water found in CaSOp2HzO crystals (in excess of the water of hydration) as a function of crystal growth rate, while Figure 2 is a similar plot from Edie and Kirwan (1973) for the system salol-thymol. In both systems the equilibrium amount of the second component should be zero, yet large quantities are found in the grown crystal. Edie and Kirwan (1973) attempted to develop a quantitative description of this trapping process. In their proposed 414
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
mechanism a crystal interface is thought to become unstable and to break up into cellular and dendritic projections under certain conditions. These dendrites subsequently trap pockets of liquid due to the growth of side branches or direct impingment of one dendrite upon another. The conditions under which a growing crystal interface becomes unstable have been well studied (Sekerka, 1968). The simplest, experimentally verified description is the constitutional supercooling criterion (Tiller et al., 1953; Tiller and Rutter; 1956) which states that, for a stable planar interface
The above criterion applies to systems with no solid-solid solubility (ko = 0) as well as solid solution formers. For steady-state crystallization of solid solution formers C d O ) may be replaced with CL(m ) / k o . Thus, under crystallization conditions where the interfacial temperature gradient, velocity, and composition are such as to produce unstable crystal growth we would expect the formation of a dendritic interface. Edie and Kirwan (1973) defined f as the fraction of interfacial solid that is truly crystalline and related to it an apparent interfacial distribution coefficient, k,, defined as the ratio of the apparent concentration of a species in the solid phase to its interfacial liquid concentration.
c, = fCC(0) + (1 - f ) C L ( O )
(2)
