Ind. Eng. Chem. Res. 1999, 38, 833-844
833
Dissolution of an Immobile Phase during Flow in Porous Media Michel Quintard Institut de Me´ canique des Fluides de Toulouse, Av. du Professeur Camille Soula, 31400 Toulouse, France
Stephen Whitaker* Department of Chemical Engineering and Material Science, University of California at Davis, Davis, California 95616
In this paper we examine the problem of the dissolution of an immobile phase during flow in porous media with the objective of deriving an evolutionary transport equation for the volume fraction of the immobile phase. Both the speed of displacement of the phase interface and the interfacial mass-transfer coefficient are represented in terms of mapping variables that are determined by the solution of two closure problems. Other problems of this type include the dissolution and phase change of solids and liquids which can occur because of concentration and temperature gradients. The dissolution of solids is of interest to geologists concerned with transport in sedimentary basins, while engineers are often confronted with freezing and melting processes in multiphase systems. The analysis presented in this paper provides a theoretical framework for the solution of these problems. 1. Introduction The process under consideration is illustrated in Figure 1 where the immobile phase is identified as the γ phase, the rigid solid as the σ phase, and the flowing fluid as the β phase. Upstream from the region that we have shown in Figure 1, the γ-phase is undoubtedly mobile and the complex process of two-phase flow plays a key role in the transport of chemical contaminants.1-7 Downstream from the region under consideration, there is a process of single-phase flow in which dispersion,8 adsorption,9 and biological activity10-12 all influence the contaminant transport process. For the region under consideration in this study, mass transfer from the immobile phase to the mobile phase dominates the process,13-15 and the volume fraction of the immobile phase as a function of time and space is the quantity of primary interest. A direct experimental study of this process for TCE has been carried out by Imhoff et al.;16 however, the interpretation of the data is not completely consistent with our analysis of the mass transport process. The problem of determining the speed of displacement during dissolution is inherently nonlinear; however, for a dilute solution of the diffusing species we show the speed of displacement of the β-γ interface given by
w‚nβγ )
1 n ‚D (I + ∇bβ)‚∇〈cAβ〉β + c°γ βγ Am 1 n ‚D ∇s (ceq - 〈cAβ〉β) c°γ βγ Am β Aβ
Here bβ and sβ are closure variables or mapping variables that must be determined by the solution of two closure problems. This result for w‚nβγ can be used to explore the dissolution process at the point scale, while at the Darcy scale it is the volume fraction of the * Author to whom all correspondence should be addressed.
Figure 1. Dissolution of an immobile phase.
γ phase as a function of time and space that is of interest. The evolutionary transport equation for the volume fraction of the γ phase is derived in section 4, and it takes the form
c°γ
∂γ eq + DAm∇〈cAβ〉β‚∇γ ) -avk(cAβ - 〈cAβ〉β) ∂t
eq - 〈cAβ〉β must be determined by Here ∇〈cAβ〉β and cAβ
10.1021/ie980212t CCC: $18.00 © 1999 American Chemical Society Published on Web 01/22/1999
834
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
solution of the macroscopic mass transport equation, while avk is determined from the closure variable, sβ. Both of these results are limited by ωAβ , 1, which allows for the linearization of the closure problem. Elsewhere17 we have considered the general N-component mass-transfer process in order to develop a useful form of the flux boundary condition at the γ-β interface; however, in this particular study, we restrict our analysis to the simplest possible process in which the γ phase consists of a single component that is slightly soluble in the β phase. This will allow us to direct our attention to the determination of the speed of displacement18 of the β-γ interface and the evolution of the volume fraction of the γ phase in space and time. This special case is described by the following equations and boundary conditions for species A:
∂cAβ + ∇‚(vAβcAβ) ) 0, in the β phase ∂t
(1)
B.C.1
eq , at the β-γ interface cAβ ) cAβ
B.C.2
cAβ(vAβ - w)‚nβγ ) -c°γw‚nβγ, at the β-γ interface (3)
B.C.3
cAβ(vAβ - w)‚nβσ ) 0, at the β-σ interface (4)
Slattery22 and Whitaker23 made use of the traditional volume average that we use in this paper. We begin by defining the superficial average of cAβ according to
(2)
The boundary conditions for all other chemical species are given by
B.C.4
cBβ(vBβ - w)‚nβγ ) 0, B ) 1, 2, ..., N, B * A, at the β-γ interface (5)
B.C.5
cBβ(vBβ - w)‚nβσ ) 0, B ) 1, 2, ..., N, B * A, at the β-σ interface (6)
eq In the first two of these boundary conditions, cAβ represents the β-phase concentration that is in equilibrium with the γ phase and c°γ represents the molar concentration of pure species A in the γ phase. Because the γ phase is assumed to be immobile, we have used the condition
vAγ‚nβγ ) 0, at the β-γ interface
Figure 2. Position vectors associated with the averaging volume.
(7)
in the construction of the flux condition given by eq 3. In this simplified version of the mass transport process, eq 3 becomes an auxiliary condition that allows us to determine the speed of displacement, w‚nβγ. In our previous study19 of this problem we directed our attention to the β-phase transport process to determine the dispersion tensor and the mass-transfer coefficient; however, the γ-phase transport problem was left unresolved.
〈cAβ〉|x )
β
β
dV
(8)
where V represents the averaging volume illustrated in Figure 2 and Vβ(x,t) is the volume of the β phase contained in the averaging volume. This definition clearly indicates that the average concentration is associated with the centroid of the averaging volume, indicated by the position vector x, while integration over Vβ(x,t) is performed in terms of the components of the relative position vector yβ. In general, we will avoid the precise notation used in eq 8 and simply express the superficial average as
〈cAβ〉 )
∫V (t)cAβ dV
1 V
(9)
β
In addition to the superficial average concentration, we will make use of the intrinsic average concentration defined by
〈cAβ〉β )
1 Vβ(t)
∫V (t)cAβ dV
(10)
β
These two concentrations are related by
〈cAβ〉 ) β〈cAβ〉β
(11)
in which β is the volume fraction of the β phase defined as
β ) Vβ(t)/V
(12)
Our objective in this paper is to illustrate how one can determine the volume fraction of the γ phase, γ, as a function of time, and this requires that we determine the derivative of γ, which can be expressed as
2. Volume Averaging Our analysis of the mass transport problem suggested by Figure 1 is based on the development of a volumeaveraged transport equation20-23 for species A. In the original papers dealing with this method, Anderson and Jackson20 and Marle,21 quite independently of each other, made use of a weighting function28 in their definition of the volume average. At the same time,
∫V (x,t)cAβ|x+y
1 V
∂β ∂γ )∂t ∂t
(13)
We begin our analysis with the superficial average of eq 1, which can be expressed as
1 V
∂c
∫V (t) ∂tAβ dV + V1 ∫V (t)∇‚(vAβcAβ) dV ) 0 β
β
(14)
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 835
To develop a governing differential equation for 〈cAβ〉β, we will make use of the general transport theorem
d dt
∂cAβ dV + c dV ) V (t) Vβ(t) Aβ β ∂t
∫
∫
∫A
βγ(t)
cAβw‚nβγ dA (15)
and the spatial averaging theorem given by
∫
1 V
∫V (t)cAβvAβ dV) + 1 ∫ c v ‚n dA V A (t) Aβ Aβ βγ
1 ∇‚(vAβcAβ) dV ) ∇‚ Vβ(t) V
(
β
(16)
Both of these results are special forms based on the assumption that the σ phase is rigid, i.e.
w‚nβσ ) 0, at the β-σ interface
(17)
and that the flux condition given by eq 4 is valid. We can use these theorems to interchange integration and differentiation in eq 14 and obtain a result that is explained in detail elsewhere24
∂ ( 〈c 〉β) + ∇‚〈vAβcAβ〉 + ∂t β Aβ 1 c (v - w)‚nβγ dA ) 0 (18) V Aβγ(t) Aβ Aβ
∫
This represents the governing differential equation for 〈cAβ〉β, and to use this equation to predict the behavior of the system illustrated in Figure 1, we need to be able to determine β, vAβ, and the interfacial flux that we identify as
}
interfacial flux per ) unit volume from the β phase to the γ phase 1 c (v - w)‚nβγ dA (19) V Aβγ(t) Aβ Aβ
∫
Directing our attention to 〈cAβvAβ〉, we impose the condition of a dilute solution of species A to obtain
cAβvAβ ) cAβvβ - DAm∇cAβ
(20)
The derivation of this result is given in the appendix, where we show that the mixture diffusivity for species A takes the form
1 DAm
B)N
)
xBβ
∑D
B)1 B*A
cAβ|x+yβ ) 〈cAβ〉β|x+yβ + c˜ Aβ|x+yβ
(24a)
vβ|x+yβ ) 〈vβ〉β|x+yβ + v˜ β|x+yβ
(24b)
where the vector x + yβ is identified in Figure 2. We can use these results and follow the work of Carbonell and Whitaker26 to obtain
βγ
{
we use Gray’s25 decompositions
(21)
AB
∫
∂ 1 c (v ( 〈c 〉β) + ∇‚(β〈vβ〉β〈cAβ〉β) + ∂t β Aβ V Aβγ(t) Aβ Aβ w)‚nβγ dA ) ∇‚〈DAm∇cAβ〉 - ∇‚〈v˜ βc˜ Aβ〉 (25) Here 〈v˜ βc˜ Aβ〉 represents the dispersive flux for which we require a closure problem that will allow us to express this flux in terms of 〈cAβ〉β. To be precise about the dispersive flux, we note that it is given explicitly by
〈v˜ βc˜ Aβ〉 )
1 V
∫V (t)v˜ β|x+y c˜ Aβ|x+y β
β
β
dV
(26)
For a dilute solution of species A and no interfacial transport of the other species present in the β phase, one can ignore variations in the mixture diffusivity to express the diffusive flux as
[
〈DAm∇cAβ〉 ) DAm〈∇cAβ〉 ) DAm ∇〈cAβ〉 + 1 V
∫A
βγ(t)
1 V
nβγcAβ dA +
∫A
βσ(t)
]
nβσcAβ dA (27)
To eliminate cAβ from this result, we can use the decomposition given by eq 24a along with the representation for the intrinsic average concentration given by eq 11 to obtain
[
〈DAm∇cAβ〉 ) DAm β∇〈cAβ〉β + 〈cAβ〉β∇β + 1 V
∫A (t)nβγ〈cAβ〉β|x+y 1 ∫ n c˜ V A (t) βγ Aβ|x+y
∫A ∫A
1 V 1 dA + β V
β
βγ
βγ
dA +
βσ(t)
|
nβσ〈cAβ〉β x+yb dA +
βσ(t)
|
]
nβσc˜ Aβ x+yβ dA (28)
When the intrinsic average concentration, 〈cAβ〉β, can be removed from the area integrals,25,27,28 one can use a result from the averaging theorem
∇β ) -
∫A
1 V
βγ(t)
nβγ dA -
∫A
1 V
nβσ dA
βσ(t)
(29)
and vβ is the mass-average velocity defined by to express eq 28 in the form
A)N
vβ )
∑ ωAvAβ
(22)
A)1
Use of eq 20 in the superficial average equation represented by eq 18 leads to
∂ ( 〈c 〉β) + ∇‚〈vβcAβ〉 + ∂t β Aβ 1 c (v - w)‚nβγ dA ) ∇‚〈DAm∇cAβ〉 (23) V Aβγ(t) Aβ Aβ
∫
To develop a useful form of the convective transport,
[
〈DAm∇cAβ〉 ) DAm β∇〈cAβ〉β +
∫A (t)nβγc˜ Aβ dA + 1 ∫ n c˜ dA] V A (t) βσ Aβ
1 V
βγ
(30)
βσ
Here it is understood that terms inside the area integrals are evaluated at x + yβ whereas terms outside the area integrals are evaluated at x where x and yβ are the position vectors illustrated in Figure 2. Use of eq 30 in eq 25 leads to the governing differential equation for 〈cAβ〉β given by
836
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
∫A (t)cAβ(vAβ 1 w)‚nβγ dA ) ∇‚[DAm(β∇(cAβ〉β + ∫A (t)nβγc˜ Aβ dA + V 1 ∫ n c˜ dA)] - ∇‚〈v˜ βc˜ Aβ〉 (31) V A (t) βσ Aβ
∂ 1 ( 〈c 〉β) + ∇‚(β〈vβ〉β〈cAβ〉β) + ∂t β Aβ V
βγ
βγ
βσ
In this result, the superficial velocity, 〈vβ〉 ) β〈vβ〉β, will be determined by Darcy’s law29,30 while the diffusive and dispersive terms on the right-hand side will be determined by the closure problem for c˜ Aβ. At this point, the accumulation and convective transport terms are in a conservative form, and for the purpose of subsequent numerical solutions, it seems best to retain this form. However, the nonconservative form is required for the development of the closure problem as we shall see in the next section.
Here we have identified all the terms that are homogeneous in c˜ Aβ, and we have used the term nonlocal to identify those terms in which the spatial deviation concentration is evaluated at points other than the centroid of the averaging volume. Because the nonlocal terms are averaged quantities, the divergence of these terms gives rise to the inverse of the large length scale and this leads to the inequalities
-1 ˜ βc˜ Aβ〉 , ∇‚(vβc˜ Aβ) β ∇‚〈v
[
-1 β ∇‚
DAm V
∫A (t)nβc˜ Aβ dA β
]
(35a)
, ∇‚(DAm∇c˜ Aβ) (35b)
Under these circumstances we can simplify eq 34 to the form
3. Closure To develop a governing differential equation for the spatial deviation concentration, c˜ Aβ, we first note that eq 1 can be expressed as
∂cAβ + ∇‚(vβcAβ) ) ∇‚(DAm∇cAβ) ∂t
(32)
This suggests that we should arrange eq 31 in the form
∂〈cAβ〉β ∂β β β + -1 〈c 〉β + -1 β ∇‚(β〈vβ〉 〈cAβ〉 ) + β ∂t ∂t Aβ -1 β c (v - w)‚nβψ dA ) ∇‚(DAm∇〈cAβ〉β) + V Aβγ(t) Aβ Aβ DAm β -1 -1 n c˜ dA + β ∇β‚DAm∇〈cAβ〉 + β ∇‚ V Aβγ(t) βγ Aβ DAm n c˜ dA - -1 ˜ βc˜ Aβ〉 (33) β ∇‚〈v V Aβσ(t) βσ Aβ
∫
[
∫
]
∫
and subtract this result from eq 32 to obtain (after some algebra)
in which the nonhomogeneous terms have been identified as sources of the c˜ Aβ field. This form of the transport equation for c˜ Aβ is more complex than that in the case of passive transport26 because of the interfacial flux and because of the time and space dependence of β. 3.1. Total Continuity Equation. To develop a more attractive form for the convective source terms, we need to make use of the total continuity equation. Without providing any constraints,31 we assume that the convective source terms in eq 36 can be successfully simplified by use of the incompressible form of the total continuity equation
∇‚vβ ) 0
(37)
This form of the total continuity equations is routinely used in the solution of liquid-phase flow problems; however, it is in the nature of an approximation to make use of this result to simplify the species A transport equation represented by eq 36. The key idea here is that certain terms have been discarded to arrive at eq 37, and while these terms may be small in the sense that they do not affect the fluid mechanical problem, they may not be small in terms of their impact on eq 36. Nevertheless, we will proceed with the use of the traditional form of the continuity equation for incompressible flow. The volume-averaged form of eq 37 is analogous to that of eq 18, and we list the
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 837
result as
the form
∂β 1 + ∇‚〈vβ〉 + ∂t V
∫A
(vβ - w)‚nβγ dA ) 0 (38)
βγ(t)
On the basis of eq 5 and the definition of the massaverage velocity given by eq 22, this result can be expressed as
∂β 1 + ∇‚〈vβ〉 + ∂t V
∫A
ωAβ(vAβ - w)‚nβγ dA ) 0
βγ(t)
(39)
This form of the volume-averaged continuity equation differs from that used by Quintard and Whitaker,19 who neglected the third term with the idea that this would have little influence on the divergence of 〈vβ〉. While it is true that this term will have a negligible effect on the fluid mechanical problem, it is prudent to retain the complete form of eq 39 to determine its impact on the species A mass transport process. Keeping in mind that 〈vβ〉 ) β〈vβ〉β, we can use eq 39 to construct the following result: -1 β β β β -1 β 〈cAβ〉 〈vβ〉 ‚∇β ) -〈cAβ〉 ∇‚〈vβ〉 - β β -1 β 〈cAβ〉 V
∫A
∂β 〈c 〉β ∂t Aβ
ωAβ(vAβ - w)‚nβγ dA (40)
βγ(t)
Substitution of this form of the total continuity equation into eq 36 provides a simplified version of the transport equation for c˜ Aβ given by
In the appendix we show that the interfacial flux can be expressed in a purely diffusive form that leads to
∫A
1 V
βγ(t)
cAβ(vAβ - w)‚nβγ dA ) -
1 V
∫A
βγ(t)
DAm∇cAβ‚nβγ dA (45)
and use of the decomposition given by eq 24a provides
∫A
1 V
βγ(t)
cAβ(vAβ - w)‚nβγ dA ) nβγ dA -
1 V
∫A
∫A
1 V
βγ(t)
DAm∇〈cAβ〉β‚
βγ(t)
DAm∇c˜ Aβ‚nβγ dA (46)
At this point it is convenient to make use of the theorem given by eq 29 to simplify our representation for the interfacial flux. To do so, we first recognize that the mass flux at the β-σ interface is zero
∫A
1 V
βσ(t)
cAβ(vAβ - w)‚nβσ dA ) 0
(47)
and we then use the diffusive representation illustrated by eq 46 to obtain
0)-
∫A
1 V
DAm∇〈cAβ〉β‚nβσ dA -
βσ(t)
1 V
For a dilute solution of species A, we have the condition
ωAβ , 1
(42)
and this leads to the following restriction associated with the interfacial flux:
〈cAβ〉
V
∫A
βγ(t)
ωAβ(vAβ - w)‚nβγ dA , -1 β V
∫A
βγ(t)
cAβ(vAβ - w)‚nβγ dA (43)
Under these circumstances the governing differential equation for the spatial deviation concentration takes
βσ(t)
DAm∇c˜ Aβ‚nβσ dA (48)
This allows us to express the interfacial flux in terms of area integrals over both the β-γ interface and the β-σ interface so that the total flux takes the form
∫A
1 V
βγ(t)
-
cAβ(vAβ - w)‚nβγ dA )
∫A
1 V
βγ(t)
DAm∇〈cAβ〉β‚nβγ dA -
nβγ dA -
-1 β β
∫A
∫A
1 V
βσ(t)
∫A
1 V
DAm∇c˜ Aβ‚
βγ(t)
DAm∇〈cAβ〉β‚nβσ dA 1 V
∫A
βσ(t)
DAm∇c˜ Aβ‚nβσ dA (49)
The disadvantage of two additional area integrals in the expression for the interfacial flux becomes an advantage when we remove DAm∇〈cAβ〉β from those integrals. This can be done when certain length-scale constraints are
838
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
valid,27,28 and when it is done, we can use eq 29 to express eq 49 as
∫ ∫
1 c (v - w)‚nβγ dA ) V Aβγ(t) Aβ Aβ 1 1 DAm∇c˜ Aβ‚nβγ dA D ∇c˜ ‚n dA + A (t) βγ V V Aβσ(t) Am Aβ βσ ∇β‚DAm∇〈cAβ〉β (50)
∫
Substitution of this representation for the volumeaveraged interfacial flux into eq 41 leads to Figure 3. Length scales associated with the averaging volume.
transport equation for c˜ Aβ to obtain
In this form of the governing differential equation for c˜ Aβ, we find accumulation, convection, diffusion, and interfacial exchange for the c˜ Aβ field. In addition, we find convective and diffusive sources that contribute to the generation of the c˜ Aβ field. 3.2. Boundary Conditions. At this point we are ready to move on to the construction of the boundary conditions that are necessary to complete our statement of the closure problem. Use of the spatial decomposition for the concentration of species A with the boundary condition given by eq 2 leads to
While the macroscopic transport process described by eq 31 is generally transient, the closure problem is normally quasi-steady on the basis of the constraint
DAmt*/lβ2 . 1
(56)
Here t* is the appropriate characteristic process time and lβ is the small length scale associated with the β phase that is illustrated in Figure 3. When the constraint indicated by eq 56 is valid, the closure problem can be expressed as
Here we have identified the difference between the equilibrium concentration and the intrinsic average concentration as a source of the c˜ Aβ field. The second boundary condition that is required for our solution of the closure problem is given by eq 4, and in the appendix we show that this condition takes the form
B.C.3
-DAm∇cAβ‚nβσ ) 0, at the β-σ interface (53)
When the decomposition given by eq 24a is used with this result, we obtain a boundary condition for the c˜ Aβ field that contains another source.
A little thought32,33 will indicate that the diffusive source in the boundary condition given by eq 54 makes a contribution to c˜ Aβ that is much, much larger than the diffusive source in the transport equation given by eq 51. Under these circumstances, we can simplify the
Here we have made use of spatially periodic boundary conditions34,35 for c˜ Aβ with the idea that there is no need
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 839
Here the area per unit volume times the mass-transfer coefficient is given by
Figure 4. Representative region of the three-phase system.
to solve the closure problem in the macroscopic region illustrated in Figure 1. Instead, we require only a solution in some representative region, such as the region illustrated in Figure 4, and in that region we have assumed that variations in the mixture diffusivity for species A can be ignored. The boundary value problem given by eq 57 contains only two nonhomogeeq - 〈cAβ〉β, and a little thought neous terms, ∇〈cAβ〉 and cAβ will indicate that the solution for c˜ Aβ will be a constant when the two nonhomogeneous terms are zero. This constant will not pass through the filters represented by the area and volume integrals in eq 31; thus, it is eq - 〈cAβ〉β, that produce the source terms, ∇〈cAβ〉 and cAβ the nontrivial part of c˜ Aβ. This suggests a solution of the form β c˜ Aβ ) bβ‚∇〈cAβ〉β + sβ(ceq A - 〈cAβ〉 )
(58)
in which bβ and sβ are referred to as the closure variables. Representations of this form are discussed in detail elsewhere,33,36 and use of this result in eq 57 leads to two boundary value problems for the two closure variables. The closure problem for the vector closure variable is given by
Problem I v˜ β + ∇‚(vβbβ) ) DAm∇2bβ - -1 β uβ
(59a)
B.C.1
bβ ) 0, at Aβγ(t)
(59b)
B.C.2
nβσ‚∇bβ ) -nβσ, at Aβσ(t)
(59c)
Periodicity:
bβ(r+li) ) bβ(r), i ) 1, 2, 3 〈bβ〉β ) 0
Average:
(59d)
∫A
1 V
βγ(t)
nβγ‚DAm∇bβ dA +
1 V
)
1 V
∫A
βγ(t)
∫A
nβγ‚DAm∇sβ dA +
∫A
1 V
nβσ‚DAm∇sβ dA
βσ(t)
nβγ‚DAm∇sβ dA
(62)
The boundary value problems for bβ and sβ can be solved by well-known numerical methods,19 and the results can be used to determine the effective coefficients that appear in the closed form of eq 31. 3.3. Closed Form of the Macroscopic Equation. We begin our development of the closed form of the macroscopic equation by using eq 50 to express eq 31 as
∂ ( 〈c 〉β) + ∇‚(〈vβ〉〈cAβ〉β) ∂t β Aβ 1 1 D ∇c˜ ‚n dA D ∇c˜ ‚n dA ) V Aβγ(t) Am Aβ βγ V Aβσ(t) Am Aβ βσ 1 n c˜ dA + ∇‚ DAm β∇〈cAβ〉β + V Aβγ(t) βγ Aβ 1 n c˜ dA - ∇β‚DAm∇〈cAβ〉β - ∇‚〈v˜ βc˜ Aβ〉 (63) V Aβσ(t) βσ Aβ
∫
∫ ∫
[ (
∫
)]
Here we have replaced β〈vβ〉β with the superficial velocity 〈vβ〉 since this will be available to us through the use of Darcy’s law. At this point we are ready to make use of the representation given by eq 58 in order to obtain the closed form of eq 63 that can be expressed as
∂ eq ( 〈c 〉β) + ∇‚(〈vβ〉〈cAβ〉β) - ∇‚[dβ(cAβ - 〈cAβ〉β)] ∂t β Aβ uβ‚∇〈cAβ〉β ) ∇‚(D/β‚∇〈cAβ〉β) - ∇β‚DAm∇〈cAβ〉β + eq - 〈cAβ〉β) (64) avk(cAβ
The dispersion tensor is defined by
(
D/β ) DAm βI +
1 V
∫A
βγ(t)
nβγbβ dA +
∫A
1 V
)
nβσbβ dA - 〈v˜ βbβ〉 (65)
βσ(t)
and the additional velocity-like coefficient is given by
nβσ‚DAm∇bβ dA (60)
βσ(t)
dβ ) DAm
(V1 ∫
nβγsβ dA +
Aβγ(t)
∫A
1 V
nβσsβ dA
βσ(t)
)
(66)
(61a)
Calculated values of the coefficients that appear in eq 64 have been presented by Quintard and Whitaker;19 however, the coefficients are sensitive to the manner in which the three phases are distributed, and much work remains to be done in terms of identifying the parameters that can be used to characterize three-phase systems.
For the scalar closure variable in the representation given by eq 58, we find the following boundary value problem:
Problem II ∇‚(vβsβ) ) DAm∇2sβ - -1 β avk
∫A
1 V
βγ(t)
(59e)
in which the constant vector uβ is given by
uβ )
avk )
B.C.1
sβ ) 1, at Aβγ(t)
(61b)
4. Volume Fraction of the Immobile Phase
B.C.2
nβσ‚∇sβ ) 0, at Aβσ(t)
(61c)
At this point we wish to make use of the auxiliary condition given by eq 3 to first determine the speed of displacement and then the volumetric rate of dissolution of the γ phase. We first arrange the auxiliary condition in the form
Periodicity:
sβ(r+li) ) sβ(r), i ) 1, 2, 3
Average: 〈sβ〉β ) 0
(61d) (61e)
840
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
B.C.2
w‚nβγ ) -
1 [c (v - w)‚nβγ], c°γ Aβ Aβ at the β-γ interface (67)
and then make use of the dilute solution approximation given in the appendix to obtain
To be clear about this point, we examine the diffusive term in eq 73 to obtain
∫A (t)nβγ‚DAm(I + ∇bβ) dA ) V1 ∫A (t)nβγDAm dA + 1 ∫ D n ‚∇bβ dA + V1 ∫A (t)DAmnβσ‚∇bβ dA V A (t) Am βγ 1 ∫ D n ‚∇bβ dA (75) V A (t) Am βσ
1 V
βγ
βγ
βγ
1 w‚nβγ ) nβγ‚DAm∇cAβ c°γ
(68)
βσ
The concentration can be decomposed according to eq 24a
w‚nβγ )
1 1 n ‚D ∇〈c 〉β + nβγ‚DAm∇c˜ Aβ c°γ βγ Am Aβ c°γ
(69)
On the basis of the boundary condition given by eq 59c, this result can be arranged in the form
∫A
1 V
1 n ‚D (I + ∇bβ)‚∇〈cAβ〉β + c°γ βγ Am 1 n ‚D ∇s (ceq - 〈cAβ〉β) (70) c°γ βγ Am β Aβ
This result can be used to explore the details of the speed of displacement of the β-γ interface. For example, if one wants to know how the fluid-fluid interface illustrated in Figure 4 evolves with time, one need only solve the closure problems given by eqs 59 and 61 and eq supply the parameters, DAm, c°γ, ∇〈cAβ〉β, and cAβ β 〈cAβ〉 , to determine the speed of displacement, w‚nβγ. To determine how the volume fraction of the γ phase evolves with time, we first make use of the general transport theorem to obtain the purely kinematic result
∂γ 1 ∂β )) ∂t ∂t V
∫A
βγ(t)
w‚nβγ dA
(71)
}
∂γ 1 1 )n ‚D (I + ∇bβ) dA ‚∇〈cAβ〉β ∂t c°γ V Aβγ(t) βγ Am 1 1 eq n ‚D ∇s dA (cAβ - 〈cAβ〉β) (72) c°γ V Aβγ(t) βγ Am β
∫
{
}
∫
nβγ‚DAm(I + ∇bβ) dA )
∫A
βσ(t)
DAmnβσ dA +
∫A
1 V
βγ(t)
{
∂γ 1 )∂t V
∫A
βγ(t)
}
eq - 〈cAβ〉β) (73) avk(cAβ
∂γ 1 ) ∂t V
∫A
cAβ(vAβ - w)‚nβγ dA
βγ(t)
DAmnβγ dA +
∫ ∫
1 D n dA + V Aβγ(t) Am βγ 1 1 D n dA + D ∇2bβ dV (77) V Aβσ(t) Am βσ V Vγ(t) Am
∫
On the basis of the governing differential equation for bβ and the geometrical relation given by eq 29, this result takes the form
∫A
1 V
βγ(t)
nβγ‚DAm(I + ∇bβ) dA )
∫V (t)[v˜ β + ∇‚(vβbβ)] dV
1 V
(78)
β
Making use of 〈v˜ β〉 ) 0, periodicity for bβ, and the divergence theorem leads to
∫A
1 V
βγ(t)
nβγ‚DAm(I + ∇bβ) dA ) -DAm∇β +
∫A
1 V
nβγ‚(vβbβ) dA +
βγ(t)
∫A
1 V
βσ(t)
nβσ‚(vβbβ) dA (79)
We can now use the condition
nβσ‚vβ ) 0, at Aβσ(t)
(80)
along with the boundary condition given by eq 59b to conclude that
nβγ‚DAm(I + ∇bβ) dA ‚∇〈cAβ〉β -
Typically, one thinks of the evolution of the volume fraction of the γ phase as being determined by the mass flux from the γ phase to the β phase, and this is certainly the case because one can use eqs 67 and 71 to express the rate of change of γ as
c°γ
βγ(t)
1 D n ‚∇bβ dA + V Aβγ(t) Am βγ 1 D n ‚∇bβ dA (76) V Aβσ(t) Am βσ
nβγ‚DAm(I + ∇bβ) dA )
and on the basis of the definition of the mass-transfer coefficient given by eq 62, this leads to
c°γ
∫ ∫
∫A
1 V
The periodicity condition given by eq 59d allows us to use the divergence theorem to simplify the last two integrals, leading to
-DAm∇β +
Substitution of eq 70 into this result leads to
{
βγ(t)
1 V
and one can use the representation for the spatial deviation concentration given by eq 58 to express the speed of displacement as
w‚nβγ )
βσ
(74)
The fact that the evolution of γ is not determined entirely by the mass-transfer coefficient in eq 73 results eq - 〈cAβ〉β) does not represent from the fact that avk(cAβ the total mass flux between the two phases.
∫A
1 V
nβγ‚DAm(I + ∇bβ) dA ) -DAm∇β
βγ(t)
(81)
This result indicates that the rate of change of the volume fraction of the γ phase is determined by
c°γ
∂γ eq ) DAm∇β‚∇〈cAβ〉β - avk(cAβ - 〈cAβ〉β) (82) ∂t
At this point we can make use of the relation
∇β ) -∇γ
(83)
to express the evolutionary equation for the volume fraction of the γ phase according to
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 841
c°γ
∂γ eq - 〈cAβ〉β) + DAm∇〈cAβ〉β‚∇γ ) -avk(cAβ ∂t
(84)
∇‚(D/β‚∇〈cAβ〉β) . ∇‚(〈vβ〉〈cAβ〉β) ∇‚(D/β‚∇〈cAβ〉β) . ∇‚(〈dβ〉〈cAβ〉β)
To determine the influence of the diffusive flux, DAm∇〈cAβ〉β, and the exchange flux on the right-hand side of this result, we consider the two extreme cases of convection-controlled mass transfer and diffusioncontrolled mass transfer. In both cases, we will examine processes in which the concentration gradient, ∇ 〈cAβ〉β, is generated by mass exchange between the γ phase and the β phase. 4.1. Dominant Convective Transport. In this case we can impose the restrictions given by eq - 〈cAβ〉β)] ∇‚(〈vβ〉〈cAβ〉β) . ∇‚[dβ(cAβ
∇‚(D/β‚∇〈cAβ〉β) . uβ‚∇〈cAβ〉β Under these circumstances, a reasonable representation of the macroscopic equation is given by
∂ eq ( 〈c 〉β) ≈ ∇‚(D/β‚∇〈cAβ〉β) + avk(cAβ - 〈cAβ〉β) ∂t β Aβ
[
(85)
∇‚(〈vβ〉〈cAβ〉β) . ∇β‚DAm∇〈cAβ〉β
∂ eq ( 〈c 〉β) + ∇‚(〈vβ〉〈cAβ〉β) ) avk(cAβ - 〈cAβ〉β) ∂t β Aβ
(86)
]
(87)
]
(88)
in which the Pe´clet number is loosely defined by
Pe ) 〈vβ〉lβ/DAm
(89)
For this condition, it is clear that the terms in eq 84 are related by eq - 〈cAβ〉β) DAm∇〈cAβ〉β‚∇γ , avk(cAβ
(90)
and the evolutionary equation for the volume fraction of the γ phase takes the form
c°γ
∂γ eq ) -avk(cAβ - 〈cAβ〉β) ∂t
L∇γ ) O(1)
(96)
eq - 〈cAβ〉β)] DAm∇〈cAβ〉β‚∇γ ) O[avk(cAβ
(97)
and the determination of the rate of change of the volume fraction of the γ phase given by eq 84 requires eq a knowledge of both ∇〈cAβ〉β and cAβ - 〈cAβ〉β.
This leads to the estimate
lβ∇γ a k(ceq - 〈cAβ〉β) Pe v Aβ
(95)
Under these circumstances, we have
1 a k(ceq - 〈cAβ〉β) 〈vβ〉 v Aβ
[
(94)
and for a diffusion-controlled process, a reasonable estimate of the gradient of the volume fraction of the γ phase is given by
If all three terms in this result are of comparable magnitude, a reasonable estimate of the gradient of the average concentration is given by
DAm∇γ‚∇〈cAβ〉β ) O
]
eq - 〈cAβ〉β)] DAm∇〈cAβ〉β‚∇γ ) O[(L∇γ)avk(cAβ
so that eq 64 simplifies to
[
L a k(ceq - 〈cAβ〉β) DAm v Aβ
Here we have used L as the characteristic length associated with significant changes in ∇〈cAβ〉β, and we have used DAm to approximate D/β. This estimate of ∇〈cAβ〉β leads to
∇‚(〈vβ〉〈cAβ〉β) . ∇‚(D/β‚∇〈cAβ〉β)
∇〈cAβ〉β ) O
(93)
and our estimate of the concentration gradient takes the form
∇〈cAβ〉β ) O
∇‚(〈vβ〉〈cAβ〉β) . uβ‚∇〈cAβ〉β
(92)
(91)
This is the basis for the traditional interpretation of a mass-transfer coefficient. It is important to remember that this result is associated with Pe´clet numbers that are large compared to 1, and when this is not the case, one needs to consider the possibility that a significant fraction of the interfacial mass flux is not proportional eq - 〈cAβ〉β. to the concentration difference, cAβ 4.2. Dominant Diffusive Transport. In this case, we impose the following restrictions on eq 64:
5. Conclusions In this work we have examined the process of the dissolution of an immobile phase caused by mass transfer to a second phase. The analysis has been limited to the simplest possible case of the dissolution of a pure phase into a dilute solution of a mobile phase. The results clearly indicate that the speed of displacement of the interface between the two phases can be determined by solution of the mass-transfer closure problem. The evolution of the volume fraction of the immobile phase is given in terms of a transport equation for γ. This transport equation contains a mass-transfer coefficient, avk, that must be determined by means of the closure problem, in addition to a concentration eq gradient, ∇〈cAβ〉β, and a concentration driving force, cAβ β - 〈cAβ〉 , that must be determined by a solution of the macroscopic mass transport equation. Nomenclature av ) Aβγ(t)/V, interfacial area per unit volume, m-1 Aβγ(t) ) area of the β-γ interface contained in the averaging volume, m2 Aβσ(t) ) area of the β-σ interface contained in the averaging volume, m2 bβ ) closure variable that maps ∇〈cAβ〉β onto c˜ Aβ, m cAβ ) molar concentration of species A in the β phase, mol/ m3
842
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
〈cAβ〉β ) intrinsic average concentration of species A in the β phase, mol/m3 〈cAβ〉 ) β〈cAβ〉β, superficial concentration of species A in the β phase, mol/m3 c˜ Aβ ) cAβ - 〈cAβ〉β, spatial deviation concentration of species A, mol/m3 c°γ ) molar concentration of pure species A in the γ phase, mol/m3 eq cAβ ) concentration of species A in equilibrium with c°γ, mol/m3 cβ ) total molar concentration of the β phase, mol/m3 dβ ) a velocity-like effective transport coefficient, m/s DAB ) binary diffusion coefficient for species A and B, m2/s DAm ) mixture molecular diffusivity for species A, m2/s D/β ) β-phase dispersion tensor, m2/s I ) unit tensor / / JAβ ) cAβuAβ , molar diffusive flux of species A in the β phase, mol/m2s k ) mass-transfer coefficient, m/s lβ ) small length scale associated with the β phase, m li ) 1, 2, 3, lattice vectors for a spatially periodic model of a porous medium L ) large-scale characteristic length, m nβR ) -nRβ, unit normal vector directed from the β phase toward the R phase (R ) σ, γ) Pe ) 〈vβ〉lβ/DAm, Pe´clet number r ) position vector, m eq sβ ) closure variable that maps cAβ - 〈cAβ〉β onto c˜ Aβ t ) time, s t* ) characteristic process time, s uβ ) a velocity-like effective transport coefficient, m/s / uAβ ) vAβ - v/β, molar diffusion velocity of species A in the β phase, m/s V ) averaging volume, m3 Vβ(t) ) volume of the β phase contained within the averaging volume, m3 vAβ ) velocity of species A in the β phase, m/s vβ ) mass-average velocity in the β phase, m/s v/β ) molar-average velocity in the β phase, m/s 〈vβ〉β ) intrinsic mass-average velocity in the β phase, m/s v˜ β ) vβ - 〈vβ〉β, spatial deviation velocity in the β phase, m/s 〈vβ〉 ) β〈vβ〉β, superficial mass-average velocity in the β phase, m/s w‚nβγ ) speed of displacement of the β-γ interface, m/s xAβ ) mole fraction of species A in the β phase x ) position vector locating the centroid of the averaging volume, m yβ ) position vector locating points in the β phase relative to the centroid of the averaging volume, m Greek Letters R ) volume fraction of the R phase, (R ) β, σ, γ) ωAβ ) mass fraction of species A in the β phase.
B)N
v/β )
∑ xBβvBβ
(A2)
B)1
and we need to relate the molar-average velocity to the mass-average velocity since a governing differential equation exists for the latter but not the former. Multiplication of eq A1 by the mass fraction and summing over all species leads to A)N
∑ A)1
A)N
ωAβvAβ )
∑ A)1
A)N
ωAβv/β +
/ ωAβuAβ ∑ A)1
(A3)
This obviously takes the form A)N
vβ ) v/β +
/ ωAβuAβ ∑ A)1
(A4)
and we can express the molar flux of species A as
To arrange the diffusive flux in a convenient form, we / in the following manner manipulate ωAβuAβ / ωBβuBβ )
)
)
/ FBβuBβ FAβ + FBβ + FC + ... + FNβ / ) MB(cBβuBβ MAcAβ + MBcBβ + ... + MNcNβ / MB(cBβuBβ )
cβ(xAβMA + xBβMB + ... xNβMN)
(A6)
Here we have used cβ to represent the total molar concentration of the β phase, and if we define the mean molecular weight as
M h ) xAβMA + xBβMB + ... + xNβMN
(A7)
we can express eq A6 in compact form according to / ωBβuBβ )
/ MB(cBβuBβ )
cβM h
(A8)
This representation allows us to express eq A5 in the form
Appendix B)N
A1. Dilute Solution Representations. Given the flux of species A in the β phase, cAβvAβ, we would like to extract a convective part determined by the massaverage velocity and a diffusive part represented in terms of the molar concentration. We begin by decomposing the species velocity in terms of the molar-average velocity and the molar diffusion velocity according to
/ - xAβ cAβvAβ ) cAβvβ + cAβuAβ
∑ B)1
MB / ) (cBβuBβ M h
(A9)
From eqs A1 and A2 we can obtain B)N
/ ))0 ∑ (cBβuBβ
(A10)
B)1 / vAβ ) v/β + uAβ
The molar-average velocity is defined by
(A1)
and this suggests that all of the diffusive fluxes are the same order of magnitude. If species A is the dominant diffusing species, this means
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 843 / / cBβuBβ ) O(cAβuAβ ), B ) 1, 2, ..., N
(A11)
This estimate does not restrict the diffusive flux of any / ; it species from being much, much less than cAβuAβ / simply means that cBβuBβ cannot be much larger than the diffusive flux of species A. Given these constraints on the diffusive fluxes, the inequality given by B)N
xAβ
∑ B)1
MB / / ) , cAβuAβ (cBβuBβ M h
xAβ , 1
(A13)
(A14)
/ cAβuAβ
)
(A15)
To develop a useful expression for J/A, we begin with the Stefan-Maxwell equation for species A
0 ) -∇xAβ +
B)Nx
∑
AβxBβ(vBβ
- vAβ)
DAB
B)1 B*A
(A16)
and use the decomposition given by eq A1 to obtain
0 ) -∇xAβ +
B)N x x (u/ Aβ Bβ Bβ
∑
/ - uAβ )
DAB
B)1 B*A
B*A
/ ) -cβDAm∇xAβ JAβ
(A17)
/ ) -DAm∇cAβ JAβ
∑
cβDAB
B)1 B*A
cAβvAβ ) cAβvβ - DAm∇cAβ
0 ) -∇xAβ +
/ xAβB)N JBβ
cβ
∑
-
1 cβ
B)1 DAB B*A
{ } B)N
xBβ
∑D
B)1 B*A
/ JAβ
0 ) -∇xAβ +
B)Nx
∑ B)1
AβcBβ(vBβ
- w)
xAβ
J/B
∑D
B)1 B*A
B)N
,
AB
B*A
so that eq A19 simplifies to
0 ) -∇xAβ -
xBβ
∑ B)1 D
1
/ JAβ
0 ) -∇xAβ +
B)Nx
∑ B)1
AβcBβ(vBβ
B)Nx
∑ B)1
{ }
cβ
xBβ
∑ B)1 D
B*A
/ JAβ
(A27)
-
BβcAβ(vAβ
- w)‚nβγ
cβDAB
(A28)
For the particular system under consideration in this work, we have a constraint on the fluxes given by
cBβ(vBβ - w)‚nβγ ) 0, B ) 1, 2, ..., N, B * A (A29) and this allows us to express eq A28 in the form
{ } B)N
B)N
cβDAB
- w)‚nβγ
B*A
AB
- w)
cβDAB
(A19)
(A20)
BβcAβ(vAβ
We now make use this result at the β-γ interface and form the scalar product with nβγ to obtain
AB
{ }{ } B)N
-
cβDAB
B*A
On the basis of the ideas leading from eq A9 to eq A14, we impose the restriction
(A26)
It is important to remember that this representation for the flux is based on the constraint given by eq A13. A1.1. Molar Flux at a Moving Interface. To develop a useful representation for the molar flux at a moving interface, cAβ(vAβ - w)‚nβγ, we return to the Stefan-Maxwell equation for species A and add and subtract the velocity w to obtain
B*A
and the sum can be written as two separate terms to obtain
(A25)
and this allows us to express the molar flux of species A according to
∑ B)1
(A18)
(A24)
Under these conditions the diffusive flux takes the form
B*A
/ - xBβJAβ
(A23)
Given the dilute solution constraint for the diffusing species indicated by eq A13, one can think of cases in which the temperature and pressure gradients are small enough so that the following restriction is satisfied.
B)Nx
B)Nx J/ Aβ Bβ
(A22)
AB
so that the molar diffusive flux can be expressed as
This obviously can be expressed as
0 ) -∇xAβ +
xBβ
∑ B)1 D
xAβ∇cβ , cβ∇xAβ
where we have represented the molar diffusive flux of species A as / JAβ
)
DAm
This allows us to simplify eq A9 to the form / cAβvAβ ) cAβvβ + JAβ
B)N
1
(A12)
is plausible whenever the following constraint is valid:
Constraint:
according to
0 ) -cβ∇xAβ + (A21)
AB
We now define the mixture diffusivity for species A
∑
xBβ
B)1 DAB B*A
cAβ(vAβ - w)‚nβγ
(A30)
At this point we can use the definition of the mixture diffusivity given by eq A22 and the restriction given by eq A24 to express eq A30 as
844
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
cAβ(vAβ - w)‚nβγ ) -nβγ‚DAm∇cAβ
(A31)
For the general case of multicomponent mass transfer, the simplification represented by eq A29 is not available and one must follow the development given by Quintard and Whitaker.17 Literature Cited (1) Abriola, L.; Pinder, G. F. A multiphase approach to the modeling of porous media contaminated by organic compounds. 1. Equation development. Water Resour. Res. 1985, 21, 19-26. (2) Corapcioglu, M. Y.; Baehr, A. L. A compositional multiphase model for groundwater contamination by petroleum products. I. theoretical considerations. Water Resour. Res. 1987, 23, 191-200. (3) Kaluarachchi, J. J. and Parker, J. C. Modeling multicomponent organic chemical transport in three-fluid-phase porous media. J. Contam. Hydrol. 1990, 5, 349-374. (4) Powers, S. E.; Loureiro, C. O.; Abiola, L. M.; Weber, W. J. Theoretical study of the significance of nonequilibrium dissolution of nonaqueous phase liquids in subsurface systems. Water Resour. Res. 1991, 27, 463-477. (5) Robinson, R. L.; Slattery, J. C. Estimation of three-phase relative permeabilities. Transp. Porous Media 1994, 16, 263-287. (6) Ahmadi, A.; Quintard, M. Large-scale properties for twophase flow in random porous media. J. Hydrol. 1996, 183, 6999. (7) Hassanizadeh, S. M.; Gray, W. G. Recent advances in theories of two-phase flow in porous media. In Flow in Porous Media; Du Plessis, J. P., Ed.; Computational Mechanics Publications: Southampton, U.K., 1997; Chapter 3. (8) Plumb, O. A.; Whitaker, S. Diffusion, adsorption, and dispersion in porous media: Small-scale averaging and local volume averaging. In Dynamics of Fluids in Hierarchical Porous Media; Cushman, J. H., Academic Press: New York, 1990; Chapter V. (9) Whitaker, S. Volume averaging of transport equations. In Flow in Porous Media; Du Plessis, J. P., Ed.; Computational Mechanics Publications: Southampton, U.K., 1997; Chapter 1. (10) Baveye, P.; Valocchi, A. An evaluation of mathematical models of the transport of biologically reacting solutes in saturated soils and aquifers. Water Resour. Res. 1989, 25, 1413-1421. (11) Kim, S.; Corapcioglu, M. Y. The role of biofilm growth in bacteria-facilitated contaminant transport in porous media. Transp. Porous Media 1997, 26, 161-181. (12) Wood, B. D.; Whitaker, S. Diffusion and reaction in biofilms. Chem. Eng. Sci. 1998, 53, 397-425. (13) Miller, C. T.; Poirier-McNeill, M. M.; Mayer, A. S. Dissolution of trapped nonaqueous phase liquids: Mass transfer characteristics. Water Resour. Res. 1990, 26, 2783-2796. (14) Powers, S. E.; Abriola, L. M.; Weber, W. J. An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: Steady-state mass transfer rates. Water Resour. Res. 1992, 28, 2691-2705. (15) Mayer, A. S.; Miller, C. T. The influence of mass transfer characteristics and porous media heterogeneity on nonaqueous phase dissolution. Water Resour. Res. 1996, 32, 1551-1567. (16) Imhoff, P. T.; Jaffe´, P. R.; Pinder, G. F. An experimental study of complete dissolution of a nonaqueous phase liquid in saturated porous media. Water Resour. Res. 1993, 30, 307-320.
(17) Quintard, M.; Whitaker, S. The mass flux boundary condition at a moving fluid-fluid interface. Ind. Chem. Eng. Res. 1995, 34, 3508-3513. (18) Slattery, J. C. Interfacial Transport Phenomena; SpringerVerlag: New York, 1990. (19) Quintard, M.; Whitaker, S. Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media. Adv. Water Resour. 1994, 17, 221-239. (20) Anderson, T. B.; Jackson, R. A fluid mechanical description of fluidized beds. Ind. Eng. Chem. Fundam. 1967, 6, 527-538. (21) Marle, C. M. Single-phase flow in porous media. Rev. Inst. Fr. Pet. 1967, 22 (10), 1471-1509. (22) Slattery, J. C. Flow of viscoelastic fluids through porous media. AIChE J. 1967, 13, 1066-1071. (23) Whitaker, S. Diffusion and dispersion in porous media. AIChE J. 1967, 13, 420-427. (24) Whitaker, S. A simple geometrical derivation of the spatial averaging theorem. Chem. Eng. Educ. 1985, 19, 18-21, 50-52. (25) Gray, W. G. A derivation of the equations for multiphase transport. Chem. Eng. Sci. 1975, 30, 229-233. (26) Carbonell, R. G.; Whitaker, S. Dispersion in pulsed systems. II. Theoretical developments for passive disperison in porous media. Chem. Eng. Sci. 1983, 38, 1795-1802. (27) Carbonell, R. G.; Whitaker, S. Heat and mass transfer in porous media. In Fundamentals of Transport Phenomena in Porous Media; Bear, J., Corapcioglu, M. Y., Eds.; Martinus Nijhoff Publishers: Dordrecht, The Netherlands, 1984; pp 123-198. (28) Quintard, M.; Whitaker, S. Transport in ordered and disordered porous media. II. Generalized volume averaging. Transp. Porous Media 1994, 14, 179-206. (29) Whitaker, S. Flow in porous media. II. The governing equations for immiscible two-phase flow. Transp. Porous Media 1986, 1, 105-125. (30) Lasseux, D.; Quintard, M.; Whitaker, S. Determination of permeability tensors for two-phase flow in homogeneous porous media: Theory. Transp. Porous Media 1996, 24, 107-137. (31) Whitaker, S. Levels of simplification: The use of assumptions, restrictions, and constraints in engineering analysis. Chem. Eng. Ed. 1988, 22, 104-108. (32) Ryan, D.; Carbonell, R. G.; Whitaker, S. A theory of diffusion and reaction in porous media. AIChE Symp. Ser. 1981, 71, 46-62. (33) Whitaker, S. Transport processes with heterogeneous reaction. In Concepts and Design of Chemical Reactors; Whitaker, S., Cassano, A. E., Eds.; Gordon and Breach Publishers: New York, 1986; pp 1-94. (34) Quintard, M.; Whitaker, S. Transport in ordered and disordered porous media: Volume averaged equations, closure problems, and comparison with experiment. Chem. Eng. Sci. 1993, 48, 2537-2564. (35) Ochoa-Tapia, J. A.; Stroeve, P.; Whitaker, S. Diffusive transport in two-phase media: Spatially periodic models and Maxwell’s theory for isotropic and anisotropic systems. Chem. Eng. Sci. 1994, 49, 709-726. (36) Quintard, M.; Whitaker, S. One- and two-equation models for transient diffusion processes in two-phase systems. Advances in Heat Transfer; Academic Press: New York, 1993; Vol. 23, pp 369-465.
Received for review April 2, 1998 Revised manuscript received August 28, 1998 Accepted September 15, 1998 IE980212T
