Directional Dispersion Coefficients in Anisotropic Porous Media Michael B. Moranville, David P. Kessler,' and Robert A. Greenkorn School of Chemical Engineering, Purdue University, West Lafayette, lndiana 47907
The effect of anisotropy of porous media on the dispersion of miscible fluids is investigated and the form of the dispersion coefficient obtained is shown to be consistent with results obtained using a single-continuum model. For transversely isotropic porous media the dispersion tensor can be characterized by six scalar invariants. The form for the directional dispersion is developed in terms of the six invariants for several sets of boundary conditions. Constraints on these invariants are developed. Although the dispersion coefficient in the direction of flow is a nonlinear function of cos2 8, the deviation from linearity is not large.
Introduction Porous media occur in great abundance and variety. All living organisms rely on their porous nature for their life functions. Soil is porous and so are most natural rocks. Ground water, oil, and natural gas all are found in porous formations, and their efficient recovery requires a knowledge of fluid dynamics in porous media. Porous media find use in virtually every aspect of the chemical processing industry from separation processes to reactor systems. Almost all porous media possess, to some degree, the nonidealities of nonuniformity, heterogeneity, and anisotropy. It is our objective here to characterize the effect of one particular nonideality, anisotropy (dependence of properties on direction), on dispersion. The effects of dispersion are important in the study of the miscible recovery of oil and in the movement in ground water systems of contaminants such as radioactive wastes, heavy metals, and salt water. Dispersion, as considered here, is the macroscopic mixing of fluids during laminar flow through a porous medium. Consider a porous bed in which a fluid containing a tracer of concentration C flows a t a steady rate as shown in Figure la. At time zero change the inlet concentration to COand monitor the fluid concentration downstream a distance L . If there were no mixing or wall effects in the bed, the outlet would show a step change in concentration; if the only mode of mixing were molecular diffusion, a narrow mixed zone would result as indicated by the dotted line in Figure Ib; however, actual breakthrough curves indicate that much more mixing occurs than can be explained simply by diffusion. This additional mixing, which is referred to as dispersion, is the result of spatial fluctuations in the velocity vector as the fluid traverses the complicated network of pores within the solid, and will occur even in a porous medium whose local average properties do not depend on location, direction, or time (an ideal medium). The usual treatment is to model porous media not as a set of two continua, that is, fluid and solid, but instead as a single continuum (fluid) in which the concentration is that of the fluid phase, the velocity is the superficial velocity, and the dispersion is characterized by a dispersion coefficient defined by adopting the mathematical form of the diffusion equation (continuity equation for a species). For a general, homogeneous, porous medium the continuity equation is (De Josselin de Jong and Bossen, 1961; Nikolaevskii, 1959; Scheidegger, 1961)
where Di, is the dispersion tensor. Several investigators (Bear, 1961; De Josselin de Jong, 1961;Guin et al., 1972; Nikolaevskii,
1959; Scheidegger, 1961) have suggested that Dij is related to the components of the mean velocity by
The fourth-order tensor a y k l in general has 81 components, which vary with medium, fluid, and velocity. Determination of such a large number of components obviously represents an intractable problem experimentally. Therefore, effort is made to draw conclusions from simple systems which nay be used to approximate real porous media. The number of components required to describe the dispersivity for a particular anisotropic medium will, of course, depend on the type of anisotropy the medium exhibits. Scheidegger (1961) has suggested two symmetry relations which should apply to all forms of anisotropy. The first of these aijkl
= aijlk
is obvious since such a reversal of indices does not affect the dispersion equation. In obtaining the second arjkl
= ajikl
Scheidegger applied Onsager's principle of microscopic reversibility. Although the argument for this second symmetry relation is rather weak, we use the relation in what follows. These two symmetry relations reduce the 81 components of arjk[to 36. However, by making use of additional symmetry properties of a particular medium, a much smaller set of scalar invariants can be developed which will describe the dispersivity for the medium. For an isotropic medium, the dispersivity tensor a [ j k ( can be described in terms of two scalar invariants so that
(3) where a1 and a11 can be functions of the medium and fluids and may vary with the magnitude of the velocity. Equation 3 has been verified by numerous experiments. Transversely isotropic media are those whose properties remain invariant under the set of transformations IS)where
(4)
The bar denotes a transformed coordinate system and ap,,ol,,,, =
dpu; detb,,,) = 1
Ind. Eng. Chem., Fundam., Vol. 16, No. 3, 1977
327
Porous Bed a:O
a = L
(0)
I . At
L Figure 2. Layered medium as an equivalent anisotropic medium.
(b)
Figure 1. Definition of dispersion. Here p , p , and u take on the values 2 and 3. It has been shown that for a transversely isotropic medium the dispersion tensor is given by [Moranville (1974) and Poreh (196511
Di; = ao6;;u
+ a1XiXju U(iM + a4 U U’ U + as[uix; + U j X i ] U ’
(5)
face is also given in the Appendix. Using this solution to compute breakthrough curves, we then can determine an effective dispersion coefficient for the layered system of Figure 2 as follows. It is also possible to consider the system of Figure 2 as an equivalent, transversely isotropic medium. Shamir and Harleman (1967) showed that for flow perpendicular to the layers, the breakthrough from a layered medium can be modeled by an equivalent dispersion coefficient. (Marle et al., 1967) did the same for flow parallel to the layers). For a step change in concentration along x 1 = 0, the dispersion equation in terms of equivalent properties becomes
where U(X) = UiXl
Here again the coefficients are scalar invariants which may depend on u and uix). Poreh (1965) has shown that for a medium having an axis of symmetry, i\, the dispersion tensor, is given by
D;j = B16ij
+ R2uiuj + BsXiXj + B i ( ~ i X+j ujXi)
The solution to eq 7 with boundary conditions \k = 0: t = O , X l > 0
9 = l:t
9 = 0: t
(6)
The Bj’s are scalar invariants which may be functions of the geometry of the medium and the fluids flowing. They also are arbitrary functions of the magnitude of the velocity ( u ) and the component of the velocity along the axis of symmetry (ux = UfiXk).
Since transversely isotropic media are a subset of all media having an axis of symmetry, it must be possible to put eq 5 in the forni of eq 6. This can be done by letting
Two-Dimensional Flow in a Semiinfinite Layered Medium Equation 5 has six scalar invariants, eq 6 has four. We propose that the ai’s in eq 5 have only slight dependence on u(x).To test this assumption we consider two-dimensional flow in the semiinfinite, layered medium shown in Figure 2. This is a system for which we can write a mathematical model as indicated in the Appendix. The solution to the dispersion equation for a step change in concentration along the upper 328 Ind. Eng. Chem., Fundam., Vol. 16, No. 3,1977
0, x1-
m
is
where
Equation 8 can be fitted to the breakthrough curves computed using the method in Appendix A to determine a value for D I ~ ( However, ~). from eq 5, using = u ( @cos ) B ( e ) , we know that we must have
D1l(e)= [(ao + a l )
It should be stressed that the equations above apply only to the particular type of anisotropy for which they were derived.
>O,r1=0
+ cos2 ~ ( e(a2 ) + a3 + + 2a5)]u(c2)(9) a4
Thus, if the ai’s were completely independent of u(x),the equivalent dispersion coefficients calculated from the breakthrough curves computed for the layered system of Figure 2 would be linear functions of cos2 Oc.1 with intercept (a0 a l ) and slope (a2 a3 a4 t 2a5). Figure 3 shows the values of D11Ie)resulting from using the data given in Table I to calculate breakthrough curves using the method in the Appendix and then fitting eq 8 to obtain Dll(e).The solid lines are the linear least-squares fits to the calculated points. From this figure we see that while strictly speaking is a nonlinear function of cos2 the deviation from linearity is not great. Thus, at least for the system of Figure 2, assuming the a,’s independent of u(x)does not lead to great errors. The value and mathematical form of directional permeability and directional dispersion coefficients in anisotropic media vary with geometry, so we shall now look at the form of the directional dispersion coefficients for various geome-
+
+ +
Sei
0
I
Sei 2 Set 3
0
0 A
/
Table I. Data Used in Determining Curves of Figure 3 ~
~~
h(" x 106, cm2
k'2) x 106,
cm*
cm/s
2
3.0 3.0
30.0 30.0
3
0.3
3.0
0.010 0.005 0.010
Set
/
no. 1
.(e),
Here
Dli = (ao + a l ) u (15) 0 2 2 = (ao + a4111 For flow in the X Z - X ~ plane with C perpendicular to the axis of symmetry, the solution is \k=
50
04
02
06
cos*
0.8
I O
0
Figure 3. Effective dispersion coefficients for layered media.
with 022
tries, and, assuming the a's to have negligible dependence on velocity as above, we can develop a set of constraints on the a's. Because of the tensorial nature of D,, it is independent of geometry. The theory leading to eq 5 does not give any information concerning the values or relative magnitudes of the a, Is. However, any diiectional dispersion coefficient must be positive. We therefore develop the dispersion coefficients in various geometries and apply the requirement of positive sign. Point Source in an Infinite Medium When the principal directions of D,, are taken as coordinate directions, the dispersion equation is
033
a*
axl
I
%
r
I
1
.I
I
I < , vi :+
n
* = 8(a"t
a*
uq--
arc,
a*
u:, - (10)
an,
to eq 19 in an infinite medium with a point source
Q located at the origin is Q
a.
+ a l + a 2 + + a4 + 2a5 > 0 + >0 + >0 + a4 > 0 a3
+ For continuous point source in a transversely isotropic medium with ri parallel to the axis of symmetry
* = 4Tt-
Q
For flow in the x 1 - x % plane with ri perpendicular to the axis of symmetry, the form of the solution is
4at-
a2
(19)
a0
a1
(20)
a0
(21) (22)
aa>O
Layered Systems. From the analysis of two-dimensional flow in a layered system above (eq 9) if a0 a1 > 0, then
+
+ a3 + + 2a5 > 0 a4
(23)
Goad's (1970) measurements in layered systems indicate that the dispersion coefficient is greater for flow perpFndicular to the axis of symmetry than for flow parallel to A. Therefore a1
+ a2 + a3 + 2a5 < 0
(24)
These last two restrictions imply that
> a1
(25)
For completeness, we now consider two other geometries which give some of the same but no additional information, but are the geometries most used to obtain experimental data. Dispersion in Transversely Isotropic Slabs Consider a semiinfinite slab taken from a transversely isotropic medium so that the axis of symmetry X makes an angle p with the normal to the slab's face, it as shown in Figure 4. If the upper face of the slab is maintained a t a pressure POand the lower face at a pressure PL, the pressure gradient will be parallel to it. The component of the Darcy velocity parallel to it will be
From eq 5
3=
(18)
a0
a4
'011022033)
(17)
Therefore, from eq 13, 15, 17
a2
- u1--
= (ao + a 4 ) ~ = uou
Sf 0
q(ni = QIA
where QIA is the flux per unit area. First we consider the permeability. The directional permeability associated with it is defined as
Ind. Eng. Chem., Fundam., Vol. 16, No. 3, 1977
329
-D ( n ) - a0 + a1 cos2p + a2 U
+ a3 2):(
P = P,
(3 cos2 p
+ a4
Figure 4. The velocity vector in a transversely isotropic slab.
where
(y2
Rearranging eq 30 and making use of eq 29 D(n) - A0 + A I cos2 P + AII cos4 p --
Scheidegger (1954) and Greenkorn et al. (1964) have shown that for an infinite slab k(n)
(26)
= nlkLJnJ
When the principal directions of k,, are taken as coordinate axes, eq 26 becomes
U
p + r2 sin2 p]
[COS~
where
+ a4) A I = (a0 + a2) + 2r(a4+ + r2(al - - 2a4) A11 = (a1 + a3 + a4 + 2a5)- %(a4 + ab) + r2(a4 - a l ) A0
= r2(ao
ab)
k ( n )= k l l cos2 p + k 2 2 sin2 P
These relations yield no new constraints.
where /3 is the angle between ri and the major axis of klJ. We now consider dispersion. If a step change in concentration occurs along the upper face of the slab, the concentration gradient will be parallel to ri and may be written
Linear Flow in a Long, Thin Cylinder
a* -- n, a* _ ax,
an
a a* --=
a2*
nlnJ 7 an
ax, ax,
(31)
a0
Now consider a long, thin cylinder of material taken from a transversely isotropic medium in the ri direction where ri makes an angle with the axis of symmetry, X. If the diameter to length ratio of the cylinder is sufficiently small, the mean velocity will be parallel to ri, there will be no net flux of material normal to ri and the dispersion equation may be written (for a coordinate system with the x 1 axis parallel to ri) as (32)
and eq 1 becomes
a*
a**
at
an
- = D(,,) -- u ( n )
a*
(27)
where
where u ( n ) = nlul = q ( n d 4 D(n) = niD,jnj a similar solution to that for permeability. The boundary conditions for this system are
The solution to this equation for a step change in concentration at x 1 = 0 is also of the form given by eq 8. Thus, from breakthrough experiments in linear models taken from a transversely isotropic medium in various directions, it is possible to determine D’(n) as a function of p. Goad (1970) performed these experiments and found that a plot of vs. p had a nearly elliptical shape. From eq 5, when ri lies in the x1-x~plane
\E=Ofort=O,O
