Du Lihui, Jiang Chunbo and
Xing Xiuying
Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China
E-mail: jcb@mail.tsinghua.edu.cn
Abstract: This paper presents a numerical model of 3D seepage flow for non-uniform medium. The governing equations of saturated-unsaturated flow are introduced and the finite element method is applied to discretize them in space. The techniques of determining the seepage faces and computing the discharge rate through a specified section are proposed. Therefore, this model describes processes in relatively smaller space and time units. It has been tested and implemented in the research of complex 3D seepage flow. The results reveal that both the seepage flow rate and the level of groundwater are affected by the scale of impervious wall.
Keywords: 3D seepage flow, impervious wall, finite element method, anti-seepage
Towards 3D seepage flow under some complex geological conditions, many problems are still not solved because their governing equations are non-linear differential ones and the complete analysis of 3D numerical models is usually restricted by computational units and methods. The SUPG[1] method is used to solve the problem of 1D unsaturated seepage; the equations of 3D saturated seepage flow are mainly utilized in isotropic and homogeneous seepage flow to assay the coupling action between seepage and stress[2]; a 3D finite element program is applied to analyze 2D steady saturated flow without seepage faces[3]. In this paper, a finite element model of seepage flow is built for non-uniform medium. Compared with boundary element method (BEM) and finite difference method (FDM), finite element method (FEM) can be applied in the irregular outer boundary as well as different medium of the domain.
In order to calculate the discharge rate through a specified section, the finite element formulae for sectional flow are obtained, using Galerkin’s method and considering the nonhomogeneity of grids as well as the influence of different permeable coefficients among different medium. The hydrological processes in the computational domain are described in comparatively smaller units of space and time. With the study on parallel computing methods further developed and parallel computers more sophisticated, it will become possible to analyze huger models of 3D seepage flow, especially for some problems with millions of free dimensions.
The governing equations of 3D seepage flow for non-uniform medium can be written as:
(1)
where
,
,
represent pressure head, porosity of pervious materials and moisture
percentage, respectively.
is equal to
, and
is equal to
.
is a constant. When
is equal to 1, it means that the domain is saturated; when
is equal to 0, it means that the domain is unsaturated. Finally
is a coefficient of permeability, which denotes the anisotropic permeable
character of the earth.
The boundary conditions are as follows:
(1) Pressure head is initialized, e.g. the water level of upstream and downstream at the embankment.
(2)
(2) The boundary of computational domain is assumed to be impervious, that is, the normal derivative of pressure head is equal to zero.
(3)
(3) Rate of flow is initialized, e.g. the rate of flow and rainfall.
(4)
where
,
stands for the outer normal vector and the known rate of flow, respectively.
(4) The solution of the seepage face is an iterative process. For the specified point at the natural boundary, it can be considered as a point at the seepage face if the following expression is fulfilled.
,
(5)
where
represents the discharge rate
through the boundary and is greater than zero when the flow direction is
outward. It is supposed that the sum of the pressure head of all these points is
equal to the potential head. Then the distribution of the pressure field is
analysed again.
Combined with the boundary conditions above, the finite element method is utilized to solve the governing differential equations. Finally the distribution of pressure head and velocity in the computational domain can be obtained.
The
discharge rate through a specified section
in space should be calculated to
test the water-sealing effect of the impermeable wall. It is assumed that, one
surface of a tetrahedral element
, is located in the specified section
. Then the formula for the discharge rate Q can be written as follows:
(6)
where
represents the total number of the
tetrahedral elements, of which one surface is located in the same side of the
section
. The numerical result reveals that, under the condition of uniform grids and
medium, the computational precision of expression (6) is relatively high.
However, under the condition of complex terrain and non-uniform medium, the
finite element grids are non-uniform, therefore the precision is low. The
following are two main reasons:(1) The
tetrahedral elements are linear, therefore the first-order derivative is a
constant in each element. Because the grids are non-uniform, the velocity in the
elements is not equal to the actual velocity through the section
. (2) Because the coefficient of permeability in different mediums is very
different, the calculation of the discharge rate is influenced by the selection
of the elements. Finally this paper presents a high-precise formula for the
discharge rate through specified section in space.
(7)
where
(
) represents the sum of the tetrahedral elements, which include the nodes of the
specified section
and are in the same side of
. And
stands for the number of the
surface of the elements.
The seepage
flow of a well, whose radius
is 75cm and water depth
is 200cm, has been investigated to
verify the numerical model built above. The computational domain is a
cylindrical region, whose outer radius
is 600cm and height
is 500cm. The coefficient of
permeability
is equal to 10-5cm/s.
The theoretical formula for the seepage discharge into the well from the medium
around is as follows:
(8)
The finite element model can be utilized to analyse the characteristics of flow in the region. In use of equation (7), the seepage discharge through the cylindrical section is about 3.144cm3/s. And in use of equation (8), the result is about 3.1734cm3/s. The difference is less than 1%. Therefore, the numerical model proves to be right.
To further evaluate the performance of numerical model presented in the previous sections, the author has utilized it to simulate 3D seepage flow around a dam foundation. The terrene composition of the foundation is shown in Fig. 1.

