 Забыли?

?

# A Fourier-Laplace transform finite element method (FLTFEM) for the analysis of contaminant transport in porous media

код для вставкиСкачать
```INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
A FOURIER-LAPLACE TRANSFORM FINITE ELEMENT
METHOD (FLTFEM) FOR THE ANALYSIS OF
CONTAMINANT TRANSPORT IN POROUS MEDIA
J. C. WANG* AND J. R. BOOKER
Department of Civil Engineering, University of Sydney, NSW 2006, Australia
SUMMARY
In many practical situations, a simpli"ed two-dimensional model is often su\$ciently accurate for practical
purposes, however on some occasions it is still necessary to consider the full three-dimensional situation.
Conventional numerical techniques such as the "nite di!erence and "nite element methods are able to deal
with contaminant transport under 3D conditions, however the computational e!ort involved in using them
is often very large. Under some circumstances it is possible to perform 3D analyses more e\$ciently.
A particular case is when the geometry is 2D but the problem is actually 3D. In such cases the analysis can be
simpli"ed by application of a Fourier transform.
In this paper, the dimensions of the contaminant transport problem will be reduced by one and the
transient mass transport equation will be transformed into a time-independent mass transport equation
through the use of a combination of integral transforms (Laplace and Fourier transforms) and a "nite
element technique. The mathematical model developed here is called Fourier}Laplace transform "nite
element method (FLTFEM). By using the FLTFEM a three-dimensional transient problem can be analysed
as a two-dimensional time-independent problem and the computational cost will thus greatly be reduced.
Copyright ( 1999 John Wiley & Sons, Ltd.
KEY WORDS: contaminant transport; Fourier transform; Laplace transform
INTRODUCTION
In many important practical situations, a simpli"ed two-dimensional model is often su\$ciently
accurate for practical purposes; however on some occasions it is still necessary to consider the full
three-dimensional situation. Conventional numerical techniques such as the "nite di!erence and
"nite element methods are able to deal with contaminant transport under 3D conditions;
however the computational e!ort involved in using them is often very large. Under some
circumstances it is possible to perform 3D analyses more e\$ciently. A particular case is when the
geometry is 2D but the problem is actually 3D. In such cases the analysis can be simpli"ed by
application of a Fourier transforms.
Although the combined use of integral transforms and approximate numerical methods to
solve linear di!erential equations can be seen in the literature,1}7 it appears most of them are only
based on the application of a Laplace transform. Booker and Leo8 were probably the "rst to use
* Correspondence to: J. C. Wang, Department of Civil Engineering, University of Sydney, Sydney, NSW 2006, Australia.
E-mail: jc.wang@civil.usyd.edu.au
CCC 0363}9061/99/141763}34\$17.50
Copyright ( 1999 John Wiley & Sons, Ltd.
Revised 26 October 1998
1764
J. C. WANG AND J. R. BOOKER
a combination of both Laplace transform and Fourier transform and boundary element methods
to study the e!ects of contaminant migration from waste repositories.
In this paper, the dimensions of the contaminant transport problem will be reduced by one and
the transient mass transport equation will be transformed into a time-independent mass transport equation through the use of a combination of integral transforms (Laplace and Fourier
transforms) and a "nite element technique. The mathematical model developed here is called
Fourier}Laplace transform "nite element method (FLTFEM). By using the FLTFEM, a threedimensional transient problem can be analysed as a two-dimensional time-independent problem
and the computational cost will thus greatly be reduced.
The numerical method (FLTFEM) developed here may be applied to assess the impact of
contamination transport (see Figure 1) due to a surface waste disposal repository or a prismatic
waste repository which the geometry being considered is the same along one spatial dimension.
One restriction of FLTFEM is that the material properties must not vary in one spatial
dimension, although they may vary in the remaining two spatial dimensions. For example, in the
case of contamination due to a crack on the bottom of a long channel (Figure 1(b)), the geometry
and the material properties are the same along the y direction, however there may be layered soils
with di!erent properties in the z direction (Figure 1(b)) or the x direction.
Figure 1. Contamination problems due to a prismatic waste repository and a surface land"ll. (a) A surface land"ll;
(b) a case of a prismatic waste repository
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1765
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
In the FLTFEM, a Gauss}Lengendre quadrature formula is used to invert the Fourier
transform while in the case of the Laplace transform, inversion is performed using a convenient
technique developed by Stehfest.9, 10
2. DEVELOPMENT OF FLTFEM FOR CONTAMINANT TRANSPORT PROBLEMS
2.1. Governing equations
The theory of a regular two- or three-dimensional fracture network within fractured media
developed by Rowe and Booker11, 12 o!ers a theoretical base to apply "nite element techniques
to analyse the contaminant migration in fractured media. Hence, the general form of the
governing mathematical equation of contaminant transport in nonfractured media and in
fractured media based on the theory of Rowe and Booker11,12 can be written as
\$ ' (D ' \$c)!V ' \$c"n
!
!
Lc
#g#q
Lt
(1)
where c is the concentration of the contaminant in the pore water of non-fractured media or the
concentration in the fractures of fractured media, D is the apparent hydrodynamic dispersion
!
tensor, V is the vector of the components of the apparent groundwater velocity (Darcy velocity),
!
\$ is the gradient operator, n is the porosity of the soil in the non-fractured media or the fracture
porsity in the fractured media, g is the sum of the rates at which contaminant is lost from the
groundwater due to sorption onto the solid matrix and the e!ects of radioactive decay and
biodegradation and q is the rate per unit volume at which contaminant is being transported into
the matrix material by matrix di!usion in the fractured media or is omitted in the non-fractured
media. The source term q in the Laplace domain has been established by Rowe and Booker11}13
and is given in the appendix.
If the principle axes of hydrodynamic dispersion are chosen to be parallel to the co-ordinates
axes, the contaminant transport equation for three-dimensional homogeneous non-fractured or
fractured media under steady groundwater #ow condition is given by the following equation:
A
B
L2c
L2c
L2c
Lc
Lc
Lc
Lc
D
#D
#D
!<
!<
!<
"nR
#j*c #q
!xx Lx2
!yy Ly2
!zz Lz2
!x Lx
!y Ly
!z Lz
Lz
(2)
where c is the concentration of the contaminant in the pore water of non-fractured media or the
concentration in the fractures of fractured media, n is the porosity of the soil in the non-fractured
media or the fracture porosity in the fractured media, R is the retardation factor
(R"R "1#oK /n) in the non-fractured media or R"R "1#b (K /n ), in the fractured
\$
\$
&
& & &
media in which b , K , and n are de"ned by Rowe and Booker),11, 12 j*, j*"j/R , in which j is
& &
&
\$
the decay constant, in the non-fractured media or j*"j/R in the fractured media, < , < , <
&
!x !y !z
are the components of the apparent groundwater velocity (Darcy velocity) in the x, y, and
z directions and D , D , D are the apparent hydrodynamic dispersion coe\$cients.
!xx !yy !zz
2.2. Transformations of the governing equations
Both Fourier and Laplace transformations are utilized in the Fourier}Laplace transform "nite
element method (FLTFEM). The Fourier integral transform is "rst applied to the governing
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1766
J. C. WANG AND J. R. BOOKER
equations in order to make a reduction of the dimensionality of the contaminant transport
problem by one. In applying such a transform, the assumption is made that the material
properties of the soil and geometry do not vary in the y direction although they may change in the
xz plane and then a Fourier transform
=
P~= c (x, y, z, t) e~*uy dy
1 =
c"
C (x, u, z, t)`*yu du
2n P
~=
C"
(3)
(4)
of the governing equation (2) becomes
A
B
L2C
L2C
LC
LC
LC
D
!(u2 D #iu< ) C#D
!<
!<
"nR
#j* C #Q
!xx Lx2
!yy
!y
!zz Lz2
!x Lx
!z Lz
Lt
(5)
where C is the Fourier transform of concentration (c), Q is the Fourier transform of q and u is the
Fourier transform parameter.
Following the application of a Fourier integral transformation, a Laplace transform is
introduced:
=
P 0 C(x, u, z, t) e~st dt
CM "
(6)
Then the contaminant transport equation in Fourier}Laplace transform space can be written as
D
L2CM
L2 CM
LCM
L CM
#D
!<
#<
"(u2 D #iu< #h) CM !nRC
!xx Lx2
!zz Lz2
!x Lx
!z Lz
!yy
!y
0
(7)
where s is the parameter of Laplace transform, CM is the concentration in Fourier}Laplace
transform space, C is the initial concentration in Fourier transforms space,
0
=
P~= c0 e~*uy dy
C "
0
G
h"
n R (s#j*)
in a non-fractured material
\$
n R (s#j*)#g6 (s#j* ) in a fractured material
.
f
&
in which g6 and j are given in the appendix.
.
It is obvious that the transformed problem in the Fourier}Laplace transform space becomes
2D problem provided that the boundary conditions can be expressed in terms of C. This will
always be the case if the boundary conditions are of the same type along any line x"constant
and z"constant.
For the special case of 3D contaminant migration from a surface land"ll through a horizontally
layered soil media in which the soil properties are the same at any horizontal location within the
layer, if a second Fourier transform is taken along the x-axis, then the method reduces to the "nite
layered method.14
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1767
2.3. Boundary and initial conditions
In order to completely de"ne the problem, it is necessary to specify the boundary and initial
conditions. Although these conditions depend on the particular problems being treated, the ones
which occur often will be described here. More complicated situations such as the speci"c #ux
boundary condition treated in the section on leakage will be discussed in detail when the related
problem is introduced.
2.3.1. Specixed concentration at boundary. As previously discussed the FLTFEM can be applied
when the boundary conditions are of the same type and the material properties do not vary in one
spatial dimension. A particular case of a boundary condition which is often considered in this
paper is introduced here. Considering a surface land"ll in which the leachate in the land"ll has
a constant concentration (c ) and zero concentration is assumed elsewhere. The schematic
B
diagram of this situation is shown in Figure 2. The boundary condition of this surface land"ll may
be written as
G
c ,
c" B
0,
D y D(H/2, D x D(¸/2, z"0
D y D'H/2, D x D'¸/2, z"0
Introducing a Fourier transform, the transformed surface boundary condition can be expressed
as
=
P~=
C "
B
H@2
P~H@2
c e~*uy dy"
A B
2c
uH
c e~*uy dy" B sin
B
u
2
(8)
Suppose that the concentration c is constant with time. It follows that taking a Laplace
B
transform of the boundary conditions in Fourier}Laplace transform space leads to
G
2c
B sin
su
CM "
B
0,
A B
uH
,
2
D x D(¸/2, z"0
D x D'¸/2, z"0
(9)
where c is the constant concentration on the boundary in real space and CM is the transformed
B
B
concentration in Fourier}Laplace transform space.
Figure 2. Boundary condition of a surface land"ll
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1768
J. C. WANG AND J. R. BOOKER
Figure 3. Schematic of an initially contaminated block
If the material properties of the soil do not vary in the y direction, then the above transformed
boundary conditions can be applied in conjunction with the FLTFEM to determine the
contaminant migration from a 3D land"ll by solving a 2D time-independent transport (7).
2.3.2. Initial concentration in the contaminated zone. Considering the case of a block of contaminated soil (Figure 3) in which the initial concentration (c ) is assumed to be spatially constant
0
within the contaminated zone and zero initial concentration elsewhere; then the initial condition
may be written as
G
c , D y D(H/2, D x D(¸ /2, D z D(¸ /2
x
z
c" 0
i
0
D y D'H/2, D x D'¸ /2, D z D'¸ /2
x
z
Taking the Fourier transform of this relation it is found that
G
2c
:H@2 c e~*uy dy" 0 sin (uH/2) when D x D(¸ /2; D z D(¸ /2
x
z
u
~H@2 0
CM "
0
0
when D x D'¸ /2; D z D'¸ /2
x
z
(10)
where c is the initial concentration in real space and C is the initial concentration in Fourier
0
0
transform space.
If the material properties of the soil do not vary in the y direction and the boundary conditions
are of the same type in the y direction, then the above transformed initial conditions can be
incorporated with equation (7) in the FLTFEM to examine the extent of contamination due to
this initially contaminated soil
2.4. The xnite element scheme in transformed space
The "nite element approximation of the transformed contaminant transport equation (7)
together with its associated boundary and initial conditions of particular problems can be
obtained by using the Galerkin "nite element technique. By using the method of weighted
residuals in the "nite element scheme, the equation for the contaminant transport problem can be
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1769
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
expressed as
PP)
C
= (x, z) D
L2 CI
L2 CI
L CI
L CI
#D
!<
!<
!xx Lx2
!zz Lz2
!x Lx
!z Lz
D
!(u2D #iu< #h) CI #n R C d)"0
!yy
!y
0
(11)
where = (x, z) are the weighting functions (i.e. interpolation functions) and ) represents the
domain of a typical section. If CI denotes the approximate solution in Fourier}Laplace transform
space, the approximate solution can be written as
m n
(12)
CI (x, z, u, s)" + + Ne (x, z) CM e (u, s)
i
i
e/1 i/1
where Ne (x, z) are interpolation functions, CM e (u, s) are the unknown nodal values of the "eld
i
i
variable in Fourier}Laplace transform space, n is the number of nodes within the element, m is the
number of elements, and the superscript e indicates a particular element.
Using integration by parts and substituting equation (12) into equation (11), the intergral in
equation (11) can be evaluated and leads to a system of linear equations which may be cast in
matrix form as
[M] MC1 N"MF N#MF N
(13)
*
&
where [M] is the global matrix including advection}dispersion, dispersion, sorption, decay
matrices which may be obtained by combining the element matrices and is given by equation (14),
[MC1 N is a vector containing the unknown nodal concentrations in Fourier}Laplace transform
space, MF N is a vector that depends on the known initial conditions which is given by equation
*
(15), and MF N is a vector due to the prescribed #ux ( f ) which is the speci"ed outward normal #ux
&
n
along the boundaries is given by equation (16). If no speci"ed #ux is encountered, the #ux is zero
and this term will vanish:
m
[M]" +
e/1
P
m
!+
e/1
m
#+
e/1
C DC DC D
C DC DC D
L Ne LNe
1
1
Lx
Lz
F
F
L Ne LNe
n
n
Lz
Lx
A%
P
L Ne LNe
1
1
Lx
Lz
F
F
L Ne LNe
n
n
Lx
Lz
A%
P
Ae
CD
Copyright ( 1999 John Wiley & Sons, Ltd.
De
!xx
0
ve
!x
0
De
!zz
0
L Ne
LNe
1 2
n
Lx
Lz
dA
L Ne
LNe
1 2
n
Lz
Lz
Ne 2 Ne
n
1
dA
0
ve
!z
Ne
1
2 Ne
n
Ne
1
F [he#De u2#iu<e ] [Ne 2 Ne] dA
!yy
!y
1
n
Ne
n
(14)
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1770
J. C. WANG AND J. R. BOOKER
m
MF N" +
*
e/1
P
Ae
m
MF N" +
&
e/1
CD
PC D
Ne
1
F [ne Re Ce ] dA
0
Ne
n
(15)
Ne
1
F [FM ] d!
n
Ne
n
(16)
!e
The above equations describe the "nite element approach for approximating the equation
governing contaminant migration by using the FLTFEM. The initial condition for a contaminant transport problem is incorporated into the transport equation (7), when a Laplace
transform is introduced and is then incorporated into the analysis through equations (13) and
(15).
The concentration is also often speci"ed on a portion of the boundary, this is easily incorporated into a "nite element analysis by specifying the transformed concentration at those boundary
nodes. Alternatively, the normal #ux entering the boundary may be speci"ed. It is then necessary
to transform the speci"ed normal #ux, this will be described in detail later. The condition of
continuity of #ux and concentration at boundaries separating di!erent materials is implicit in the
"nite element procedure.
2.5. Numerical inversion methods
The solutions obtained directly from the "nite element method are in the Fourier}Laplace
transform space. Hence, numerical inversion techniques are necessary to invert those solutions in
transformed space into the solutions in real time and real space. The inversion of the Laplace
transform can be performed using a convenient algorithm developed by Stehfest9,10 while in
the case of Fourier transform inversion, a conventional numerical integration technique,
The inversion of Laplace transform has been given in detail in References 9, 10 and 15. Hence,
only the inversion of Fourier transform will be discussed here. The concentration of contaminant
(c) is related to the transformed concentration (C) by equation (4) which is restated here for each
of reference:
1
c"
2n
=
P~= C (x, u, z, t) e`*yu du
(17)
The integral in equation (17) may be evaluated numerically by the use of Gauss quadrature as
follows. First the in"nite integral is truncated and then broken into a number of subintervals so
that
1
c"
2n
=
P~=
C (x, u, z, t) e`*yu du :
1 k/`p~1
+
C (x, u, z, t) e`*yu du
2n
K/~p
(18)
where ¸ is the length of the subinterval.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1771
Each of the integration on these subintervals is then approximated using Gaussian quadrature
so that
m
C (x, u, z, t) e`*yu du : + = C (x, u , z, t) e`*yuK
(19)
K
K
k¸
K/1
where = are the Gauss weights, u are the Gauss co-ordinates and m is the total number of
K
K
Gauss points.
Details of Gauss}Lengendre quadrature may be found in References 16}18. Generally, using
more Gauss points will produce more accurate results. In order to assure the accuracy of
numerical integration technique the numerical solutions obtained in this paper was corroborated
by using large integration range and more Gauss point. Experience in dealing with the contaminant migration problems considered here shows that both the range of integration and the
number of Gauss points required to achieve a speci"ed accuracy of solution depend on the
problem being considered and the time at which the solution is being evaluated. Generally,
a larger range of integration and number of Gauss points will be required when the solution is
evaluated at larger times.
P
(k#1)¸
3. VERIFICATION
The technique (FLTFEM) developed in this paper was veri"ed by evaluating four test examples
including two 2D and two 3D cases of contaminant migration in the non-fractured media. In all
test examples the soil is assumed to have the same geometry and material properties in the
y direction and thus the Fourier transforms can be taken along the y-axis to reduce the
dimensionality of problems being treated. As mentioned earlier, in the FLTFEM,
a Gauss}Lengendre quadrature is used to invert the Fourier transform and the Stehfest inversion
method is used to invert the Laplace transform. The numerical solutions obtained by FLTFEM
were compared with the analytical solutions derived by Cleary and Ungs,19 Leij and Dane,20 and
Booker.21
3.1. Two-dimensional contaminant transport problems
The FLTFEM is "rst applied to contaminant migration under two-dimensional conditions. As
was mentioned previously, through the application of FLTFEM, these 2D contaminant transport problems can be solved by using a 1D time-independent transformed mass transport
equation in Fourier}Laplace space.
3.1.1. Example for constant concentration surface source. Suppose that a surface land"ll (Figure
4(a)) of a "nite width of ¸ is the y direction but with a relatively larger length in the x direction
overlays a homogeneous, isotropic porous medium. Under these conditions, there will be no
variation of concentration in the x direction. A undirectional steady #ow is also assumed to have
a vertical seepage velocity < in the z direction and the values of the dispersion coe\$cients in that
;
direction and orthogonal to it are respectively denoted as D and D . In this example, the
z
y
contaminant is assumed to be non-reactive and non-decaying (R "1 and j"0). The concentra\$
tion of the contaminant is assumed to be spatially and temporally constant in the land"ll and
initially zero elsewhere.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1772
J. C. WANG AND J. R. BOOKER
Figure 4. Schematic diagram for 2D example of constant concentration surface source. (a) 2D problem; (b) Transformed
1D problem
In order to be speci"c, the width of the land"ll source (¸) is assumed to be 100 m. The
transport parameters D "1 m2/day, D "0.1 m2/day, and < "0.1 m/day are used. Since
z
y
z
the geometry and material properties of soil are assumed to be the same in the y direction, the
Fourier transform can be taken along the y-axis in the numerical approach and this twodimensional contaminant transport problem is then reduced to a problem involving only one
spatial dimension (see Figure 4(b)). A one-dimensional "nite element mesh consisted of 111 linear
bar elements and each element of length of 2 m was used to analyse this two-dimensional
problem.
The transformed concentration, which was given by equation (9) in Section 2.3, is speci"ed at
the boundary node of the ground surface because a constant concentration (c ) is assumed at the
B
surface boundary. A no-#ux condition was assumed at the lower boundary in the FEM. However,
the remote boundary (222 m from the ground surface) was found to be deep enough to avoid any
noticeable boundary e!ects.
The numerical results of FLTFEM were compared with the analytical solution of Cleary and
Ungs.19 By using FLTFEM, the concentration of a contaminant at particular positions and at
speci"ed times can be quickly evaluated. In this case it was found that 60 G points, one
subinterval and ¸"4 (the range of integration of each subinterval used in Gauss}Lengendre
quadrature) was adequate to obtain an accurate solution because the solutions examined were
found to be almost the same when even more Gauss points, a larger value of ¸ or more
subintervals were used. Comparisons of FLTFEM numerical results and the analytical solutions
are shown in Figure 5. It can be seen that numerical results of concentration distributions are
nearly identical to those obtained analytically.
3.1.2. Example for an initially contaminated rectangular zone. The purpose of this example is to
test the ability of the FLTFEM to model contaminant transport from an initially contaminated
region. For this reason, groundwater contamination due to waste liquid seeping from the bottom
of a waste repository into a homogeneous aquifer is considered. It is assumed that this gives rise
to an initially contaminated area of a spatially constant concentration in the aquifer and that
there is no migration above and below the aquifer (see Figure 6).
An uniform steady groundwater #ow with seepage velocity < "1 m/y is also assumed in the
x
direction of the x-axis and the transport parameters are taken to be D "10 m2/y and
x
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1773
Figure 5. Two-dimensional examples of constant concentration surface source
Figure 6. A schematic showing a vertical section of the aquifer and the surface impoundment
D "1 m2/y. Both the length and the width of the initially contaminated area are assumed to be
y
10 m. A plan of the details are shown schematically in Figure 7(a).
Since the material properties of the aquifer are assumed to be constant in the y direction the
FLTFEM converts the 2D governing equation of contaminant transport into a 1D transport
equation. Figure 7(b) shows the one-dimensional "nite element mesh of 240 linear bar elements,
which have an equal length of 1 m, used to analyse this problem.
A mesh of "nite extent has to be used since it is not possible to model an in"nite boundary, no
#ux boundaries are assumed at the boundaries of "nite element domain. However the 1D "nite
element mesh (total length of 240 m) used to analyse this problem was found to be su\$ciently
large to avoid any noticeable boundary e!ects. The treatment of the initial condition is incorporated into the analysis through equations (7), (10), (13) and (15).
In this case, 60 Gauss points and two subintervals with ¸"4 (the range of integration of each
subinterval) were found to be adequate in the Gauss}Lengendre quadrature since no signi"cant
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1774
J. C. WANG AND J. R. BOOKER
Figure 7. Schematic diagram for 2D example of initial concentration rectangle source. (a) 2D problem; (b) Transformed
1D problem
Figure 8. Two-dimensional example of initial concentration rectangle source
di!erence was found in the solutions even when more Gauss points, a larger value of ¸ or more
subintervals were used. Figure 8 shows the comparison of the numerical solutions obtained by the
FLTFEM and the analytical solutions from Booker.21 It can be seen that the numerical results
obtained using the FLTFEM provide a very good "t to the analytical solutions as well.
3.2. Three-dimensional contaminant transport problems
The accuracy of the FLTFEM has been demonstrated by the accuracy of the numerical
solutions it provided to the 2D contaminant transport test examples. On some occasions,
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1775
however, it is necessary to consider contaminant migration under three-dimensional condition.
Therefore, it is desirable to extend the veri"cation of the FLTFEM to some 3D contaminant
transport problems. As previously discussed, by applying the Fourier transform together with the
Laplace transform, a transformed 2D time-independent mass transport equation can be used to
solve some particular three-dimensional transient contaminant migration problems. This leads to
signi"cant savings in computational cost.
3.2.1 3D case of a landxll with a leachate which has a constant concentration. In order to further
verify the FLTFEM under 3D conditions, a waste disposal impoundment (land"ll) constructed
overlaying a homogeneous soil, shown in Figure 9, is considered. It is assumed that the waste is
continuously dumped into the land"ll and that the concentration within the land"ll remains
constant over the period of interest. It is also assumed that the groundwater #ow is in the vertical
direction and has a constant seepage velocity < . This 3D contaminant transport problem has
z
been solved analytically by Leij and Dane.20
The contaminant is also assumed to be non-reactive and non-decaying (R "1 and j"0). The
\$
Fourier transform is taken along the y-axis as it was in the previous 2D problems. Thus this
three-dimensional consideration of contaminant transport from a constant concentration land"ll
source can then be reduced to a problem involving only two spatial dimensions (shown in
Figure 10).
Figure 9. A surface land"ll
Figure 10. Schematic diagrams of 3D case of a land"ll with constant concentration condition and the transformed 2D
problem. (a) 3D real problem; (b) Transformed 2D problem
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1776
J. C. WANG AND J. R. BOOKER
The transformed 2D problems is symmetric about the z-axis. Therefore, a simplifed half plane
(xz) domain (Figure 10(b)) was used to analyse this problem and a no-#ux boundary was assumed
as the z-axis (x"0 m) because the problem is symmetric about it. It was assumed that there was
a constant concentration (c ) in the land"ll because the land"ll contained a very large amount of
B
contaminant. Zero concentration was also assumed on the ground surface since it was assumed
that periodic rainfall #ushed any surface contamination away. Since it is not possible to model an
in"nite boundary a mesh of "nite extent was used. A zero #ux was assumed at the other
boundaries of the "nite element domain. In order to be speci"c, the transport parameters D "1
z
m2/day, D "D "0.1 m2/day, and < "0.1 m/day were used in this 3D test example and both
x
y
z
¸ and ¸ were assumed to be 10 m.
x
y
In the "nite element approach, a two-dimensional rectangular "nite element mesh with total
number of 300 linear quadrilateral elements (Figure 11) was used. The mesh was re"ned at the
edge of the land"ll in order to model the high concentration gradients which occur there. The
transformed concentration (CM ) which was discussed in Section 2.3 was speci"ed at the boundary
B
nodes of land"ll (z"0) and is found to be
G
A B
2c
u¸
B sin
y
su
2
CM "
B
0
when
0(x(5 m
elsewhere
No perceptible boundary e!ects were found during the period of interest indicating that the
positioning of the boundary was satisfactory. In this case, the solutions obtained by using 60
G points, two subintervals and ¸"4 (the range of integration of each subinterval) in the
Gauss}Lengendre quadrature were found to be satisfactory because the solutions obtained were
virtually the same when more Gauss points, a larger value of ¸ or more subintervals were used.
Comparisons of the numerical result obtained by using the FLTFEM and the analytical solution
are shown in Figures 12(a)}14(a). Excellent agreement is observed between the numerical and
Figure 11. The "nite element mesh used in FLTFEM for 3D constant concentration surface source
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1777
Figure 12. Three-dimensional example of constant concentration surface source at y"5 m plane. (a) Concentration
distribution curves at di!erent values of x; (b) concentration contours at y"5 m plane (time"100 days)
analytic solutions at di!erent values of y. Figures 12(b)}14(b) also show the concentration
contours of half of the xz plane but at di!erent values of y. Inspecting these "gures, it can be seen
that the concentration pro"les at diferent values of x are obvious di!erent. This is because of
land"ll was considered to be of "nite length in the x direction. A signi"cant di!erence is also
observed at the same value of x (x"4 or 5 m) but di!erent values of y (i.e. y"5 m in Figure 12,
y"0 m in Figure 13 and y"6 m in Figure 14) since the land"ll considered also has a "nite length
in the y direction.
3.2.2. Case of a block of initially contaminated soil. The next example considers the extent of
contamination due to a three-dimensional block of contaminated soil which was assumed to be at
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1778
J. C. WANG AND J. R. BOOKER
Figure 13. Three-dimensional example of constant concentration surface source at y"0 m plane. (a) Concentration
distribution curves at di!erent values of x; (b) concentration contours at y"0 m plane (time"300 days)
a certain depth (h) beneath the ground surface in a deep homogeneous soil (see Figure 15).
A spatially uniform initial concentration was assumed within this contaminated block and the
initial concentration is assumed to be zero elsewhere. It is also assumed that there is zero
concentration on the surface of ground because periodic rainfall #ushes any surface contamination away.
The dimension of the contaminated soil and the speci"c values of parameters used are shown in
Figure 15. A uniform steady groundwater #ow of seepage velocity < "10 m/y was assumed in
x
the x direction and the values of the dispersion coe\$cients adopted were D "10 m2/y and
xx
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1779
Figure 14. Three-dimensional example of constant concentration surface source at y"6 m plane. (a) Concentration
distribution curves at di!erent values of x; (b) concentration contours at y"6 m plane (time"300 days)
D "D "1 m2/y, respectively. The length (¸ ) and width (¸ ) of this initially contaminated soil
yy
zz
x
y
are 10 m and the height (¸ ) of it is assumed to be 5 m.
z
By introducing a Laplace transform and a Fourier transform which was taken along the
y direction, this 3D problem was then transformed into a simpler 2D time-independent problem
which is shown schematically in Figure 16. A two-dimensional "nite element mesh shown in
Figure 17 with a total number of 528 linear quadrilateral elements was used. The nodal
concentrations are zero at the surface and a no-#ux condition was used at the other boundaries of
"nite element domain.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1780
J. C. WANG AND J. R. BOOKER
Figure 15. Schematic diagram of 3D case of a block of contaminated soil with initial concentration condition
Figure 16. 2D domain solved in FLTFEM for a 3D problem of a block of initially contaminated soil
The treatment of initial condition has been discussed previously and the transformed initial
concentration (C ) in Fourier transform space, is given below:
0
2c
u¸
0 sin
y
when !¸ /2(x(¸ /2,!¸ /2(z(¸ /2
x
x
z
z
2
C " u
0
0
elsewhere
G
A B
This condition is incorporated into the transport equation (7) when a Laplace transform is
introduced and is then incorporated into the analysis through equations (13) and (15).
Numerical results shows that using 60 G points, three subintervals and ¸"4 (the range of
integration of each subinterval) are su\$cient to obtain an accurate solution in this case considering that nearly identical solutions were obtained when more Gauss points, a larger value of ¸ or
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1781
Figure 17. The "nite element mesh used in FLTFEM for a 3D problem of a block of initially contaminated soil
Figure 18. Three-dimensional problem of a block of initially contaminated soil
more subintervals were used. It is also found that there is no noticeable e!ect due to the zero-#ux
boundaries used. Comparisons of numerical results and analytical solutions from Booker21 are
given in Figures 18}20. It can be seen that at di!erent times and di!erent values of y the
simulation results obtained by the FLTFEM are in very close agreement with those obtained
analytically.
4. LEAKAGE FROM A WASTE REPOSITORY
Often, leachate migration from sanitary land"lls is an important source of groundwater pollution.
As well, other types of land disposal units including domestic septic tanks, waste disposal through
wells and deeply buried were repository, etc., are potential sources of contamination. Thus, it is
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1782
J. C. WANG AND J. R. BOOKER
Figure 19. Three-dimensional problem of a block of initially contaminated soil
Figure 20. Three-dimensional problem of a block of initially contaminated soil
Figure 21. Schematic diagram of leakage problem in a waste repository
necessary to be able to predict the likely impact of the contamination due to such buried water
repository and to assess if remediation action is required.
Suppose leakage is observed from a long deeply buried cylindrical pipeline (shown schematically in Figure 21) which was constructed for conveying chemical liquid or waste substances. The
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1783
long cylindrical pipe is assumed to be buried deep in the non-fractured porous media. It is
assumed that the pipe develops a break of width (H) at one of the pipe joints and that a constant
rate of contaminant mass #ux ( f ) is released from the circumference of the pipeline into the
0
surrounding soil. Since the pipe is assumed to be very long and that the material properties of soil
do not vary in one spatial dimension (the y direction in this case) thus the FLTFEM can be
applied to assess the impact of contamination due to this leakage.
Two cases of leakage from this cylindrical pipeline are investigated in this section. The "rst one
is the case of leaking of the contaminant from the pipe occurring over a long period (¹ PR).
\$
The second is the case of leakage which occurs over a period of "nite duration (¹ ). The governing
\$
equations and boundary condition for the problem are described below.
The governing equations: The governing equation is equation (2). As in the cases considered
previously, the assumptions that the same geometry and material properties along the y-axis are
made and then, by introducing the Fourier and Laplace transforms, the transport equation in
transformed space is given by equation (7).
Boundary condition at the section of leakage. The mass #ux of contaminant #owing out from the
cylindrical pipeline is assumed to be temporally and spatially constant during the period (¹ ) of
\$
leaking. In other words, there is a constant rate of #ux ( f ), for the duration ¹ , over the section
0
\$
(H) where the leakage occurs (see Figure 21). This boundary condition can be expressed as:
G
f"
f
0
0
D y D(H/2, r"R and if 0(t(¹ over the section of the leakage
\$
Dy D'H/2, r"R
where r is the radial distance from the centre of the pipe (r"Jx2#z2), and R is the radius of the
pipe.
The contaminant #ux on the circumference of the pipe, in Fourier}Laplace transform space,
then becomes
A B
2f
uH
FM " 0 sin
(1!e~sT\$)
su
2
(20)
where FM is the radial contaminant mass #ux in Fourier}Laplace transform space, f is the radial
0
contaminant mass #ux in real space, s is the Laplace transform parameter, u is the Fourier
transform parameter, and ¹ is the duration of leakage.
\$
If the leaking of the contaminant from the pipe occurs over a long period (¹ PR), the
\$
contaminant #ux emanating from the circumference of the pipe in Fourier}Laplace transform
space is given by
A B
2f
uH
FM " 0 sin
su
2
(21)
4.1. Case 1: leakage occurred over an inxnite duration
In this case, it is assumed that a long cylindrical pipeline having a diameter of 2 m is
constructed in an idealized isotropic clayey soil and that is develops a break of 0)2 m wide.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1784
J. C. WANG AND J. R. BOOKER
A constant rate of contaminant mass #ux ( f "1000 g/m2/a) is released into the surrounding soil
0
over the section where the leakage occurs. The lacking of contaminant is assumed to occur over
an in"nite duration (¹ "R).
\$
Also it is assumed that there is a uniform advection velocity (< "!0)02 m/a) in the
x
x direction and the hydrodynamic dispersion coe\$cient is 0)05 m2/a in all directions because the
groundwater velocity in clayey soil is relatively small. The porosity of soil is assumed to be 0)4 and
the contaminant is conservative (K "0) and does not decay (j"0).
\$
The material properties of soil and the geometry considered do not vary along the y direction
and thus a Fourier transform can be taken along the y-axis. This enables the three-dimensional
problem of leakage from a long cylindrical pipe to be transformed into a 2D time-independent
problem which can be evaluated using a 2D "nite element mesh.
Figure 22 shows the actual 3D situation and the 2D cross-section of the domain. A half-plane
2D "nite element domain was used in the FLTFEM. The 2D problem is symmetric about the
x-axis and thus a non-#ux condition was assumed at the boundary z"0 m. The 2D "nite element
mesh of 440 linear quadrilateral elements shown in Figure 23 was used to analyse this problem.
Figure 22. Schematic diagram of 3D leakage problem and its 2D domain used in the FLTFEM
Figure 23. The 2D "nite element mesh used in FLTFEM for 3D leakage problem from a cylindrical waste repository
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1785
A constant #ux ( f ) was released into soil from the circumference of the cylindrical pipe. Thus, it is
0
necessary to transform the constant #ux into the nodal #uxes through equations (16) and (21)
these nodal #uxes are applied at the boundary nodes of pipe. A no #ux boundary condition was
also adopted at the other boundaries in this case since it is not possible to model an in"nite
domain. The simulation results were examined to ensure that there was no noticeable boundary
e!ect due to the "nite size of the mesh adopted. In this case a satisfactory solution was obtained
by using 100 G points, 1 subinterval and ¸"100 (the range of integration of each subinterval) in
the FLTFEM since no signi"cant di!erence in the solution was found by using a larger value of ¸,
more Gauss points or more subintervals.
The concentration distributions in the radial direction at di!erent angles for di!erent values of
y are shown in Figures 24 and 25. The semi-analytical solutions for this problem developed by
Booker and Leo8 are also presented in the "gures. It can be seen the numerical solutions are
almost identical to those of semi-analytical solutions.
Inspection of the simulation results (Figures 24 and 25) shows that the highest peak concentration is found in the direction of h"03 (the upstream side of pipe) and followed by the h"903 (the
top side of pipe) at the two cross-sections (y"0 and 0.1 m). The smallest peak concentration is
Figure 24. Solutions of leakage problem in a waste repository for in"nite leaking duration (¹ "R) at y"0 m plane
\$
Figure 25. Solutions of leakage problem in a waste repository for "nite leaking duration (¹ "R) at y"0.1 m plane
\$
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1786
J. C. WANG AND J. R. BOOKER
found along the direction of h"1803 (the downstream side of pipe) which has the same direction
as that (negative x) of groundwater #ow. This is because the migration of the contaminant in the
opposite direction of the #ow (along the direction of h"03) was retarded by the groundwater
#ow and the contaminant is taken away along the direction of h"1803 from the pipe. Thus, the
contaminant is transported further away from the waste repository in the direction of h"1803
(the downstream side of pipe) than in the direction of h"03 (the upsteam side of pipe). No
matter what direction (h) is considered, the peak concentration is always found at the boundary
of pipe.
4.2. Case 2: leakage occurred over a period of xnite duration
Now, exactly same leakage problem as discussed above (see also Figure 22) was again
considered, except that the leaking of contaminant is now assumed to occur at a constant rate of
#ux ( f "1000 g/m2/a) but over a "nite period of time (¹ "10 y). The transport parameters
0
\$
adopted in this case are identical to those used in the previous case and the "nite element mesh
shown in Figure 23 was also used.
As would be expected, the simulation results at 10 y time are found to be exactly same as the
results obtained in the previous case because during this period (10 y) the two situations are
identical. The solution was evaluated at both 12 and 20 y which are some time after the leaking
ceased and both numerical and semi-analytic results are given in Figures 26 and 27.
Similar behaviour to that found in the previous case is also observed in di!erent radial
directions. The highest peak concentration is found in the direction of h"03 (the upstream side of
pipe) and followed by the h"903 (the top side of pipe). The smallest peak concentration is found
along the direction of h"1803 (the downstream side of pipe) which has the same direction as that
of groundwater #ow. It is also observed that the contaminant is transported further away in the
direction of h"1803 (the downstream side of pipe) and the migration of the contaminant in the
opposite direction of the #ow is retarded by the groundwater #ow.
After an additional 10 y (t"20 y), Figures 26 and 27 show that the concentration of the
contaminant becomes smaller as time increases and the peak concentration in the direction of
h"1803 (the downstream side of pipe) is found not to be at the boundary of pipe. This is because
Figure 26. Solutions of leakage problem in a waste repository for "nite leaking duration (¹ "10 yr) at y"0 m plane
\$
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1787
Figure 27. Solutions of leakage problem in a waste repository for "nite leaking duration (¹ "10 yr) at y"0.1 m plane
\$
leaking of the contaminant ceased 10 y previously and after that time the concentration at the
boundary of the pipe starts to decrease since contaminant is transported away from the pipe and
is not replenished by the leak.
Based on the results presented, it is evident that good agreement is obtained between FLTFEM
solutions and the semi-analytical solutions of Booker and Leo8. It has been demonstrated by this
case that the FLTFEM can also be used to analyse a contaminant transport problem which
involves a change of the environmental conditions by using the multi-stage approach.22
5. 3D WASTE REPOSITORY IN FRACTURED MEDIA
The migration of contaminant in non-fractured media under two- or three-dimensional situation
has been discussed in previous sections. Sometimes the migration of a contaminant occurs in
a "ssured porous media. For example, if a waste disposal repository was constructed overlaying
a "ssured media, it may cause contamination of the "ssured material beneath it sometime after its
construction because of aging, design shortcomings or acts of nature. Therefore, a 3D waste
disposal repository placed over a fractured material, shown schematically in Figure 28, is
examined in this section by using FLTFEM.
Although, in many important practical situations such as land"ll, a simpli"ed 2D or 1D model
is often su\$ciently accurate for practical purposes because the longitudinal dimension of the
contaminant source under consideration is often relatively large compared with the other
dimensions or the depth of penetration of contaminant from a land"ll into the soil beneath it over
the period of interest. It is also interesting to know what the di!erence will be between a 3D
analysis and a simpli"ed 2D or 1D analysis of a surface waste repository.
It is assumed that waste is dumped into a waste repository, shown in Figure 28(a), continuously
and that this leads to a spatially uniform and temporally constant concentration of contaminant
within the repository. A constant concentration (c ) was thus assumed in this waste repository.
B
Initially, it is assumed that there is no contamination of the "ssured material. Also zero
concentation was assumed on the ground surface because it was assumed that periodic rainfall
#ushed any surface contamination away.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1788
J. C. WANG AND J. R. BOOKER
Figure 28. Schematic diagrams of a waste repository under 3D and 2D considerations. (a) 3D surface waste disposal
repository; (b) 3D real problem; (c) 2D problem
Figure 29. 2D fracture system
The fracture system begin considered is shown in Figure 29 and the material properties of
fractured media are listed below.
two sets of fractures:
set 1 and 3
fracture spacings:
H "H "0)05 m
1
3
h "h "0)001 m
1
3
< "n < "0 m/a
!x
pxx x
< "n < "0 m/a
!y
pyy y
< "n < "0)01 m/a
!z
pzz z
fracture openings:
Darcy velocity:
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
Dispersivity:
a "1 m, a "0)1 m
L
T
fracture di!usion coe\$cient:
D "0)01 m2/a
o
matrix di!usion coe\$cient:
D "0)0036 m2/a
.
matrix porosity:
n "0)1
.
1789
The 3D model of this waste repository is shown in Figure 28(b). A simpli"ed 2D model in the xz
plane shown in Figure 28(c) and an idealized 1D model were also used to analyse this problem.
The FLTFEM was used to implement the 3D analysis while the numerical model (LTFEM)
developed by Booker and Wang23 for the fractured porous media was used to analyse the
pollutant migration under 2D and 1D conditions.
The problem domain (see Figure 28(c), and 2D "nite element mesh (see Figure 30) which has
1800 linear quadrilateral elements were used in LTFEM for the analysis of contaminant
migration under 2D condition and a 1D mesh consisted of 100 linear bar elements each of length
of 1 m was used for the analysis under 1D condition. The material properties do not vary in the
y direction, so by applying Fourier transform along the y-axis, the 3D contamination problem
can be treated as a 2D problem. Thus, the exactly same problem domain (see Figure 28(c)) and the
2D "nite element mesh (see Figure 30) can be used in FLTFEM for the analysis of contaminant
migration under 3D conditions.
As mentioned above, a contant concentration was assumed in the waste repository and zero
concentration was assumed on the ground surface. A no-#ux condition was assumed at the other
boundaries of the "nite element domain. However, when the numerical results were examined to
perceptible boundary e!ects were found. The simulation results show that 60 G points, "ve
subintervals and ¸"4 (the range of integration of each subinterval) are satisfactory to obtain an
accurate solution because no signi"cant di!erence was found in the solution by using more Gauss
points, a larger value of ¸, or more subintervals.
Figure 30. Finite element mesh
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1790
J. C. WANG AND J. R. BOOKER
Figure 31. Concentration contours at time"50 yr of 1D analysis
Figure 32. Concentration contours at time"50 yr of 2D analysis
The concentation contours at time"50 years for the analysis using both 2D and 1D models
are presented in Figures 31 and 32. Also Figures 33 and 34 show the concentration contours
obtained for the analysis using the 3D model at di!erent locations of y"0 m, (the centre of the
waste repository) and y"5 m (the edge of the waste repository).
Figures 35 and 36 show the concentration contours at time"100 y obtained by using both 2D
and 1D models. Also the concentration contours obtained by the use of the 3D model at di!erent
locations of y"0 m (the centre of the waste repository) and y"5 m (the edge of the waste
repository) are presented in Figures 37 and 38.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1791
Figure 33. Concentration contours at y"0 m (¹"50 yr) of 3D analysis
Figure 34. Concentration contours at y"5 m (¹"50 yr) of 3D analysis
Examining Figures 31}33 shows that the contaminant plume at the plane of y"0 m obtained
by using the 3D analysis is very close to the results obtained by using the 2D analysis and the
simulation results of the contaminant migration along the depth by the use of 1D model is also
close to those of the analyses of using 2D and 3D models at the position near the centre of the
waste repository. This indicates that the extent of the contaminant beneath the centre line
(y"0 m) of the waste repository may be analysed by using a simpli"ed 2D model or even using
a simple 1D model without sacri"cing the accuracy in this case at 50 yr.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1792
J. C. WANG AND J. R. BOOKER
Figure 35. Concentration contours at time"100 yr of 1D analysis
Figure 36. Concentration contours at time"100 yr of 2D analysis
However, examination of the contaminant plume on the plane of y"5 m (the edge of the waste
repository) shows that there are signi"cant three-dimensional e!ects. The edge e!ects can be seen
when comparing results of the 1D and 2D analyses (see Figures 31 and 32). As would be expected,
investigation of the simulation results (shown in Figures 35}38) at a larger time (100 y) of these
three models shows that the di!erence between a 3D model and a simpli"ed 2D and 1D model
becomes greater at locations further from the centre of the waste repository and the edge e!ects
become more signi"cant as time increases (see Figures 32 and 36 or Figures 33 and 37).
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
1793
Figure 37. Concentration contours at y"0 m (time"100 yr) of 3D analysis
Figure 38. Concentration contours at y"5 m (time"100 yr) of 3D analysis
6. CONCLUSIONS
The Fourier}Laplace "nite element method (FLTFEM) developed in this paper has been
demonstrated to be an accurate method for analysing some particular mulit-dimensional contaminant migration problems such as those involving a surface waste disposal repository and
a long prismatic waste repository in which the geometry and material properties do not vary
along one spatial dimension. The analysis uses a Fourier transform to reduce the dimensions of
the multi-demensional contaminant transport problem by one, while using Laplace transform to
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1794
J. C. WANG AND J. R. BOOKER
simplify a transient contaminant transport equation to a time-independent di!erential equation
involved only spatial derivatives. The "nite element scheme is used to approximate the transformed governing equations. Once the solutions in Fourier}Laplace transform space are obtained, the solutions in real space and time can be attained through the numerical inversion
techniques. The inversion of Laplace transform can be performed using a convenient algorithm
developed by Stehfest 9,10 while in the case of the inversion of Fourier transform, a conventional
numerical integation technique, Grass}Lengendre is used.
With the combination use of Fourier and Laplace transforms, the Fourier}Laplace transform
"nite element method (FLTFEM) greatly simpli"es a 3D contaminant transport problem and
this certainly results in a computational bene"t.
APPENDIX: EXPRESSION OF THE LAPLACE TRANSFORM OF q AND g FUNCTION
The Laplace transform of the function g is used in the theory of Rowe and Booker11,12 to
simulate the di!usive transport of contaminant from the fractures into the solid matrix in the
fractured media. The expression for g6 and qN , which is the Laplace transform of the rate (q) at
which the contaminant migrates from the fractures into the solid matrix per unit volume of
fracture-matrix system, are brie#y introduced below.
Considering that the molecular di!usion is the dominating process in the solid matrix and the
advection is neglected, it is found that the concentration (c ) in the matrix satis"ed the equation
.
Lc
n D +2c "(n #o K ) .#n j c
. .&
.
.
. . Lt
m . .
(22)
where c is concentration in the solid matrix, D is the di!usion coe\$cient of the solid matrix, n
.
.&
.
is the porosity of the solid matrix, o is the dry density of the solid matrix, K is the distribution
.
.
coe\$cient of the solid matrix, and j is the sum of the radioactive decay and biological
.
degradation constants in the solid matrix.
The concentration of the surface of the matrix block will be equal to that in the "ssure system
(c ), so
&
G
c on the surface of the solid matrix
c " &
.
0 initially
By taking Laplace tansform, equation (22) becomes
n D +2 cN "n R (s#j* ) cN
. .
. .&
.
. .
(23)
where R "1#o K /n is the retardation coe\$cient of the solid matrix, j* "j /R and
.
. .
.
. . .
cN "cN on the surface of the matrix.
.
&
Now, the Laplace transform of q may be de"ned as
: FM
dS
.
qN "! ./
: d<
.
(24)
where !F denotes the normal component of #ux entering the matrix, : FM dS is the integral
./
./ .
along the surface S of the fracture, and : d< is an integral over the volume < of the solid
.
.
.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1795
ANALYSIS OF CONTAMINANT TRANSPORT IN TRANSPORT MEDIA
matrix. Using the divergence theorem, it is found
: n D +2 cN d<
n R (s#j* ) : cN d<
. ." . .
. . .
qN " . .&
(25)
: d<
: d<
.
.
It was shown by Rowe and Booker11,12 that the Laplace transform of q can always be expressed
in the form
(26)
qN "(s#j* ) g6 cN
.
&
Referring to the equations (A.4) and (A.5), it is found that the Laplace transform of g is given by
: cN d<
. .
g6"n R
(27)
. . : cN d<
& .
Case 1: A single set of xssures: Assuming there is only one set of "ssures (set 1) and considering
di!usive transport into the matrix in the x direction, the concentration within the matrix satisi"es
the equation
L2cN
."n R (s#j* ) cN
n D
. .& Lx2
. .
. .
(28)
where
cN "cN
.
&
x"\$H
1
It it thus found that
cosh kx
cN "cN
.
& cosh kH
1
(29)
tah kH
1
kH
1
(30)
and
g6"n R
. .
where
k2"R (s#j* )/D .
. .&
.
Also g6 may be expressed as a series expansion,
C
D
1
=
(s#j*)
m
g6"n R 1!2 +
(31)
. .
(s#j* #a2 (D /R ) (a H )2
i
.& . i 1
.
i/1
Case 2. Two sets of xssures: Assuming there are two sets of "ssures (sets 1 and 3) and
considering di!usive transport matrix in the x and z directions, Laplace transform of g is given by
Rowe and Booker.11, 12
C
D
1
1
= =
(s#j* )
.
g6"n R 1!4 + +
. .
(s#j* #(a2#c2) (D /R ) (a H )2 (c H )2
k
.& .
i
.
i 1
k 3
i/1 k/1
where a "(i!1/2)n/H and c "(k!1/2) n/H .
i
1
k
3
Copyright ( 1999 John Wiley & Sons, Ltd.
(32)
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
1796
J. C. WANG AND J. R. BOOKER
Case 3. Three sets of xssures: Finally, if there are three sets of "ssures (sets 1, 2 and 3) and it is
observed that the transport distances of interest along the fracture system are large compared to
the dimensions of the block, then Laplace transform of g is found by Rowe and Booker:11,12
C
D
1
1
1
= = =
(s#j* )
.
g6"n R 1!8 + + +
. .
[s#j* #(a2#b2#c2) (D /R )] (a H )2 (b H )2 (c H )2
.
i
j
k
.& .
i 1
j 2
k 3
i/1 k/1 j/1
(33)
where a "(i!1/2)n/H , c "(k!1/2) n/H and b "( j!1/2)n/H .
i
1 k
3
j
2
REFERENCES
1. J. A. Liggett and L. F. Liu, ¹he Boundary Integral Equation Method for Porous Media Flow, George Allen & Unwin.,
U.K., 1983.
2. H. T. Chen, T. M. Chen and C. K. Chen, &Hybrid Laplace transform/"nite element method for one-dimensional
transient heat conduction problems' Comput. Meth. Appl. Mech. Engng., 63(1), 83}95 (1987).
3. Y. S. Yu, M. Heidari and W. T. Wang, &A semi-analytical method for solving unsteady #ow equations', Paper
presented at Proc. ASCE National Conf. on Hydraulic Eng., Am. Soc. Civ. Eng., Williamsburg, 1987.
4. H. T. Chen and C. K. Chen &Application of hybrid Laplace transform/"nite di!erence method to transient heat
conduction problems', Numer. Heat ¹ransfer, 14, 343}356 (1988).
5. E. A. Sudicky, &The Laplace transform Galerkin technique: a time-continuous "nite element theory and application to
mass transport in groundwater', =ater Resour. Res., 25(8), 1833}1846 (1989).
6. C. K. Chen, T. M. Chen and J. W. Cleaver, &New hybrid laplace transform/"nite element method for two-dimensional
transient heat conduction problem', Numer. Heat ¹ransfer, 20, 191}205 (1991).
7. G. J. Moridis and D. J. Reddel, &The Laplace transform "nite di!erence method for simulation of #ow through porous
media', =ater Resour. Res., 27(8), 1873}1884, 1991.
8. J. R. Booker and C. J. Leo, &A FLT boundary integral equation method for analysis of contaminant leakage', in
Valliappan, Pulmano, Tin-Loi (eds), Computational Mechanics- from Concepts to Computations, Proc. 2nd AsianPaci,c Conf. on Computational Mechanics, Sydney, Australia, 1993.
9. H. Stehfest, &Algorithm 368, Numerical inversion of Laplace transforms', Commun. Assoc. Comput. Mach. 13(1), 47}49
(1970).
10. H. Stehfest, &Remark on Algorithm 368, Numerical inversion of Laplace transforms', Commun. Assoc. Comput. Mach.
13(10), 624}625 (1970).
11. R. K. Rowe and J. R. Booker &A semi-analytic model for contaminant migration in a regular two- or threedimensional fractured network: conservative contaminants', Int. J. Numer. Anal. Meth. Geomech., 13(5), 531}550
(1989).
12. R. K. Rowe and J. R. Booker &Contaminant migration in a two- or three-dimensional fractured network: reactive
contaminants', Int. J. Numer. Anal. Meth. Geomech., 14(6), 401}425 (1990).
13. R. K. Rowe and J. R. Booker, &Modelling of contaminant movement through fractured or jointed media with parallel
fractures', Proc. of 6th int. Conf. on Numerical Methods in Geomechnics, 1988, pp. 855}862.
14. R. K. Rowe and J. R. Booker, &A "nite layer technique for calculating three-dimensional pollutant migration in soil',
Geotechnique, 36(2), 205}214 (1986).
15. J. C. Wang and J. R. Booker, &A Laplace transform "nite element method (LTFEM) for the analysis of contaminant
transport in porous media', Research Report (R756), Department of Civil Engineering, University of Sydney, 1997.
16. M. Abramowitz and I. A. Stegun, &Handbook of Mathematical Function, Dover Publications, New York, 1965.
17. B. Carnahan, H. A. Luther and J. O. Wilkes, 0Applied Numerical Methods', Wiley, New York, 1969.
18. A. Ralston and P. Rabinowitz, A First Cource in Numerical Analysis, 2nd edr, McGraw-Hill, New York, 1978.
19. R. W. Cleary and M. J. Ungs &Groundwater pollution and hydrology, mathematical modles and computer programs',
Rep. 79-=R-15, Water Resour. Program, Princeton Univ., Princeton, NJ, 1978.
20. F. J. Leij and J. H. Dane, &Analytical solutions of the one-dimensional advection equation and two- or threedimensional dispersion equation', =ater Resour. Res., 26(7), 1475}1482 (1990).
21. J. R. Booker, &Analytic and semi analytic methods for the analysis of contaminant migration in soils', Proc. 8th Int.
Conf. on Computer Methods and Adcances in Geomechanics, Morgantown, West Virginia, USA, 1994, pp. 3}19.
22. J. C. Wang, &A study of contaminant transport in porous media', Ph.D. ¹hesis, Department of Civil Engineering,
University of Sydney, 1997.
23. J. R. Booker and J. C. Wang &A Laplace transform "nite element method (LTFEM) for the analysis of contaminant
transport in fractured porous media', Int. J. Numer. Anal. Meth. Geomech., (1999) submitted.
Copyright ( 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Anal. Meth. Geomech., 23, 1763}1796 (1999)
```
###### Документ
Категория
Без категории
Просмотров
6
Размер файла
1 160 Кб
Теги
porous, elements, contaminants, finite, transform, transport, method, fourier, analysis, media, fltfem, laplace
1/--страниц
Пожаловаться на содержимое документа