- Email: [email protected]
Numerical studies of transport in porous media
Ann. Nucl. Eneryy, Vol. 20, No. 3, pp. 203-224, 1993 Printed in Great Britain. All rights reserved
0306-4549/93 $6.00+0.00 Copyright © 1993 Pergamon Press Ltd
NUMERICAL STUDIES OF TRANSPORT IN POROUS MEDIA
R. L. BUCKLEY and S. K. LOYALKA Mechanical and Aerospace Engineering and Particulate Systems Research Center University of Missouri-Columbia, Columbia, MO 65211 (USA) and M. M. R. WILLIAMS Electrowatt Engineering Services (UK) Ltd. Horsham, Sussex, United Kingdom
(Received September 8, 1992) ABSTRACT
Williams has recently developed a new formulation for flow through porous media by using a modification of the kinetic equations of particle transport. We explore in this paper several numerical techniques for solution of this equation. The methods explored include the upwind, the Smolarkiewicz method [19831, and the NND schemes of Zhang and Zhuang [1992]. The numerical methods were benchmarked against an exact twodirectional solution (the forward-backward model) and it was found that of these finite difference methods, the NND schemes, developed to accommodate shock waves, reduce the inherent diffusion error the most while avoiding spurious oscillations. The numerical techniques were further explored for a multi-directional case for which exact solutions cannot be obtained, and it was again found that the NND schemes outperformed the other methods. INTRODUCTION A porous medium consists of a solid phase and interconnected pores (or voids). These voids, which may be either liquid or gaseous in nature, generally have no distinct pattern (i.e.the geometric and hydraulic properties are spatially variable, or heterogeneous), which makes the problem of describing the transport in a porous medium a tedious one. Porous media transport depends on many factors: the physico-chemical interactions between the solid material and solutes (void-fluid), the often unsteady and non-uniform flow conditions, the complicated geohydrological structures of the aquifer, as well as the mechanisms by which the solute spreads [Dagan, 1987]. The phenomenon of transport of flow through a porous medium is relevant in numerous areas. "Groundwater," for example, a term used to describe the flow of water through the Earth's crust, has long been an area of intense study. The movement of pollutants through ground water is of increasing importance in considering alternatives for disposing of nuclear waste. In a more general sense, those working in such fields as hydrology, chemical engineering, environmental engineering, and soil physics [Dagan, 1987] are concerned with tracking the motion of the fluid contained within the voids in an accurate fashion. This knowledge is important in many applications, such as the groundwater application above, where 2O3
204
R . L . BUCKLEYet al.
pollutants are created as a byproduct of the fission reaction, as well as the accurate knowledge of deposition of pollutants into the surrounding environment (whether it be through complicated urban building structures, or the intricate branching of vegetation, for instance). Due to the natural intricacies found on the microscopic scale, engineers and scientists are often concerned with the average values of the flow or transport variables. In other words, the material is regarded as a spatially random medium, which requires stochastic analyses. Such statistical analyses require a simplification of the equations describing the flow phenomenon. The starting point for analytically describing the complicated porous media flow of a species is the advective-dispersion equation [Williams, 1992a; Bear, 1979]:
0C & =
-V. J - tC
(1)
J = - D V C + UC
(2)
where
with C {= C(r, t)} denoting concentration, D {= D(r)} being the dispersivity tensor, U {= U(r) } representing flow velocity, and ;Ldescribing the decay term. This equation gives the transport due to overall advective motion as well as that resulting from the variable paths taken by the fluid flowing through the porous medium. It is assumed that the motion due to diffusion obeys Hck's Law. The dispersion tensor has been determined experimentally, and from mathematical models, to take the form [Bear, 1979]:
Dij =aTSijU +(aL--aT)
U
(3)
where aT and aL are the transverse and longitudinal dispersion lengths, respectively (which are functions of the porous medium geometry), and Ui are the orthogonal components of the flow velocity. Now, if the dispersion lengths are to be considered fundamental parameters of an experimental system, their values should be independent of measurement size or time. However, in the experimental determination of D, it has been found that at, is dependent on the distance for which measurements are taken (i.e. the value of aL depends on the "size" of the experiment) [Williams, 1992b]. This indicates that dispersivity increases in proportion to the distance of particle flow travel and the prescription (3) is not adequate for description of flow in both small and large scale porous media. It is therefore of considerable interest to investigate eqn (1) without use of the prescriptions (2) and (3) Often, a stochastic (random) averaging technique based on the continuum concept is used to solve eqn (I). One obtains a mean value equation which describes the scale dependence of aL. While this approach provides a solution that is non-Fickian, it is also expensive to implement in practical situations. Obtaining random properties for the medium is difficult, and the results acquired from these methods are commonly used only as indicators of general trends.
Numerical studies of lranspor! in porous media
205
Williams has thus proposed a new model based on the discrete kinetic theory of particle transport and not on continuum concepts used in dispersion theory. The stochasticity of the process is incorporated by assurmng the particles are subject to random collisions. While the model equation can be solved analytically for simple cases, in general it would be necessary to employ numerical techniques to solve the equation. Thus, our goal in this research is to explore various numerical techniques to solve the new equation and to examine the strengths and weaknesses of these methods. NEW FORMULATION AND ANALYTICAL SOLUTIONS The new formulation [Williams, 1992a, 1992b] is based on an equation for the species concentration dependent on space, direction, velocity, and time (r, [1, v, and t, respectively). After various simplifications to the transport equation, [Williams, 1992a], the following integro-differential equation is obtained for the 'angular concentration' C(r, ~ , O:
[R~+vfl.V+vZ(r,~)+)lJ~]C(r, fl, t) (4)
+1
= f(r, fJ)fdff
~(r, f f )C(r, f f ,t)+ S(r,O,t)
-1
where S represents a source term and R is the retardation factor. The new concepts involve f and ~. The scattering function, ~, represents a reaction rate. In essence, it is the chance that a particle can be taken out of the flow media (such as by deposition onto the solid surface, or via a chemical reaction). The anisotropy function,f, on the other hand, represents a redistribution of the flow field. It answers the question of where the flow travels when a branch or fracture intersection is encountered. Thus, in onedimensional transport, we have: [~+
d vlt--~ + v ~,(X,la ) + Z ]C(x,lz, t) = vf (x,lt )!÷tdlt" ~.(x,.')C(x,lt',t) + S(x,lt,t) (5)
where/~ = cos 0 and 0 is the angle that the particle makes with respect to the x-axis as it v flows through the structure. The velocity v is now ~. Williams has made substantial progress in obtaining analytical solutions to the various approximations of eqn (5), but for realistic geological situations, numerical methods are required. We begin by assuming ~ a n d f a r e spatially constant. A complete description of a porous medium would involve an anisotropy function of the form [Williams, 1992a] K
f(/z)= Y~akt0z-/~D k=l
with
(6)
206
R . L . BUCKLEYet al. K
~a t = 1
(7)
k=l
where k refers to the direction taken by the particle and 'a' with this probability, 0 < a < 1. We also assume the source term is:
(8)
S(x,//,t) = So(X,t)f(gt) The resulting coupled first order partial-differential equations are:
( t9
+ v//k
/
'
°~ +vz~(//t)+/q, Ct(x,t)=vat~-r,(//j)Cj(x,t)+So(x,t)at. Tx :=1
where k = 1,2 .... ,J. By assigning / 3 t = v / / t , g / = v , ~ ( / / / ) + ~ , , and at) = vat ~(//j), we obtain a more compact form:
(9)
St(x,t)=atSo(x,t)
'
~+flt
+~'t Ct(x,t)= Y.°tkjCj(x,t)+St(x,t)
(10)
j=l subject to the initial and boundary conditions
Ck(X,O) = gk(X) such thatx ;~ 0, and Ck(0,/) = hk(t), where the functions g and h are known. An analytical solution can be obtained for eqn (10) if all cz are set to 0, which removes any coupling of the equations describing flow in various directions. We choose to examine a special case of eqn (6) which retains the directional coupling. Referred to as the forward-backward model, this scheme is commonly used in neutron transport and involves only two directions,//x and//2. As the name would imply, this simple scheme assumes the species flow is restrained along a crack inclined at an angle 01 (between 0 and 90 °) relative to the fluid flow direction, with the possibility of back-flow in the opposite direction at an angle of 02 (90 ° < 02 < 180°). This leads to the anisotropy function being written as [Williams, 1992b]: f ( / / ) = a S ( / / - / / t ) + (1 -
a)6(//+//2).
(11)
We seek a solution of the form:
Cr(x,t) = C~(x,t) + C2(x,t).
(12)
Assuming a source term as in eqn (8), Williams [1992b] has shown that one can obtain the coupled set of equations:
o3 + (I - a)~q+~]C1(x,t)_aZ2C2(x,t )= aSo (x, t) [~o3+//~-~
(13a)
Numericalstudies of transport in porous media d B [~-laz-~+a~2+~.]Ca(x,t)-(1-a)~,~Cl(x,t)=(1-a)S,(x,t)
207 (13b)
where it should be noted that ~1 = ,,~(/.tl), .,~z = ~(-/t2). In addition, velocity has been absorbed such that vt--~t, 2/v--*Jt, and vC-~C. These equations can be combined to yield a generalized Telegrapher's equation:
~'--~+t--~-~+ sc-~t]C.(x,t)=[D--~2-U ~x - ~ ]Co(X,t)+Qo(x,t)
(14)
where .St = a,S2+(1- a)2~t +Z, ,~,'~
1
v~, r " t = I~! -1~2
'
Z,r+Z Zr #1#2 v,
and U = ~ { v~--"(., -/22) + a ~2//, - (1 - a).~t/.t2 } • Now, if Co(x, t)= Cl(x, t), then the source term is:
whereas, if
Co(x, t)= C2(x, t), then
Qo(x,t)= r ( a - a ) [ ° + ..~+(/~+~l)}]So(X,t). ,
(16)
Summing eqns (15) and (16), or solving eqn (14) with Co(x, t)-- CT(X, t), yields
Qo(X,t)=[--~r{~-V*ff-~}+l]So(X,t) where
(17)
208
R.L.
BUCKLEY et al.
(1 - a)/~l
V" = a/z 2 -
for the total concentration. Williams has obtained an exact solution to eqn (14) by assuming a spatially infinite medium, x e (-00,-0), a source of the form
(18)
So(X,t)= go<~(x)S(t),
with S(t) = 1, t > 0; 0, otherwise; and using Green's function and Laplace-Fourier transforms. Assuming ~ = ~1 = ~2, and utilizing the following:
I:~1
( ~2
1
1
U.e'~2
r
7
"D
4D2.)J
! ~=
+
+ ~
2x and C=
2rd)+Ul £2 + 4 D ' r '
with
• =xZ
~=vZt
D= vD° Z
U=vUo
~= ~° v£
v~=vq
as well as exp
v
- --
~j£2 + 4DoXo 1
g(x,t)-tt, and
v+Jk v;JJ
l = £_.~0
Z #=VZ#o
V ° = yV;
C
=
V#~,C o
Numerical studiesof transport in porous media
209
1
tk
v; j
k
v2 2J
he obtained for CT(x,t):
{/
v:l
X ~o l + - ~ - + l sr_ l t - - ~•- lIe x p lr-c.~l .-TTl+
v; ) k
'
ls(7_Fo)exp(_Cffo)A(~,F°)d~ ° L v; j ._
v; j
}
if o < .~ < v,%
Cr(x,t )
=
{iv.
_
_,
Z "Co 1----e-°]s(,+x-~--]expl cox1+ ~s(t-t°)exp(-Cffo)fr(~,~°)dE ° if-vj
v;3
v;/k
LV;j
}19,
<:~_< O,
where
fr(X't)=El-'Co(Co+(U°-Co¢)V: 2D° )]Io[fl°g(x,t)] (20) +--'Co[3o 1 P(x,t)+ 1+ 2 v-~j
l~[fl, g(x,t)]
and lo[flog(X,t)] and ll[flog(X,t)] are Bessel functions of the first kind of order zero and s(z) -- 1 if z _>0, and equal to zero, otherwise.
one, respectively. Also,
Extending the analysis for Cl(X,0 and C2(x,t) we obtain:
fs0-io/eXp(-qo4(~.~o,aio}'
[ZJa*o(l+~Llt-=]expr-~-+ "" - • -Co~l+ I[~ ~j~ v~j LV~,~ ifO < ~
5 v~i,
q(x,t) =
Z a'¢o 1-
t+
i f - v 2 i _ < ~_< O, with
exp
+ s(~-~o)eXp(-Co~o)ft(~,~, )
R.L. BUCKLEYet al.
210
(22)
+~a'Coflo((ll~2ip(x.t)+{lk + MVo )P(x,t)J _.'-!J,,,ao~,x.,,, z Lk v'~) and, similarly: tI .~ [(1- a)'ro(1./"Vo) 27.+k- s(t-=]od \ v~,)
]
..-.=1 -Co~
L v~, j
+ (?- t-o)exp(-Co~o)B(~,~o)di-° if0 < ~ <
v~?,
C2(x,t) =
v~ )ik"
v'~)
LV~j (23)
Z
+ ~_~(?-L)exp(-cJ°)f2(X,to)dto
i f - v ~ ? < ~ < 0, with
f2(x,t)=[(1-a)~o(~,+Z,1)-(1-a)~o(Co-(U°-c°£°)~lll 2D° Io[flo(g(x,t)] (24)
+1(1-a)'¢°~°[(a+lal)P(x't)+lx-~-~'l+l-'l-'~]ll[fl°g(x't)]'v-~) ~, v; )P(x,t)J NUMERICAL TECHNIQUES The necessity for developing efficient numerical schemes for the transport equation is evident in th:' analytical solutions are feasible only for special cases. In this section, we examine several numerical techniques based on finite differencing. The numerical results are compared with the exact solutions for the forward-backward model as described by eqns (19), (21), and (23). A. UPWIND METHOD We first describe the simple upwind scheme. For clarity, we will assume v to be constant. When applied to the wave equation,
N u m e r i c a l studies o f t r a n s p o r t in p o r o u s m e d i a
211
(25)
OC(x,t) ~ v °~C(x't) = 0
Ot
dx
one utilizes backward (or forward) spatial differencing and forward temporal differencing to obtain
c~i "÷' -C7 I- v C , -C?_, At
= 0
Ax
C.'_+' - C 7 ~-v C7+' - C7 = 0 At Ax
i f v > 0,
(26a)
if v < 0
(26b)
where i refers to the spatial step and n to the temporal step. Often, the terms sought in these numerical computations are those containing the (n+l) time step, since they represent the unknowns in time, and the time step is generally the outermost loop used in finite difference codes. It is assumed that from initial conditions, all values at the (n) time level are known. Eqn (26) has a truncation error of O[Ax, At] and is thus first order accurate. We call this form explicit since only one term, t,-,,t+l ~ , is unknown at each time step. This equation is stable provided the Courant number, r/=[~-~
(27)
The method suffers from numerical dissipation (amplitude) error which is typical of methods with first order accuracy. However, its extreme simplicity lends itself to easy programming. A variation of the upwind method described has been utilized to model eqn (10) as given below: (t~n+l
r,n
t'i,k.~i,k
n
n
~ , on (C;,k--C;-1,k ~a_an t',n+'
J~ ,,,,,n ¢.-,n+l 'K = Z.~ mi,kj~'i,j
(28a)
if ]],'.k > O, and qk
-qk
.
J •
~,
-.n
t~n+l
+
Sink
(28b)
j=l
if .8;:.
o.
The added subscript k denotes the direction taken by the particle. Note the use of concentration associated with a "new" time step for the decay term and the coupling term. This is a modification to the usual procedure which has been devised to avoid instability in the solution.
212
R.L.
BUCKLEY
et
al.
Due to the presence of the coupling term, we have an implicit scheme based on direction with a system of J equations and J unknowns. For instance, in the forward-backward model, we have two directions. Thus, in eqn (28), J = 2 and we have the coupled set of equations: ("~n+l ~i.l
--
C/~.I'
"4" ~i,1
Ca+l
CT, l _ C ~ n _ i , l
n
At ~1
i.2
~i,2
At
. o.
+ t..'~.2
n n+l "I- Ai,ICI,I --
Ax Cin2
t.-,,~t
~-1,2
'
~ ..,.
n ,.i. Si, 1 _
• ["~n'l'll_#,.iCn (",.+I ai,ll,.,i, 1 _ _i,12~1,2
~_,.,. r.+;
.+i
-
(29a)
c..+l
(29b)
+ ,.,,:,,...:,+_+ = --Sl.~- -+.2~-+.,. + ~i,221-'i,2
Ax
if flu, fl,'a > 0, and similarly for 3~j, 3,'.2 < 0. We thus have two unknowns denoted by the (n + 1) superscripts. Solution to a general set of equations involves using matrix operations. Thus, eqn (28) can be solved numerically using the following matrix representation: [A] (y) = (c)
(30)
where [A] is the matrix of known coefficients, (e) is a vector containing known parameters, and (y) is the vector of unknowns (i.e., the "new" concentrations). In general, these components can be written as: tt
tl
II
1 + Ai: - ai,u At n
--O[i•21
-oci,;2
-ot~,t~
1 n n At + Ai '2 - Cti.22
-Ocl,z~
[A] =
-a~m"
" :
: a
n
--~i,J1
--O~i,J2
c?., + s:. ' At
c?,2 (c) =
+ s .2 - / L + s?,, -
II
-o~i•ls
tl
-O~i'as (31)
1 + A,'.3
At
:
.
.
.
.
~i,3J
.. n
--O[i,j3
..
1 At
n
4"
n
t~i, J -- OCi,jj
c?.,-c_?_,,, ,'i.~
n
...
Ax
n n c?,2 - c i-1,2
(32) C,'.~ - C?_ u
--
if fl > 0, and similarly for non-positive values• Finally, (y) is
Numerical studies of transport in porous media
213
(33) i
"
Lc,-.;'j Gaussian elimination and back substitution are then used to obtain (y). This procedure is carried out first through all the spatial sections, x, and directions,/~, for a given time step. Once this is complete, the concentrations are updated and we proceed to a new time step. Though this modified implicit upwind scheme is more expensive computationally than an explicit method, it is necessary to avoid spurious results (i.e. negative concentrations). For the forward-backward model, we have to make use of the anisotropy functional form given by eqn (11). This is expressed in eqn (28) via the source term, which is assumed to be as given in eqn (8) where
So(x,t ) = -So~(x)S(t). For the numerical solution, we can thus take:
Sf,k =Sola~k,l~,i+(1-a)Sk.2S~_l,i]S(t)
(34)
where d/above is the Kronecker delta. B. SMOLARKIEWICZ
Smolarkiewicz [1983] noted the dissipative problems associated with the upwind model and utilized Taylor series expansions about the point (xi, tn) to closely approximate the hyperbolic wave equation, eqn (25), to obtain: OC(x,t) ~
o3t
o3[v(x,t)C(x,t)] o3x
- o3[ K ~ OC(x,t)] -~
Ox
(35)
where
(IvlAx-
K/,,~ =
vZAt)
(36)
is a diffusivity term which helps to eliminate the errors of the upwind scheme. This is essentially a predictive step. To correct this, he reversed the solution of the differential equation as
3C(x,t) ~- O[v~C(x,t)] = 0 Ot ~x where
(37)
R.L. BUCKLEYet al.
214
{
- K ~ . . BC(x,t) v~= C(x,t) Bx 0
if
C(x,t) > 0
if
C(x,t) = 0
(38)
is the diffusion velocity. In addition, it is convenient to introduce a function: 1 At
+
F(n~,n,÷~,~):~[(~
I~)ni+(~-l~)ni+t].
(39)
In order to use this scheme in our problem, we def'me:
(
'
= )gi,kCi,k -- E ai.kjCi,j + S j=l
At
(40)
where )~, a, and S represent decay, coupling and source terms, respectively. Note the use of "old" concentrations with the various terms. Instability does not arise, thus leading to an explicit scheme which is computationally cheaper. Utilizing the notation describing the upwind scheme, the predictive step can then be expressed as: n
F
n
n
n
n
n
~
n
where fin _1 n + ,I i+~/2,k -- "2 (fli,k fli+l,k) /~,,
1
.
(42a)
n
i-~.k : "~(fli-l.k + fli,k ).
(42b)
Essentially, fli+¼.k represents the velocity at a distance midway between mesh point (i) and (i+1) for a given direction k. After the predictive step, we apply
C i,k n+l = C i *. k -
F
*,k,Ci+l,k, *
,i+y2, k - F
-1.k,Ci,k,
,i_y2, k
+ Zi.k
(43)
to obtain the concentration at a new time step. The forward-difference "anti-diffusion" (or anti-dissipation) velocity in the above is expressed as:
fl~.i+~.k :
(Ci:, *C;+,.k + e)Ax
(44a)
Numerical studies of transport in porous media
215
where • is - 1 0 "15 to ensure one does not encounter a computational division-by-zero error. A similar expression is written for the backward-difference anti-diffusion velocity as:
+
+
(44b)
The maximum Courant number, I/, is given by the usual ratio where the velocity in question is represented by the absolute value of the maximum anti-diffusion velocity and the truncation error is of O(Ax2, At2). C. NND SCHEMES
Zhang and Zhuang [1992] have developed a series of numerical schemes which were originally designed to accommodate the complex flow fields associated with shock waves. Numerically modelling such shocks leads to spurious oscillations near or in the shock region. In the past, modeling attempts have been concerned with utilizing mixed-dissipative f'wstorder schemes containing free parameters near the shock wave, and second-order models elsewhere, with resulting unsatisfactory oscillations. Thus, Zhang and Zhuang have suggested a scheme which is nonoscillatory, contains no free parameters, and is dissipative (as opposed to dispersive, which implies oscillations), which they logically refer to as NND. To illustrate the scheme, consider the exact Navier-Stokes equation governing a steady flow with constant viscosity and containing a normal shock wave: ~'
c3C
. a2C
(45)
= o
subject to:
x ~ -**,C --.>C,,.,-~-- ~ O, where ~ is the kinematic viscosity and "o." refers to freestream conditions. The finitedifference form of this equation can be expressed as: 0C
dt
OC
.32C
~ B~C
. 03C
,. B4C
(46)
where the terms on the right-hand side of the expression are the truncation error associated with the discrete representation of the continuous function. To totally suppress the oscillations associated with numerically modelling eqn (45), Zhang found that the coefficient, ~'3, must be selected such that the Second Law of Thermodynamics is satisfied. A detailed analysis is made with the conclusion that if ~3 > 0 upstream of the shock and ~'3 < 0 downstream of the shock, the oscillations can be suppressed. In physical terms, if ~3 < 0 everywhere, then the increasing entropy condition
216
R.L. BUCKLEYet
al.
is not met in the upstream portion of the shock, and oscillations result in this region. The converse is true for ~3 > 0. Zhang and Zhuang further consider the equation: OC + alp(C) = 0.
at
(47)
ax
If we let p =/3C, then we have:
Oc
+
o~
= 0.
(48)
We define the function p as:
p = p+ + p- =/3+C +/3-C where the (+) quantities denote upstream locations and (-) quantities signify downstream locations with respect to the shock wave. Now,
2
2
which yields:
W+
(49)
+-L-x=0.
The following are then applied in all NND-schemes: (1) Upstream (a) 2nd order upwind difference in place of op+
Jx
(2)
(b) 2rid order central difference in place of 0pOx Downstream (a) 2nd order central difference in place of op+
ax
(b) 2nd order upwind difference in place of Op-
ax"
Placing this in the finite difference form of the equation, we have:
(50) where
Numerical studies of transport in porous media 1
217
P~+l- 1 rain m°d(APi-+~,APi-+~)
(51a)
h/_~ = P+-I + 1 rain mod(Ap+_3~,Ap+_/~)+ p~-- 1 min mod(ApT_~,Ap~ ~)
(51b)
h/+~ = p+ + -) min mod(APL~'AP++~ + 2
with ± - P~:, APi-~ ± = P~ - P~-l, etc... @ i±+ ~ ----Pi+l and
[Ixl if Ixl <
•
mlnmod(x,y)
:ly,I
lY[ x and y of like sign,
iflxl
(52)
> Lvl,x and y of like sign, ifx and y differ in sign.
These choices then ensure compliance with the Second Law of Thermodynamics. The schemes which result from the basic formulation given by eqns (47) to (52) are second order accurate in space (except at extremum points), and are dissipative since ~ < 0. In addition, the authors prove the semi-discrete NND schemes possess the TVD, or total variation diminishing property. Civen rc":Xlq
,,-crl
,then
,:0.
We
l
examine two of these schemes. L
NND-1
We use Euler first order time differencing such that
af. I n = C i,k n+l - C i,k n
(53)
where i denotes the spatial node, k the considered direction, and n the temporal step, as before. Modifications similar to that performed in Smolarkiewicz's method with the use of eqn (40) lead to the NND- 1 scheme being written in finite-difference form as:
Ci,k "+l=C"i , k - - ~At(h,, h" k i+y2'k-- i-Y2'k
)
"
+ ~,i,k
(54)
where h is as given by eqn (51) for a given direction, k. The maximum Courant number, r/, is 2/~ while the truncation error associated with this method is of O(Ax 2, At). Note the use of an explicit form of solution with regard to direction.
218
R.L.
2.
BUCKLEY
et al.
bIND-2:
This is a two-step predictor-corrector scheme. We have, C~,~-C n -
At (h"
-n
~+ n
(55a)
1[ n +C~, t A t ( h . * )] 1 * cn+li,k = "~ Ci, k -- " ~ , i+~/2,k -- hi_~2,k + ~Zi.k
(55b)
where the * quantities are temporary predictor concentration values used in the correetor step. In this instance, 77 = 1 and T.E. = O(Ax2,At2). The second order accuracy in time makes this scheme an improvement over NND-1. R E S U L T S AND D I S C U S S I O N S A. FORWARD-BACKWARD MODEL Simulations for the four numerical schemes presented have been performed and compared to the analytical solution. The parameters selected are as follows: a = 0.9, ~1 --- ~r2 = 1,
[09] [o,1 [091
/al = 1£2 = 0.9, ~, -- 0, and v = 1 leading to fl = -0.9 •
= 0.9 ' S = 0.1 6(x), and
a = [001 009]. That is, eqn (28) is explicitly:
oaC~(x,t) + 0.9 BC~(x,t) + 0.1C~(x, t) = 0.9C2(x,t) + 0.96(x)S(t) Ot Ox
(56a)
o3C2(x, t) _ 0.9 0C2(x't) + 0.9C2(x,t) = 0.1C~ (x,t) + O.18(x)S(t) Ot dx
(56b)
where the source term with regard to//(x) is as described by eqn (34). 1.2
Exact
-
1.0 =
. 0.8-
- .........
Upwind
----
NND-I
- .......
Smolar.
1.5 Exact
~
0.62 o
Upwind
o
0.4-
+
1.0
. . . .
NND-I
........
$molar•
~ 0.5
0.2-
_J
0.0
! -1.0
-0.5
0.0
!
!
0.5
1.0
0.0
' '
-10
I
I
I
I
-5
0
5
l0
X
X
Figure l a : F o r w a r d Direction C o m p a r i s o n at t = 1.0 sec
Figure l b : F o r w a r d Direction C o m p a r i s o n at t = 10.0 sec
s t u d i e s o f t r a n s p o r t in p o r o u s m e d i a
Numerical
219
Results were obtained for the four schemes (upwind, Smolarkiewicz, NND-1, and NND2) at elapsed times of 1.0, 10.0 and 100.0 seconds and compared with the exact results to illustrate the dissipative nature of the numerical schemes with time. Figure 1 shows the nature of the forward solution (C1), while Figure 2 describes the nature of the backward solution (C2). We note that for this particular case both of the NND methods lead to similar results, so the accompanying figures contain only the results of NND- I. To enable a meaningful comparison, we used the same timestep, At, for all methods. (The timestep was varied only with respect to the total elapsed time). We discuss the nature of our results as obtained by each method below. 1. E X A C T
The exact solution has a jump in concentration (resulting from the clef'tuition of the source term) at x = 0 as indicated in Figs. 1 and 2. As expected, the exact solution exhibits sharp breaks at or around the jump for both the forward and backward directions. 1.5 -
F~tswt
° 11 o
~
Upwind
mmm~
NND-I
; 1.o
Smolar.
iI
~ 0.5
J
0.0 - 100
-50
!
|
50
100
Figure lc: Forward Direction Comparison at t = 100.0 sec As the total time increases, one can see that an asymptotic solution is approached. For the
/
0.15
I= o
=
0.10
4=#
= 0.05 o
0.20
Upwind
Upwind . . . .
NND-1
........
Smolar.
0.15 a tw
NND-1
O U 0.05
|
-0.5
.... -.---.-.
". O.lO
0.00
-1.0
Exact
0.0 x
0.5
1.0
Figure 2a: Backward Direction Comparison at t = 1.0 see
ii •
0.00 -10
i -5
i 0 x
i 5
|
10
Figure 2b: Backward Direction Comparison at t -- 1 0 . 0 s e e
forward direction, the maximum concemration occurs after the positive x-region is entered. It is interesting that while for small times the backward direction peaks just before x ffi 0
220
BUCKLEY et al.
R.L.
and then drops off, for large times the maximum once again occurs asymptotically at x just greater than 0 before decreasing at the larger distances. This phenomenon can be explained by the form of the exact solution as time approaches infinite. 2. UPWIND For this particular case, the implicit upwind scheme proved to be computationally the quickest. However, this is true only because our problem involves only two directions. For a more general problem containing 5 directions, for example, the use of Gaussian elimination and backsolving would become increasingly time consuming. Examining the figures, one will notice a 'rounding' of the edges near sharp comers as well as a drop in number concentration at the largest distances when compared to the exact solution. This indicates a dissipative (or diffusive) error resulting from the discretization of the finitedifferencing scheme. Again, such behavior is expected for the simple first order modification of the upwind scheme. 0.15 .........
Exact Upwind
....
NND-1
........
Smolar.
°"
e~
0.10 L., .,J gl
= 0.05 O
0.00
,
-100
-50
,
-J
.
0
50
100
X
Figure 2c: B a c k w a r d Direction C o m p a r i s o n at t -- 100.0 sec
3. SMOLARKIEWICZ This method incorporates an "anti-diffusive" velocity to alleviate dissipative problems. While the dissipation was reduced (in comparison to upwind), unsatisfactory dispersive (oscillatory) errors about the shock resulted. A particular advantage of the Smolarkiewicz scheme over the upwind method is that the former is less sensitive to varying timesteps. In other words, one can use a larger timestep (hence, less computational costs) and retain stability. 4. NND Figures 1 and 2 reveal the same dissipative errors occurring with the NND schemes as in the other numerical methods. However, the results are clearly superior. Not only is the dissipative error reduced, but the oscillations inherent with the Smolarkiewicz scheme are also avoided. Varying the value for the timestep, At, revealed that the NND-2 scheme was less sensitive to changes than NND-1. In fact, it was the least sensitive of all numerical techniques studied here. On the other hand, the NND-1 method was most effected by an alteration of the timestep, which is not surprising considering its need for a smaller Courant number than the other numerical schemes for stability. Although in this example the NND-1 model performed as well as NND-2, it should be noted that if the need for larger timesteps arose,
Numerical studies of transport in porous media
221
NND-1 would be the first to suffer instability in its solution. This is also not surprising since it is only In'st order accurate in time. B. MULTI-DIRECTIONAL TEST CASE A separate test case involving three directions (and no analytical solution) was developed to ascertain the utility of the new formulation and to examine the performance of the numerical methods. In this instance, the particle has the option of travelling at an angle of 10% 80°, or 135 ° with respect to the flow direction. The scattering cross-section, £ ~ ) , has been assumed to equal the absolute value of/~. The distribution function,f(/t) is given by eqn (6) where K = 3. The following values for probability of particle flow direction and scattering cross-section were assigned: Direction 1 2 3
St 0.9848 0.1736 -0.7071
2~(4u) 0.9848 0.1736 0.7071
a 0.4 0.5 0.1
F 0.9848 ] f l = / 0 . 1 7 3 6 1, / / L-o.7071J
With v = 1 and A = O, this corresponds to the numerical values:
7098481
l
70.3939 0.0694
z = / 0 . , 7 3 6 | , and a=10.4924 0.0868 0 ~ , ~ , / Assuming Ck(x,O)= 1 . 0 a n d n o [0.7071J L0.0985 0.0174 0.0707J source introduction, we obtain the results shown in Figs. 3 and 4. In Fig. 3, we show 0.50
1.50
0.40 @
..............
Upwind NND-1
/ A
....
NND-2
~
........ Smolar.
0.30
~!i
1.25 ~
/S
1.00
!~!iNND-I Upwmd NND-2 Smolar.
~0.75~
0.20
0.50"
¢J
0.10
. . . . .
0.00 0.0
0.2
0.4
0.6
;'~~~7
~0.25 ,
0.8
X
Figure 3a: Time = 1.00, Direction 1
].0
0.00
0.0
I
I
0.2
0.4
X
I
I
l
0.6
0.8
1.0
Figure 3b: Time = 1.00, Direction 2
results of the various numerical techniques at the final time of t = 1.00 sec, while in Fig. 4, a comparison of the NND-2 scheme and the simple upwind modification (for all directions) for final elapsed times of 0.50, 1.00, and 1.50 see is illustrated. Several inferences may be deduced from the figures. From Fig. 3, it is evident that the upwind scheme gives results that vary considerably from the other numerical techniques. The more diffusive nature of this scheme is depicted by the 'flatter' concentration curves.
222
R . L . BUCKLEY et al. 0.60 O m I.
~
030 -'
a
r
O.2O
lI O
0.10 -~ 0.00 0.0
,
,
0.2
0.4
X
s
,
,
0.6
0.8
1.0
Figure 3c: Time = 1.00, Direction 3 From Fig. 4, at the smallest elapsed time o f 0.50 sec, we see a peak value near x = 0 for the third direction. This makes sense when one considers the angle for this direction is
.
1~0]S
'-~
1.25
Dk. 1--Upw. --""o---
Dir. I--NND-2 Dir. 2--Upw.
0.75 o
0.50
1"1
it. t¢ ~- ~" I0,
---~'-"
Dir. 2--NND-2
....
Dir. 3--Upw.
- - ~t -
Dir. 3--NND-2
0.25 ~ 0.011
0.0
0.2
0.4
Figure 4a:
x
0.6
0.8
1.0
Three.Directional Variation at Time = 0.50 se¢
greater than 90 ° . With flow going in the opposing direction, particles will tend to get washed outmore quickly with distance. Similarly, for the 1 and 2 directions, the maximum values occur at downstream distances.
f
Dir. l--Upw. Dir. 1--NND-2 Dir. 2--Upw.
0.75
~ - ~ ......
o,o I "',j,-. 0.00
.
.
.
.
----
--
0.0
0.2
0.4
X
0.6
0.8
1.0
Figure 4b: Three-Directional Variation at Time = 1.1111sec
Dir. 2--NND-2 Dh'. 3--Upw. Dir. 3--NND-2
Numerical studies of transport in porous media
223
As the elapsed time increases, we see a gradual drop in concentration in all directions, with a more prominent decline in the In'st and third directions. This is expected since no source is introduced. With regard to directions 1 and 3, their loss is quicker than in 2 since their respective reaction cross-sections, Z, are larger and the probability of the particles taking these directions (via the parameter a) is smaller than for the second direction. As time increases further (not shown), we note a gradual shifting in the second direction's peak toward larger downstream distances. This indicates that the majority of the contaminant is being washed away in the direction of fluid flow. 1.50 oo o~
1.25 O al IN
Dir. l--Opw. Dir. 1--NND-2
1.00
Dir. 2--Upw.
0.75 /£I
o
........ • ......
0.50
Dir. 2--NND-2 Dir. 3--Upw. Dir. 3--NND-2
0.25 0.00 0.0
0.2
0.4
Figure 4e:
X
0.6
0.8
1.0
T h r e e - D i r e c t i o n a l Variation at T i m e = 1.50 sec
CONCLUSIONS We have discussed the importance of accurately modeling flow through porous media in the area of mass transport. The inadequacies of current simplifications to the dispersiveconvective equation have led Williams to develop an equation based not on continuum concepts but rather on the kinetic theory of particle transport. Due to the mathematical complexity of the resulting one-dimensional integro-differential equation, analytical solutions are difficult to produce. However, Williams has devised one for the simple case of flow in only two directions, the forward-backward model. We have explored in this report several numerical schemes which apply to more general physical applications. Comparisons with the analytical solution to the forward-backward case have been performed. Based on these comparisons, it is our assessment that the modified NND schemes are most promising due to their avoidance of oscillations and their characteristically smaller dissipative error. While the upwind scheme was computationally most efficient for the simple two-direction problem, it is apparent that the explicit forms of Smolarkiewicz and NND methods are warranted for more complicated systems involving several directions. A CKNO WLED GMENTS
Research performed by Robert Buckley was under appointment to the Environmental Restoration and Waste Management Fellowship program administered by Oak Ridge Associated Universities for the U. S. Department of Energy.
224
R.L. BUCKLEYet al. REFERENCES Anderson, D. E., Tannehil, J. C. and R. H. Pletcher (1984), Computational Fluid Mechanics and Heat Transfer, Hemisphere Publishing Corporation, 40-3. Bear, J. (1979), Hydraulics of Groundwater, McGraw Hill. Dagan, Gedeon (1987), "Theory of Solute Transport by Groundwater," Annual Review of Fluid Mechanics, 19, 183-215. Smolarkiewicz, P. K. (1983), "A Simple Positive Definite Advective Scheme with Small Implicit Diffusion," Monthly Weather Review, 111, 479-86. Williams, M. M. R. (1992a), "A new model for describing the lxansport of radionuclides through fractured rock, Part I: General Theory, "Report to CEC under contract F12W-CT91-0087. Also Annals of Nuclear Energy 19 (1992). Williams, M. M. R. (1992b), "A new model for describing the transport of radionuclides through fractured rock, Part II: Numerical results and extension to overlapping fracture sets," Report to CEC under contract F12W-CT91-0087. Zhang, H. and F. Zhuang (1992), "NND Schemes and Their Applications to Numerical Simulation of Two- and Three-Dimensional Flows," Advances in Applied Mechanics, 29, 193-256.
0306-4549/93 $6.00+0.00 Copyright © 1993 Pergamon Press Ltd
NUMERICAL STUDIES OF TRANSPORT IN POROUS MEDIA
R. L. BUCKLEY and S. K. LOYALKA Mechanical and Aerospace Engineering and Particulate Systems Research Center University of Missouri-Columbia, Columbia, MO 65211 (USA) and M. M. R. WILLIAMS Electrowatt Engineering Services (UK) Ltd. Horsham, Sussex, United Kingdom
(Received September 8, 1992) ABSTRACT
Williams has recently developed a new formulation for flow through porous media by using a modification of the kinetic equations of particle transport. We explore in this paper several numerical techniques for solution of this equation. The methods explored include the upwind, the Smolarkiewicz method [19831, and the NND schemes of Zhang and Zhuang [1992]. The numerical methods were benchmarked against an exact twodirectional solution (the forward-backward model) and it was found that of these finite difference methods, the NND schemes, developed to accommodate shock waves, reduce the inherent diffusion error the most while avoiding spurious oscillations. The numerical techniques were further explored for a multi-directional case for which exact solutions cannot be obtained, and it was again found that the NND schemes outperformed the other methods. INTRODUCTION A porous medium consists of a solid phase and interconnected pores (or voids). These voids, which may be either liquid or gaseous in nature, generally have no distinct pattern (i.e.the geometric and hydraulic properties are spatially variable, or heterogeneous), which makes the problem of describing the transport in a porous medium a tedious one. Porous media transport depends on many factors: the physico-chemical interactions between the solid material and solutes (void-fluid), the often unsteady and non-uniform flow conditions, the complicated geohydrological structures of the aquifer, as well as the mechanisms by which the solute spreads [Dagan, 1987]. The phenomenon of transport of flow through a porous medium is relevant in numerous areas. "Groundwater," for example, a term used to describe the flow of water through the Earth's crust, has long been an area of intense study. The movement of pollutants through ground water is of increasing importance in considering alternatives for disposing of nuclear waste. In a more general sense, those working in such fields as hydrology, chemical engineering, environmental engineering, and soil physics [Dagan, 1987] are concerned with tracking the motion of the fluid contained within the voids in an accurate fashion. This knowledge is important in many applications, such as the groundwater application above, where 2O3
204
R . L . BUCKLEYet al.
pollutants are created as a byproduct of the fission reaction, as well as the accurate knowledge of deposition of pollutants into the surrounding environment (whether it be through complicated urban building structures, or the intricate branching of vegetation, for instance). Due to the natural intricacies found on the microscopic scale, engineers and scientists are often concerned with the average values of the flow or transport variables. In other words, the material is regarded as a spatially random medium, which requires stochastic analyses. Such statistical analyses require a simplification of the equations describing the flow phenomenon. The starting point for analytically describing the complicated porous media flow of a species is the advective-dispersion equation [Williams, 1992a; Bear, 1979]:
0C & =
-V. J - tC
(1)
J = - D V C + UC
(2)
where
with C {= C(r, t)} denoting concentration, D {= D(r)} being the dispersivity tensor, U {= U(r) } representing flow velocity, and ;Ldescribing the decay term. This equation gives the transport due to overall advective motion as well as that resulting from the variable paths taken by the fluid flowing through the porous medium. It is assumed that the motion due to diffusion obeys Hck's Law. The dispersion tensor has been determined experimentally, and from mathematical models, to take the form [Bear, 1979]:
Dij =aTSijU +(aL--aT)
U
(3)
where aT and aL are the transverse and longitudinal dispersion lengths, respectively (which are functions of the porous medium geometry), and Ui are the orthogonal components of the flow velocity. Now, if the dispersion lengths are to be considered fundamental parameters of an experimental system, their values should be independent of measurement size or time. However, in the experimental determination of D, it has been found that at, is dependent on the distance for which measurements are taken (i.e. the value of aL depends on the "size" of the experiment) [Williams, 1992b]. This indicates that dispersivity increases in proportion to the distance of particle flow travel and the prescription (3) is not adequate for description of flow in both small and large scale porous media. It is therefore of considerable interest to investigate eqn (1) without use of the prescriptions (2) and (3) Often, a stochastic (random) averaging technique based on the continuum concept is used to solve eqn (I). One obtains a mean value equation which describes the scale dependence of aL. While this approach provides a solution that is non-Fickian, it is also expensive to implement in practical situations. Obtaining random properties for the medium is difficult, and the results acquired from these methods are commonly used only as indicators of general trends.
Numerical studies of lranspor! in porous media
205
Williams has thus proposed a new model based on the discrete kinetic theory of particle transport and not on continuum concepts used in dispersion theory. The stochasticity of the process is incorporated by assurmng the particles are subject to random collisions. While the model equation can be solved analytically for simple cases, in general it would be necessary to employ numerical techniques to solve the equation. Thus, our goal in this research is to explore various numerical techniques to solve the new equation and to examine the strengths and weaknesses of these methods. NEW FORMULATION AND ANALYTICAL SOLUTIONS The new formulation [Williams, 1992a, 1992b] is based on an equation for the species concentration dependent on space, direction, velocity, and time (r, [1, v, and t, respectively). After various simplifications to the transport equation, [Williams, 1992a], the following integro-differential equation is obtained for the 'angular concentration' C(r, ~ , O:
[R~+vfl.V+vZ(r,~)+)lJ~]C(r, fl, t) (4)
+1
= f(r, fJ)fdff
~(r, f f )C(r, f f ,t)+ S(r,O,t)
-1
where S represents a source term and R is the retardation factor. The new concepts involve f and ~. The scattering function, ~, represents a reaction rate. In essence, it is the chance that a particle can be taken out of the flow media (such as by deposition onto the solid surface, or via a chemical reaction). The anisotropy function,f, on the other hand, represents a redistribution of the flow field. It answers the question of where the flow travels when a branch or fracture intersection is encountered. Thus, in onedimensional transport, we have: [~+
d vlt--~ + v ~,(X,la ) + Z ]C(x,lz, t) = vf (x,lt )!÷tdlt" ~.(x,.')C(x,lt',t) + S(x,lt,t) (5)
where/~ = cos 0 and 0 is the angle that the particle makes with respect to the x-axis as it v flows through the structure. The velocity v is now ~. Williams has made substantial progress in obtaining analytical solutions to the various approximations of eqn (5), but for realistic geological situations, numerical methods are required. We begin by assuming ~ a n d f a r e spatially constant. A complete description of a porous medium would involve an anisotropy function of the form [Williams, 1992a] K
f(/z)= Y~akt0z-/~D k=l
with
(6)
206
R . L . BUCKLEYet al. K
~a t = 1
(7)
k=l
where k refers to the direction taken by the particle and 'a' with this probability, 0 < a < 1. We also assume the source term is:
(8)
S(x,//,t) = So(X,t)f(gt) The resulting coupled first order partial-differential equations are:
( t9
+ v//k
/
'
°~ +vz~(//t)+/q, Ct(x,t)=vat~-r,(//j)Cj(x,t)+So(x,t)at. Tx :=1
where k = 1,2 .... ,J. By assigning / 3 t = v / / t , g / = v , ~ ( / / / ) + ~ , , and at) = vat ~(//j), we obtain a more compact form:
(9)
St(x,t)=atSo(x,t)
'
~+flt
+~'t Ct(x,t)= Y.°tkjCj(x,t)+St(x,t)
(10)
j=l subject to the initial and boundary conditions
Ck(X,O) = gk(X) such thatx ;~ 0, and Ck(0,/) = hk(t), where the functions g and h are known. An analytical solution can be obtained for eqn (10) if all cz are set to 0, which removes any coupling of the equations describing flow in various directions. We choose to examine a special case of eqn (6) which retains the directional coupling. Referred to as the forward-backward model, this scheme is commonly used in neutron transport and involves only two directions,//x and//2. As the name would imply, this simple scheme assumes the species flow is restrained along a crack inclined at an angle 01 (between 0 and 90 °) relative to the fluid flow direction, with the possibility of back-flow in the opposite direction at an angle of 02 (90 ° < 02 < 180°). This leads to the anisotropy function being written as [Williams, 1992b]: f ( / / ) = a S ( / / - / / t ) + (1 -
a)6(//+//2).
(11)
We seek a solution of the form:
Cr(x,t) = C~(x,t) + C2(x,t).
(12)
Assuming a source term as in eqn (8), Williams [1992b] has shown that one can obtain the coupled set of equations:
o3 + (I - a)~q+~]C1(x,t)_aZ2C2(x,t )= aSo (x, t) [~o3+//~-~
(13a)
Numericalstudies of transport in porous media d B [~-laz-~+a~2+~.]Ca(x,t)-(1-a)~,~Cl(x,t)=(1-a)S,(x,t)
207 (13b)
where it should be noted that ~1 = ,,~(/.tl), .,~z = ~(-/t2). In addition, velocity has been absorbed such that vt--~t, 2/v--*Jt, and vC-~C. These equations can be combined to yield a generalized Telegrapher's equation:
~'--~+t--~-~+ sc-~t]C.(x,t)=[D--~2-U ~x - ~ ]Co(X,t)+Qo(x,t)
(14)
where .St = a,S2+(1- a)2~t +Z, ,~,'~
1
v~, r " t = I~! -1~2
'
Z,r+Z Zr #1#2 v,
and U = ~ { v~--"(., -/22) + a ~2//, - (1 - a).~t/.t2 } • Now, if Co(x, t)= Cl(x, t), then the source term is:
whereas, if
Co(x, t)= C2(x, t), then
Qo(x,t)= r ( a - a ) [ ° + ..~+(/~+~l)}]So(X,t). ,
(16)
Summing eqns (15) and (16), or solving eqn (14) with Co(x, t)-- CT(X, t), yields
Qo(X,t)=[--~r{~-V*ff-~}+l]So(X,t) where
(17)
208
R.L.
BUCKLEY et al.
(1 - a)/~l
V" = a/z 2 -
for the total concentration. Williams has obtained an exact solution to eqn (14) by assuming a spatially infinite medium, x e (-00,-0), a source of the form
(18)
So(X,t)= go<~(x)S(t),
with S(t) = 1, t > 0; 0, otherwise; and using Green's function and Laplace-Fourier transforms. Assuming ~ = ~1 = ~2, and utilizing the following:
I:~1
( ~2
1
1
U.e'~2
r
7
"D
4D2.)J
! ~=
+
+ ~
2x and C=
2rd)+Ul £2 + 4 D ' r '
with
• =xZ
~=vZt
D= vD° Z
U=vUo
~= ~° v£
v~=vq
as well as exp
v
- --
~j£2 + 4DoXo 1
g(x,t)-tt, and
v+Jk v;JJ
l = £_.~0
Z #=VZ#o
V ° = yV;
C
=
V#~,C o
Numerical studiesof transport in porous media
209
1
tk
v; j
k
v2 2J
he obtained for CT(x,t):
{/
v:l
X ~o l + - ~ - + l sr_ l t - - ~•- lIe x p lr-c.~l .-TTl+
v; ) k
'
ls(7_Fo)exp(_Cffo)A(~,F°)d~ ° L v; j ._
v; j
}
if o < .~ < v,%
Cr(x,t )
=
{iv.
_
_,
Z "Co 1----e-°]s(,+x-~--]expl cox1+ ~s(t-t°)exp(-Cffo)fr(~,~°)dE ° if-vj
v;3
v;/k
LV;j
}19,
<:~_< O,
where
fr(X't)=El-'Co(Co+(U°-Co¢)V: 2D° )]Io[fl°g(x,t)] (20) +--'Co[3o 1 P(x,t)+ 1+ 2 v-~j
l~[fl, g(x,t)]
and lo[flog(X,t)] and ll[flog(X,t)] are Bessel functions of the first kind of order zero and s(z) -- 1 if z _>0, and equal to zero, otherwise.
one, respectively. Also,
Extending the analysis for Cl(X,0 and C2(x,t) we obtain:
fs0-io/eXp(-qo4(~.~o,aio}'
[ZJa*o(l+~Llt-=]expr-~-+ "" - • -Co~l+ I[~ ~j~ v~j LV~,~ ifO < ~
5 v~i,
q(x,t) =
Z a'¢o 1-
t+
i f - v 2 i _ < ~_< O, with
exp
+ s(~-~o)eXp(-Co~o)ft(~,~, )
R.L. BUCKLEYet al.
210
(22)
+~a'Coflo((ll~2ip(x.t)+{lk + MVo )P(x,t)J _.'-!J,,,ao~,x.,,, z Lk v'~) and, similarly: tI .~ [(1- a)'ro(1./"Vo) 27.+k- s(t-=]od \ v~,)
]
..-.=1 -Co~
L v~, j
+ (?- t-o)exp(-Co~o)B(~,~o)di-° if0 < ~ <
v~?,
C2(x,t) =
v~ )ik"
v'~)
LV~j (23)
Z
+ ~_~(?-L)exp(-cJ°)f2(X,to)dto
i f - v ~ ? < ~ < 0, with
f2(x,t)=[(1-a)~o(~,+Z,1)-(1-a)~o(Co-(U°-c°£°)~lll 2D° Io[flo(g(x,t)] (24)
+1(1-a)'¢°~°[(a+lal)P(x't)+lx-~-~'l+l-'l-'~]ll[fl°g(x't)]'v-~) ~, v; )P(x,t)J NUMERICAL TECHNIQUES The necessity for developing efficient numerical schemes for the transport equation is evident in th:' analytical solutions are feasible only for special cases. In this section, we examine several numerical techniques based on finite differencing. The numerical results are compared with the exact solutions for the forward-backward model as described by eqns (19), (21), and (23). A. UPWIND METHOD We first describe the simple upwind scheme. For clarity, we will assume v to be constant. When applied to the wave equation,
N u m e r i c a l studies o f t r a n s p o r t in p o r o u s m e d i a
211
(25)
OC(x,t) ~ v °~C(x't) = 0
Ot
dx
one utilizes backward (or forward) spatial differencing and forward temporal differencing to obtain
c~i "÷' -C7 I- v C , -C?_, At
= 0
Ax
C.'_+' - C 7 ~-v C7+' - C7 = 0 At Ax
i f v > 0,
(26a)
if v < 0
(26b)
where i refers to the spatial step and n to the temporal step. Often, the terms sought in these numerical computations are those containing the (n+l) time step, since they represent the unknowns in time, and the time step is generally the outermost loop used in finite difference codes. It is assumed that from initial conditions, all values at the (n) time level are known. Eqn (26) has a truncation error of O[Ax, At] and is thus first order accurate. We call this form explicit since only one term, t,-,,t+l ~ , is unknown at each time step. This equation is stable provided the Courant number, r/=[~-~
(27)
The method suffers from numerical dissipation (amplitude) error which is typical of methods with first order accuracy. However, its extreme simplicity lends itself to easy programming. A variation of the upwind method described has been utilized to model eqn (10) as given below: (t~n+l
r,n
t'i,k.~i,k
n
n
~ , on (C;,k--C;-1,k ~a_an t',n+'
J~ ,,,,,n ¢.-,n+l 'K = Z.~ mi,kj~'i,j
(28a)
if ]],'.k > O, and qk
-qk
.
J •
~,
-.n
t~n+l
+
Sink
(28b)
j=l
if .8;:.
o.
The added subscript k denotes the direction taken by the particle. Note the use of concentration associated with a "new" time step for the decay term and the coupling term. This is a modification to the usual procedure which has been devised to avoid instability in the solution.
212
R.L.
BUCKLEY
et
al.
Due to the presence of the coupling term, we have an implicit scheme based on direction with a system of J equations and J unknowns. For instance, in the forward-backward model, we have two directions. Thus, in eqn (28), J = 2 and we have the coupled set of equations: ("~n+l ~i.l
--
C/~.I'
"4" ~i,1
Ca+l
CT, l _ C ~ n _ i , l
n
At ~1
i.2
~i,2
At
. o.
+ t..'~.2
n n+l "I- Ai,ICI,I --
Ax Cin2
t.-,,~t
~-1,2
'
~ ..,.
n ,.i. Si, 1 _
• ["~n'l'll_#,.iCn (",.+I ai,ll,.,i, 1 _ _i,12~1,2
~_,.,. r.+;
.+i
-
(29a)
c..+l
(29b)
+ ,.,,:,,...:,+_+ = --Sl.~- -+.2~-+.,. + ~i,221-'i,2
Ax
if flu, fl,'a > 0, and similarly for 3~j, 3,'.2 < 0. We thus have two unknowns denoted by the (n + 1) superscripts. Solution to a general set of equations involves using matrix operations. Thus, eqn (28) can be solved numerically using the following matrix representation: [A] (y) = (c)
(30)
where [A] is the matrix of known coefficients, (e) is a vector containing known parameters, and (y) is the vector of unknowns (i.e., the "new" concentrations). In general, these components can be written as: tt
tl
II
1 + Ai: - ai,u At n
--O[i•21
-oci,;2
-ot~,t~
1 n n At + Ai '2 - Cti.22
-Ocl,z~
[A] =
-a~m"
" :
: a
n
--~i,J1
--O~i,J2
c?., + s:. ' At
c?,2 (c) =
+ s .2 - / L + s?,, -
II
-o~i•ls
tl
-O~i'as (31)
1 + A,'.3
At
:
.
.
.
.
~i,3J
.. n
--O[i,j3
..
1 At
n
4"
n
t~i, J -- OCi,jj
c?.,-c_?_,,, ,'i.~
n
...
Ax
n n c?,2 - c i-1,2
(32) C,'.~ - C?_ u
--
if fl > 0, and similarly for non-positive values• Finally, (y) is
Numerical studies of transport in porous media
213
(33) i
"
Lc,-.;'j Gaussian elimination and back substitution are then used to obtain (y). This procedure is carried out first through all the spatial sections, x, and directions,/~, for a given time step. Once this is complete, the concentrations are updated and we proceed to a new time step. Though this modified implicit upwind scheme is more expensive computationally than an explicit method, it is necessary to avoid spurious results (i.e. negative concentrations). For the forward-backward model, we have to make use of the anisotropy functional form given by eqn (11). This is expressed in eqn (28) via the source term, which is assumed to be as given in eqn (8) where
So(x,t ) = -So~(x)S(t). For the numerical solution, we can thus take:
Sf,k =Sola~k,l~,i+(1-a)Sk.2S~_l,i]S(t)
(34)
where d/above is the Kronecker delta. B. SMOLARKIEWICZ
Smolarkiewicz [1983] noted the dissipative problems associated with the upwind model and utilized Taylor series expansions about the point (xi, tn) to closely approximate the hyperbolic wave equation, eqn (25), to obtain: OC(x,t) ~
o3t
o3[v(x,t)C(x,t)] o3x
- o3[ K ~ OC(x,t)] -~
Ox
(35)
where
(IvlAx-
K/,,~ =
vZAt)
(36)
is a diffusivity term which helps to eliminate the errors of the upwind scheme. This is essentially a predictive step. To correct this, he reversed the solution of the differential equation as
3C(x,t) ~- O[v~C(x,t)] = 0 Ot ~x where
(37)
R.L. BUCKLEYet al.
214
{
- K ~ . . BC(x,t) v~= C(x,t) Bx 0
if
C(x,t) > 0
if
C(x,t) = 0
(38)
is the diffusion velocity. In addition, it is convenient to introduce a function: 1 At
+
F(n~,n,÷~,~):~[(~
I~)ni+(~-l~)ni+t].
(39)
In order to use this scheme in our problem, we def'me:
(
'
= )gi,kCi,k -- E ai.kjCi,j + S j=l
At
(40)
where )~, a, and S represent decay, coupling and source terms, respectively. Note the use of "old" concentrations with the various terms. Instability does not arise, thus leading to an explicit scheme which is computationally cheaper. Utilizing the notation describing the upwind scheme, the predictive step can then be expressed as: n
F
n
n
n
n
n
~
n
where fin _1 n + ,I i+~/2,k -- "2 (fli,k fli+l,k) /~,,
1
.
(42a)
n
i-~.k : "~(fli-l.k + fli,k ).
(42b)
Essentially, fli+¼.k represents the velocity at a distance midway between mesh point (i) and (i+1) for a given direction k. After the predictive step, we apply
C i,k n+l = C i *. k -
F
*,k,Ci+l,k, *
,i+y2, k - F
-1.k,Ci,k,
,i_y2, k
+ Zi.k
(43)
to obtain the concentration at a new time step. The forward-difference "anti-diffusion" (or anti-dissipation) velocity in the above is expressed as:
fl~.i+~.k :
(Ci:, *C;+,.k + e)Ax
(44a)
Numerical studies of transport in porous media
215
where • is - 1 0 "15 to ensure one does not encounter a computational division-by-zero error. A similar expression is written for the backward-difference anti-diffusion velocity as:
+
+
(44b)
The maximum Courant number, I/, is given by the usual ratio where the velocity in question is represented by the absolute value of the maximum anti-diffusion velocity and the truncation error is of O(Ax2, At2). C. NND SCHEMES
Zhang and Zhuang [1992] have developed a series of numerical schemes which were originally designed to accommodate the complex flow fields associated with shock waves. Numerically modelling such shocks leads to spurious oscillations near or in the shock region. In the past, modeling attempts have been concerned with utilizing mixed-dissipative f'wstorder schemes containing free parameters near the shock wave, and second-order models elsewhere, with resulting unsatisfactory oscillations. Thus, Zhang and Zhuang have suggested a scheme which is nonoscillatory, contains no free parameters, and is dissipative (as opposed to dispersive, which implies oscillations), which they logically refer to as NND. To illustrate the scheme, consider the exact Navier-Stokes equation governing a steady flow with constant viscosity and containing a normal shock wave: ~'
c3C
. a2C
(45)
= o
subject to:
x ~ -**,C --.>C,,.,-~-- ~ O, where ~ is the kinematic viscosity and "o." refers to freestream conditions. The finitedifference form of this equation can be expressed as: 0C
dt
OC
.32C
~ B~C
. 03C
,. B4C
(46)
where the terms on the right-hand side of the expression are the truncation error associated with the discrete representation of the continuous function. To totally suppress the oscillations associated with numerically modelling eqn (45), Zhang found that the coefficient, ~'3, must be selected such that the Second Law of Thermodynamics is satisfied. A detailed analysis is made with the conclusion that if ~3 > 0 upstream of the shock and ~'3 < 0 downstream of the shock, the oscillations can be suppressed. In physical terms, if ~3 < 0 everywhere, then the increasing entropy condition
216
R.L. BUCKLEYet
al.
is not met in the upstream portion of the shock, and oscillations result in this region. The converse is true for ~3 > 0. Zhang and Zhuang further consider the equation: OC + alp(C) = 0.
at
(47)
ax
If we let p =/3C, then we have:
Oc
+
o~
= 0.
(48)
We define the function p as:
p = p+ + p- =/3+C +/3-C where the (+) quantities denote upstream locations and (-) quantities signify downstream locations with respect to the shock wave. Now,
2
2
which yields:
W+
(49)
+-L-x=0.
The following are then applied in all NND-schemes: (1) Upstream (a) 2nd order upwind difference in place of op+
Jx
(2)
(b) 2rid order central difference in place of 0pOx Downstream (a) 2nd order central difference in place of op+
ax
(b) 2nd order upwind difference in place of Op-
ax"
Placing this in the finite difference form of the equation, we have:
(50) where
Numerical studies of transport in porous media 1
217
P~+l- 1 rain m°d(APi-+~,APi-+~)
(51a)
h/_~ = P+-I + 1 rain mod(Ap+_3~,Ap+_/~)+ p~-- 1 min mod(ApT_~,Ap~ ~)
(51b)
h/+~ = p+ + -) min mod(APL~'AP++~ + 2
with ± - P~:, APi-~ ± = P~ - P~-l, etc... @ i±+ ~ ----Pi+l and
[Ixl if Ixl <
•
mlnmod(x,y)
:ly,I
lY[ x and y of like sign,
iflxl
(52)
> Lvl,x and y of like sign, ifx and y differ in sign.
These choices then ensure compliance with the Second Law of Thermodynamics. The schemes which result from the basic formulation given by eqns (47) to (52) are second order accurate in space (except at extremum points), and are dissipative since ~ < 0. In addition, the authors prove the semi-discrete NND schemes possess the TVD, or total variation diminishing property. Civen rc":Xlq
,,-crl
,then
,:0.
We
l
examine two of these schemes. L
NND-1
We use Euler first order time differencing such that
af. I n = C i,k n+l - C i,k n
(53)
where i denotes the spatial node, k the considered direction, and n the temporal step, as before. Modifications similar to that performed in Smolarkiewicz's method with the use of eqn (40) lead to the NND- 1 scheme being written in finite-difference form as:
Ci,k "+l=C"i , k - - ~At(h,, h" k i+y2'k-- i-Y2'k
)
"
+ ~,i,k
(54)
where h is as given by eqn (51) for a given direction, k. The maximum Courant number, r/, is 2/~ while the truncation error associated with this method is of O(Ax 2, At). Note the use of an explicit form of solution with regard to direction.
218
R.L.
2.
BUCKLEY
et al.
bIND-2:
This is a two-step predictor-corrector scheme. We have, C~,~-C n -
At (h"
-n
~+ n
(55a)
1[ n +C~, t A t ( h . * )] 1 * cn+li,k = "~ Ci, k -- " ~ , i+~/2,k -- hi_~2,k + ~Zi.k
(55b)
where the * quantities are temporary predictor concentration values used in the correetor step. In this instance, 77 = 1 and T.E. = O(Ax2,At2). The second order accuracy in time makes this scheme an improvement over NND-1. R E S U L T S AND D I S C U S S I O N S A. FORWARD-BACKWARD MODEL Simulations for the four numerical schemes presented have been performed and compared to the analytical solution. The parameters selected are as follows: a = 0.9, ~1 --- ~r2 = 1,
[09] [o,1 [091
/al = 1£2 = 0.9, ~, -- 0, and v = 1 leading to fl = -0.9 •
= 0.9 ' S = 0.1 6(x), and
a = [001 009]. That is, eqn (28) is explicitly:
oaC~(x,t) + 0.9 BC~(x,t) + 0.1C~(x, t) = 0.9C2(x,t) + 0.96(x)S(t) Ot Ox
(56a)
o3C2(x, t) _ 0.9 0C2(x't) + 0.9C2(x,t) = 0.1C~ (x,t) + O.18(x)S(t) Ot dx
(56b)
where the source term with regard to//(x) is as described by eqn (34). 1.2
Exact
-
1.0 =
. 0.8-
- .........
Upwind
----
NND-I
- .......
Smolar.
1.5 Exact
~
0.62 o
Upwind
o
0.4-
+
1.0
. . . .
NND-I
........
$molar•
~ 0.5
0.2-
_J
0.0
! -1.0
-0.5
0.0
!
!
0.5
1.0
0.0
' '
-10
I
I
I
I
-5
0
5
l0
X
X
Figure l a : F o r w a r d Direction C o m p a r i s o n at t = 1.0 sec
Figure l b : F o r w a r d Direction C o m p a r i s o n at t = 10.0 sec
s t u d i e s o f t r a n s p o r t in p o r o u s m e d i a
Numerical
219
Results were obtained for the four schemes (upwind, Smolarkiewicz, NND-1, and NND2) at elapsed times of 1.0, 10.0 and 100.0 seconds and compared with the exact results to illustrate the dissipative nature of the numerical schemes with time. Figure 1 shows the nature of the forward solution (C1), while Figure 2 describes the nature of the backward solution (C2). We note that for this particular case both of the NND methods lead to similar results, so the accompanying figures contain only the results of NND- I. To enable a meaningful comparison, we used the same timestep, At, for all methods. (The timestep was varied only with respect to the total elapsed time). We discuss the nature of our results as obtained by each method below. 1. E X A C T
The exact solution has a jump in concentration (resulting from the clef'tuition of the source term) at x = 0 as indicated in Figs. 1 and 2. As expected, the exact solution exhibits sharp breaks at or around the jump for both the forward and backward directions. 1.5 -
F~tswt
° 11 o
~
Upwind
mmm~
NND-I
; 1.o
Smolar.
iI
~ 0.5
J
0.0 - 100
-50
!
|
50
100
Figure lc: Forward Direction Comparison at t = 100.0 sec As the total time increases, one can see that an asymptotic solution is approached. For the
/
0.15
I= o
=
0.10
4=#
= 0.05 o
0.20
Upwind
Upwind . . . .
NND-1
........
Smolar.
0.15 a tw
NND-1
O U 0.05
|
-0.5
.... -.---.-.
". O.lO
0.00
-1.0
Exact
0.0 x
0.5
1.0
Figure 2a: Backward Direction Comparison at t = 1.0 see
ii •
0.00 -10
i -5
i 0 x
i 5
|
10
Figure 2b: Backward Direction Comparison at t -- 1 0 . 0 s e e
forward direction, the maximum concemration occurs after the positive x-region is entered. It is interesting that while for small times the backward direction peaks just before x ffi 0
220
BUCKLEY et al.
R.L.
and then drops off, for large times the maximum once again occurs asymptotically at x just greater than 0 before decreasing at the larger distances. This phenomenon can be explained by the form of the exact solution as time approaches infinite. 2. UPWIND For this particular case, the implicit upwind scheme proved to be computationally the quickest. However, this is true only because our problem involves only two directions. For a more general problem containing 5 directions, for example, the use of Gaussian elimination and backsolving would become increasingly time consuming. Examining the figures, one will notice a 'rounding' of the edges near sharp comers as well as a drop in number concentration at the largest distances when compared to the exact solution. This indicates a dissipative (or diffusive) error resulting from the discretization of the finitedifferencing scheme. Again, such behavior is expected for the simple first order modification of the upwind scheme. 0.15 .........
Exact Upwind
....
NND-1
........
Smolar.
°"
e~
0.10 L., .,J gl
= 0.05 O
0.00
,
-100
-50
,
-J
.
0
50
100
X
Figure 2c: B a c k w a r d Direction C o m p a r i s o n at t -- 100.0 sec
3. SMOLARKIEWICZ This method incorporates an "anti-diffusive" velocity to alleviate dissipative problems. While the dissipation was reduced (in comparison to upwind), unsatisfactory dispersive (oscillatory) errors about the shock resulted. A particular advantage of the Smolarkiewicz scheme over the upwind method is that the former is less sensitive to varying timesteps. In other words, one can use a larger timestep (hence, less computational costs) and retain stability. 4. NND Figures 1 and 2 reveal the same dissipative errors occurring with the NND schemes as in the other numerical methods. However, the results are clearly superior. Not only is the dissipative error reduced, but the oscillations inherent with the Smolarkiewicz scheme are also avoided. Varying the value for the timestep, At, revealed that the NND-2 scheme was less sensitive to changes than NND-1. In fact, it was the least sensitive of all numerical techniques studied here. On the other hand, the NND-1 method was most effected by an alteration of the timestep, which is not surprising considering its need for a smaller Courant number than the other numerical schemes for stability. Although in this example the NND-1 model performed as well as NND-2, it should be noted that if the need for larger timesteps arose,
Numerical studies of transport in porous media
221
NND-1 would be the first to suffer instability in its solution. This is also not surprising since it is only In'st order accurate in time. B. MULTI-DIRECTIONAL TEST CASE A separate test case involving three directions (and no analytical solution) was developed to ascertain the utility of the new formulation and to examine the performance of the numerical methods. In this instance, the particle has the option of travelling at an angle of 10% 80°, or 135 ° with respect to the flow direction. The scattering cross-section, £ ~ ) , has been assumed to equal the absolute value of/~. The distribution function,f(/t) is given by eqn (6) where K = 3. The following values for probability of particle flow direction and scattering cross-section were assigned: Direction 1 2 3
St 0.9848 0.1736 -0.7071
2~(4u) 0.9848 0.1736 0.7071
a 0.4 0.5 0.1
F 0.9848 ] f l = / 0 . 1 7 3 6 1, / / L-o.7071J
With v = 1 and A = O, this corresponds to the numerical values:
7098481
l
70.3939 0.0694
z = / 0 . , 7 3 6 | , and a=10.4924 0.0868 0 ~ , ~ , / Assuming Ck(x,O)= 1 . 0 a n d n o [0.7071J L0.0985 0.0174 0.0707J source introduction, we obtain the results shown in Figs. 3 and 4. In Fig. 3, we show 0.50
1.50
0.40 @
..............
Upwind NND-1
/ A
....
NND-2
~
........ Smolar.
0.30
~!i
1.25 ~
/S
1.00
!~!iNND-I Upwmd NND-2 Smolar.
~0.75~
0.20
0.50"
¢J
0.10
. . . . .
0.00 0.0
0.2
0.4
0.6
;'~~~7
~0.25 ,
0.8
X
Figure 3a: Time = 1.00, Direction 1
].0
0.00
0.0
I
I
0.2
0.4
X
I
I
l
0.6
0.8
1.0
Figure 3b: Time = 1.00, Direction 2
results of the various numerical techniques at the final time of t = 1.00 sec, while in Fig. 4, a comparison of the NND-2 scheme and the simple upwind modification (for all directions) for final elapsed times of 0.50, 1.00, and 1.50 see is illustrated. Several inferences may be deduced from the figures. From Fig. 3, it is evident that the upwind scheme gives results that vary considerably from the other numerical techniques. The more diffusive nature of this scheme is depicted by the 'flatter' concentration curves.
222
R . L . BUCKLEY et al. 0.60 O m I.
~
030 -'
a
r
O.2O
lI O
0.10 -~ 0.00 0.0
,
,
0.2
0.4
X
s
,
,
0.6
0.8
1.0
Figure 3c: Time = 1.00, Direction 3 From Fig. 4, at the smallest elapsed time o f 0.50 sec, we see a peak value near x = 0 for the third direction. This makes sense when one considers the angle for this direction is
.
1~0]S
'-~
1.25
Dk. 1--Upw. --""o---
Dir. I--NND-2 Dir. 2--Upw.
0.75 o
0.50
1"1
it. t¢ ~- ~" I0,
---~'-"
Dir. 2--NND-2
....
Dir. 3--Upw.
- - ~t -
Dir. 3--NND-2
0.25 ~ 0.011
0.0
0.2
0.4
Figure 4a:
x
0.6
0.8
1.0
Three.Directional Variation at Time = 0.50 se¢
greater than 90 ° . With flow going in the opposing direction, particles will tend to get washed outmore quickly with distance. Similarly, for the 1 and 2 directions, the maximum values occur at downstream distances.
f
Dir. l--Upw. Dir. 1--NND-2 Dir. 2--Upw.
0.75
~ - ~ ......
o,o I "',j,-. 0.00
.
.
.
.
----
--
0.0
0.2
0.4
X
0.6
0.8
1.0
Figure 4b: Three-Directional Variation at Time = 1.1111sec
Dir. 2--NND-2 Dh'. 3--Upw. Dir. 3--NND-2
Numerical studies of transport in porous media
223
As the elapsed time increases, we see a gradual drop in concentration in all directions, with a more prominent decline in the In'st and third directions. This is expected since no source is introduced. With regard to directions 1 and 3, their loss is quicker than in 2 since their respective reaction cross-sections, Z, are larger and the probability of the particles taking these directions (via the parameter a) is smaller than for the second direction. As time increases further (not shown), we note a gradual shifting in the second direction's peak toward larger downstream distances. This indicates that the majority of the contaminant is being washed away in the direction of fluid flow. 1.50 oo o~
1.25 O al IN
Dir. l--Opw. Dir. 1--NND-2
1.00
Dir. 2--Upw.
0.75 /£I
o
........ • ......
0.50
Dir. 2--NND-2 Dir. 3--Upw. Dir. 3--NND-2
0.25 0.00 0.0
0.2
0.4
Figure 4e:
X
0.6
0.8
1.0
T h r e e - D i r e c t i o n a l Variation at T i m e = 1.50 sec
CONCLUSIONS We have discussed the importance of accurately modeling flow through porous media in the area of mass transport. The inadequacies of current simplifications to the dispersiveconvective equation have led Williams to develop an equation based not on continuum concepts but rather on the kinetic theory of particle transport. Due to the mathematical complexity of the resulting one-dimensional integro-differential equation, analytical solutions are difficult to produce. However, Williams has devised one for the simple case of flow in only two directions, the forward-backward model. We have explored in this report several numerical schemes which apply to more general physical applications. Comparisons with the analytical solution to the forward-backward case have been performed. Based on these comparisons, it is our assessment that the modified NND schemes are most promising due to their avoidance of oscillations and their characteristically smaller dissipative error. While the upwind scheme was computationally most efficient for the simple two-direction problem, it is apparent that the explicit forms of Smolarkiewicz and NND methods are warranted for more complicated systems involving several directions. A CKNO WLED GMENTS
Research performed by Robert Buckley was under appointment to the Environmental Restoration and Waste Management Fellowship program administered by Oak Ridge Associated Universities for the U. S. Department of Energy.
224
R.L. BUCKLEYet al. REFERENCES Anderson, D. E., Tannehil, J. C. and R. H. Pletcher (1984), Computational Fluid Mechanics and Heat Transfer, Hemisphere Publishing Corporation, 40-3. Bear, J. (1979), Hydraulics of Groundwater, McGraw Hill. Dagan, Gedeon (1987), "Theory of Solute Transport by Groundwater," Annual Review of Fluid Mechanics, 19, 183-215. Smolarkiewicz, P. K. (1983), "A Simple Positive Definite Advective Scheme with Small Implicit Diffusion," Monthly Weather Review, 111, 479-86. Williams, M. M. R. (1992a), "A new model for describing the lxansport of radionuclides through fractured rock, Part I: General Theory, "Report to CEC under contract F12W-CT91-0087. Also Annals of Nuclear Energy 19 (1992). Williams, M. M. R. (1992b), "A new model for describing the transport of radionuclides through fractured rock, Part II: Numerical results and extension to overlapping fracture sets," Report to CEC under contract F12W-CT91-0087. Zhang, H. and F. Zhuang (1992), "NND Schemes and Their Applications to Numerical Simulation of Two- and Three-Dimensional Flows," Advances in Applied Mechanics, 29, 193-256.