Fig. 1 Geological composition at the foundation of a dam
Along the river, there are three high permeable layers, whose obliquity is about 35 degrees and width is about 110m from the center of the riverbed. Their coefficient of permeability is about 20 to 50 Lu (1Lu is equal to 10–5cm/s). Around them are some low permeable layers, whose coefficient of permeability is about 2 to 5 Lu. What can be seen from the figure is that, the geological condition is very complex and the grids of the elements are non-uniform. The impervious wall is built in the section across the axis of the dam. The terrene composition of the section is shown in Fig. 2(a).

(a) (b)

(c) (d)
Fig. 2 Construction schemes of impervious wall in different operating conditions
In the section, there are three horizontal zones, which are all high permeable layers. The construction scheme without impervious wall is defined as operating condition 0. When the thickness of impervious wall is 4m, according to different width, the sectional figures are shown in Figure 2(b), 2(c) and 2(d), respectively. Their construction schemes are defined as operating condition 1, 2 and3, respectively. When the thickness is 2m and the width is 1.5H, the construction scheme is defined as operating condition 4. The counting parameters under different operating conditions are shown in Table 1.
Table 1 Scopes of impervious wall in different operating conditions
|
Operating condition |
0 |
1 |
2 |
3 |
4 |
|
Thickness (m) |
4 |
4 |
4 |
4 |
2 |
|
Width (m) |
0H |
0.5H |
1.0H |
1.5H |
1.5H |
|
Area (m2) |
|
2570 |
4250 |
6700 |
6700 |
Then the discharge rate
through different zones across the axis of the dam is obtained, as shown in
Table 2. The relevant positions for the calculation of flow are shown in Fig. 3.
Table 2 Amount of leakage in different zones of the section (L/min)
|
Operating condition |
0 |
1 |
2 |
3 |
4 |
|
Q4 |
36.6 |
37.2 |
37.7 |
38.8 |
37.8 |
|
Q3 |
43.5 |
52.0 |
57.2 |
71.7 |
65.7 |
|
Q2 |
31.3 |
33.9 |
53.0 |
19.0 |
21.3 |
|
Q1 |
118.5 |
93.3 |
60.9 |
72.5 |
87.7 |
|
Sum |
229.9 |
216.4 |
208.8 |
202.0 |
212.5 |

Fig. 3 Positions for the calculation of flow in the section
From table 2, some conclusion can be drawn as follows:
(1) The amount of leakage decreases with the scope of the impervious wall increasing.
(2) As the influence of the water permeability of the materials in different zones, the amount of leakage Q4 increases with the scope of the impervious wall increasing. The reason is that, after the zones 1-3 are dealt with seepage prevention, the permeability in these zones decreases largely and then the relative permeability in zone 4 increases.
(3). The anti-seepage effect is corresponded with the quantity of anti-seepage material. In operating condition 4, the width is 1.5H and the thickness is 2m. The amount of leakage through the section across the axis of the dam is about 212.5L/min and the total volume of anti-seepage material is about 14400m3. In operating condition 1 and 2, the average amount of leakage is about 212.6L/min and the total volume is about 14180m3. Therefore, if the volume of anti-seepage is equal, the anti-seepage effect is alike.
The 3D seepage model built in this paper can be applied in the research of anisotropic underground water flow in heterogeneous medium. In this model, the problems of saturated and unsaturated flow are combined to be solved, therefore the spatial position of free surface can be defined quickly and conveniently. In this paper, the author develops an iterative solver for the problem of 3D seepage flow and presents his own opinion about the solution of discharge rate through specified section. The model is used to solve the complex 3D seepage problem and calculate the amount of leakage. The relation between the scope of the impervious wall and anti-seepage effect is discussed.
Acknowledgements
This research was supported by the State key 973-Project for fundamental Research and Development Program(Grant No. G1999043607 ) and the National Science Foundation (Grant No. 59979013).
[1] Zhou Xueyu, Xie Chunhong, Qian Xiaoxing, SUPG finite element method of unsaturated flow. Journal of Hydraulics, 1996, 37-42.
[2] Du Yanling, Research of inconstant seepage governing equations. Journal of Hydraulics, 1993, 62-68.
[3] Griffiths D V , Fenton G A. Three-dimensional seepage through spatially random soil. J of Geotech Engineering, 1997, 153-160.
[4] Li L , Barry D A , Pattiaratchi C B. Numerical modeling of tidal-induced beach water table fluctuations. Coastal Engineering, 1997, 30: 105-123.
[5] Lee K K , Leap D I. Simulation of a free surface using boundary-fitted coordinate system method. J Hydrology, 1997, 196:297-309.
[6] Masters I , Usmani A S , Cross J T , et al. Finite element analysis of solidification using object-oriented and parallel techniques. Int J for Numerical Method in Engineering, 1997, 40: 2891-2909.