A new numerical method for the analysis of monolithic seepage problems with complex drainage systems in a groundwater recharge area for a hydropower station in China

Abstract

After construction of a dam impounding water in a reservoir, a monolithic seepage field develops in the surrounding rock mass. Here, a new finite element method is proposed for determining the shape and characteristics of the 3D monolithic seepage field including the free surface, considering complex drainage systems consisting of densely-spaced drainage holes and drainage galleries. To this end, the previously proposed virtual flux method is improved by a refined numerical integration scheme and a regularized Heaviside function for distinguishing the subregions below and above the free surface within a particular finite element. Leakage and overflow drainage holes are modeled as internal boundaries. The proposed numerical method is verified by an academic example, for which the analytical solution is available. Finally, the numerical simulation of the seepage field developing in the vicinity of a high dam and underground power house, constructed in the context of a hydropower plant project in China is used to show its application to a problem in engineering practice.

Headnote

Article history:

Manuscript received April 11,2019

Revised manuscript accepted August 28, 2019

Available online December 20, 2019

Keywords: monolithic seepage problem, drainage system, free surface, finite element method

Abstract

After construction of a dam impounding water in a reservoir, a monolithic seepage field develops in the surrounding rock mass. Here, a new finite element method is proposed for determining the shape and characteristics of the 3D monolithic seepage field including the free surface, considering complex drainage systems consisting of densely-spaced drainage holes and drainage galleries. To this end, the previously proposed virtual flux method is improved by a refined numerical integration scheme and a regularized Heaviside function for distinguishing the subregions below and above the free surface within a particular finite element. Leakage and overflow drainage holes are modeled as internal boundaries. The proposed numerical method is verified by an academic example, for which the analytical solution is available. Finally, the numerical simulation of the seepage field developing in the vicinity of a high dam and underground power house, constructed in the context of a hydropower plant project in China is used to show its application to a problem in engineering practice.

(ProQuest: ... denotes formulae omitted.)

1. INTRODUCTION

In recent years, a large number of high dams, large reservoirs and pumped storage power stations, usually accompanied by underground powerhouses and other caverns, have been (and are) under construction in China. In these construction projects, extremely complex seepage fields develop, especially when complex geological conditions are encountered and/or drainage galleries with arrays of densely-spaced drainage holes are arranged.

In general, the seepage flow in the vicinity of a high dam and an underground powerhouse are analyzed separately. However, after construction of a high dam and impoundment of the reservoir, a three-dimensional monolithic seepage field is formed in the dam, in the vicinity of the underground powerhouse and in the surrounding rock mass. This monolithic seepage field may be approximately regarded as the combination of free surface seepage flow, mainly in the dam body and the surrounding rock mass and seepage flow in the vicinity of galleries with arrays of densely-spaced drainage holes. The behaviour of the monolithic seepage field is directly related to the safe operation of the power station. Hence, the accurate numerical analysis of the complex monolithic seepage field, e.g. by means of the finite element method (FEM), is of great importance.

Seepage flow with a free surface is usually analyzed on the basis of a fixed finite element mesh as originally proposed by Neuman (NEUMAN, 1973), because in contrast to adaptive mesh methods, unfavourable distortions of finite elements are avoided. On the basis of the fixed mesh method, different numeric techniques have been proposed, such as the variable permeability method (BATHE et al., 1979), the residual flow method (DESAI & LI, 1983), the initial flow method (ZHANG et al., 1988), the virtual flux method (VFM) (ZHU & SU, 1991; ZHU, 1997), the variational inequality method (JANSEN et al., 1988), the refined Gauss point method (WANG & HUANG, 1997), the variational inequality method (ZHENG et al., 2005b) and related methods (SHU et al., 2007; ALT, 1980; WANG, 1998; ZHENG et al., 2005a; BORJA & KISHNANI, 1991; GABRIELLA et al., 2016; AZUSA et al., 2018). However, for determining the free surface with higher accuracy, most of the aforementioned numerical methods require a very fine discretization, especially in the vicinity of the free surface, which is computationally expensive for 3D problems in engineering practice.

The FE analysis of seepage flow considering galleries with arrays of densely-spaced drainage holes involves two difficulties: one is the FE mesh generation of arrays of drainage holes with small diameters and dense spacing and the other is the determination of the free surface in the vicinity of the complex drainage system. To overcome these difficulties, either the equivalent permeability method or the direct method can be employed. In the equivalent permeability method, the drainage holes and the surrounding rock mass are homogenized as a complex medium with an equivalent anisotropic permeability (CHEN et al., 2004; JING et al., 2005; CHEN et al., 2010). An advantage of this method is the simple generation of the FE mesh. However, it is not easy to determine the equivalent permeability tensor for the finite elements affected by drainage holes, and thus the drainage facilities cannot be described precisely. In the direct method, the drainage holes are regarded as internal drainage boundaries of the seepage field, which are modeled by a semi-analytical approach (FIPPS et al., 1986) or a substructure technique (ZHU & ZHANG, 1997; CHEN et al., 2008). The direct method allows considering drainage holes with a relatively high precision, but generation of the FE mesh is difficult and it is also computationally expensive.

In this study, an improved virtual flux method (IVFM) for seepage flow with a free surface and a direct method for modeling seepage flow in regions with densely-spaced drainage holes is proposed. After validation of the method by a academic example, it was employed for determining the monolithic seepage field formed in the vicinity of a high dam and underground powerhouse, constructed in the context of a hydropower plant project in China.

2 METHODS

2.1. DESCRIPTION OF THE FREE SURFACE SEEPAGE PROBLEM

Figure 1. shows the free surface seepage flow through a dam body. The dam is characterized by the domain W=W1 U W2, however, seepage flow is actually restricted to the domain W1 Thus, the domain W1 can be referred to as the actual domain while the domain W2 can be denoted as the virtual domain. No seepage flow occurs through the free surface Γ3, which is the common boundary

of domains W1 and W2. The total waterhead at a point in the domain W1 is defined as:

... (1)

where x3 denotes the vertical coordinate, hc = p/gw the pressure water head, p the pore water pressure, and gw the unit weight of water.

The seepage velocity in the domain W1 can be calculated by Darcy's law:

... (2)

where k is the permeability tensor and V represents the gradient operator.

The seepage flow satisfies the continuity equation:

... (3)

According to Figure 1. the following boundary conditions are prescribed:

(1) The Dirichlet type boundary condition i.e. the waterhead boundary condition

... (4)

where Γi is the part of the boundary with the prescribed waterhead h.

(2) The Neumann type boundary condition i.e. the flux boundary condition

... (5)

where n and qn are the outward unit normal vector to the corresponding part of the boundary and the respective prescribed flux in the direction normal to the boundary, respectively, and Γ2 is the part of the boundary with the prescribed flux.

(3) The Cauchy type boundary conditions

... (6)

on the free surface Γ3 = EG.

(4) The composite boundary conditions on the seepage surface T4

... (7)

2.2. FE FORMULATION OF FREE SURFACE SEEPAGE PROBLEMS

2.2.1. FE discretization of free surface seepage problems

For the free surface seepage problem shown in Fig. 1, the domain W = W1 U W2 is divided into n finite elements. Within a particular finite element e, the total water head h is approximated by the vector of nodal water head values he and a vector of shape functions Ne as h = (Ne)T he (MAO et al., 1999). The corresponding pressure water head he is then computed from Eq. (1). hc is positive within W1 and negative within W2. Hence, W1 is denoted as the actual domain of seepage flow, whereas W2 is denoted as the virtual domain. Accordingly, finite elements in the domain W1 are denoted as actual finite elements, whereas finite elements in the domain W2 are denoted as virtual elements. Finite elements passed through by the free surface, which are characterized partly by h > 0 and partly by h < 0, are referred to as transition elements.

Substituting Eq. (2) into Eq. (3), gives the governing differential equation. Discretization of the weak form of the governing differential equation yields the FE balance equation in the seepage domain W1

... (8)

In (8) K1 denotes the global permeability matrix for the domain W1, obtained by assembling the element permeability matrices ke located in domain W1, h denotes the vector of the unknown nodal values of the total waterhead, and q1 represents the vector of nodal fluxes in W1 computed from prescribed values of the waterhead on Γ1 and Γ1. As the domain W1 and the seepage surface Γ4 are unknown in advance, the solution of the free surface seepage problem cannot be obtained directly from Eq.(8). Hence, an iterative solution scheme is proposed. It is characterized by improving the virtual flux method (VFM) (ZHU & SU, 1991; ZHU, 1997) by a refined Gauss point method combined with a regularized Heaviside function (CHEN et al., 2008) for determining the free surface.

2.2.2. Improved virtual flux method (IVFM)

As the domain W1 is unknown in advance, Eq. (8) is replaced by the respective FE balance equation for domain W

... (9)

where K represents the global permeability matrix for domain W and q the corresponding flux vector computed from the prescribed values of the waterhead.

The basic idea of the VFM is to determine the virtual domain W2 in an iterative manner and to gradually deduct the virtual flux contribution of the nodes in domain W2 until there is no flux across the interface between W1 and W2. The boundary conditions on the free surface Γ3 in Eq. (6) are then met automatically and the free surface can be identified by interpolating the water head according to the requirement of h = x3 within the transition elements.

Similarly, in Eq. (8), the FE balance equation for the virtual domain W2 is given as:

... (10)

where K2 is the global permeability matrix for the domain W2 and q2 denotes the corresponding nodal flux vector.

From Eq. (8), Eq. (9) and Eq. (10) it then follows that

... (11)

Substituting Eq. (11) into Eq. (8) yields the FE balance equation within the framework of the IVFM:

Kh = q - q2 + K2h (12)

Making use of Eq. (9) and Eq. (12), the iterative form of the IVFM can be written as

... (13)

with i denoting the iteration step; hi and qi are the nodal waterhead and the nodal flux vectors for the domain W at iteration step i, whereas K2 and q2 denote the permeability matrix and the flux vector for the virtual domain W2 at the respective iteration step.

In the IVFM, the solution of the free surface seepage problem is obtained by gradually modifying domain W2 and, hence, q2 and K2 in each iteration step. Based on the first estimate of the waterhead, computed from h0 in the first part of the Eq. (13), an initial estimate of the virtual domain W2 and, hence initial estimates of K2 and q2, denoted as K2 and q2 in the second part of Eq. (13), is determined. Subsequently, the second part of Eq. (13) is solved for hi, yielding an improved estimate of the waterhead and, consequently, an improved estimate of the domain W2. The iteration is continued until the convergence criterion for the nodal waterhead vector

... (14)

is met, where er denotes an error threshold value.

q2 represents the nodal flux vector in domain W2, which is contributed by the prescribed boundary values of the waterhead as well as by the prescribed boundary values of the flux. Hence, in general, non-zero items in q2 are only computed for transition elements with prescribed waterheads on Γ1 and Γ4. Since the contribution from the latter is small compared to q( q!2 could be neglected.

Since the virtual domain consists of virtual elements and fractions of transition elements, K2 is obtained by assembling the permeability matrices ke of both the virtual elements and the respective virtual fractions of the transition elements. Whereas calculating ke for the virtual elements is straightforward, it is not for the virtual fractions of the transition elements.

2.2.3 Permeability matrix ke of a transition element

The permeability matrix ke of a particular finite element e is obtained by Gauss-Legendre integration, yielding for a three-dimensional iso-parametric finite element

... (15)

where ng is the number of integration points in each direction and Wi, Wj and Wm and Wm are the respective integration weights; Fe(⅞.) = [(Be)TkeBe det Je] is the integrand, evaluated at the integration point ,nj ,Çm), with Be as the geometric matrix, containing the derivatives of the shape functions and detJe as the determinant of the Jacobian matrix.

To compute the permeability matrix ke for a transition element, a refined Gauss point method (WANG & HUANG, 1997) combined with a regularized Heaviside function (CHEN et al., 2008) is implemented in the framework of the IVFM.

(1) The refined Gauss point method

A Gauss quadrature point in a transition element is either assigned to the virtual domain if its waterhead hc is negative or to the actual domain if hc is positive. For determining the contribution of the element permeability matrix ke to the global permeability matrix K2 for a transition element, the integration of ke is then carried out only for the Gauss points within the virtual domain. The accuracy of calculating ke can be improved by increasing the number of Gauss points (SHU et al., 2007), as shown in Fig. 2. In the developed code, up to 7 Gauss points in each direction may be used for a 3D iso-parametric finite element. The increase of the number of Gauss points has an effect similar to refining the mesh in the vicinity of the free surface.

However, depending on the computed value of the waterhead hc at a particular Gauss point, the subregion of a finite element associated with this Gauss point is still entirely assigned to either the actual domain or the virtual domain. This shortcoming is shown for example in Fig. 3. If, for example, the free surface passes through the region between Gauss point I and node 1, the actual domain in the shaded part 1 would be wrongly treated as

a virtual domain due to negative values of hc at all Gauss points of this element. Similarly, the virtual domain in the shaded part 2 would be wrongly treated as an actual domain due to positive values of hc at all Gauss points of this element.

(2) The regularized Heaviside function

To overcome the aforementioned shortcoming, a regularized Heaviside function (SU et al., 1999; ZHANG & WU, 2005)

... (16)

is introduced with the parameters e1 < 0 and e2 > 0 (Fig. 4(a)). For a 2D finite element Fig.4(b) illustrates the application of the regularized Heaviside function for the two cases shown in Fig. 3 and the definition of the regularization parameters e1 and e2.

By implementing the above regularized Heaviside function, the permeability matrix ke of the virtual domain in a transition element is computed as

... (17)

where ... is the water head at the Gauss point Zi ,4j Zm ).

2.2.4. Determination of the seepage surface

As expressed by Eq. (7), on the seepage surface Γ4 the flux is non-negative and the waterhead is equal to the vertical coordinate x3 of the respective point on the boundary. According to the equivalent nodal flux method proposed (ZHU & ZHANG, 1997), the total flux on the seepage surface is equivalent to the sum of the nodal values of the flux on the respective surface. Thus, the boundary condition for the seepage surface qn > 0 is replaced by the condition qn > 0 for the node i, located on the seepage surface.

Since the seepage surface Γ4 = DE in Fig. 1 is unknown in advance, it is determined in an iterative manner. In the first iteration step, it is assumed that all nodes on the possible seepage surface DF in Fig. 1 satisfy the water head requirement h = x3. Then the fluxes at the nodes along DF are calculated. If, for a particular node qn > 0 holds, then the respective node is considered to belong to the seepage surface Γ4 in the next iteration step.

2.3. CONSIDERATION OF DENSELY SPACED DRAINAGE HOLES

Drainage holes in hydraulic or geotechnical engineering are commonly classified into leakage drainage holes and overflow drainage holes according to the respective seepage behaviour.

Figure 5. shows the drainage galleries and arrays of drainage holes arranged at the upstream side of an underground powerhouse. Drainage holes of array 1 are drilled in upward direction from the roof of drainage gallery 1, while drainage holes of array 2 and array 3 are drilled between the upper gallery and the lower gallery. In these drainage holes, shown in Fig. 6(a), the water leakage from the hole walls is discharged into the lower galleries. In this study, they are classified as leakage drainage holes satisfying the boundary conditions

... (18)

... (19)

A leakage drainage hole is inactive if it is entirely located above the free seepage surface. In this study, drainage galleries are also treated as leakage drainage holes.

In Fig. 5, the drainage holes of array 4 are drilled in a downward direction from the bottom of the lowest gallery and they are classified as overflow drainage holes since water flows freely from the top of these holes when the holes are filled with water.

In this case, shown in Fig.6(b), the waterhead for nodes located on the hole wall is equal to that of the nodes at the orifice of the hole, i.e.

... (20)

where h' represents the waterhead at node a or a'.

The power of the modern computer allows the direct modelling of densely-spaced drainage holes of small diameter. In this study, the wall of each drainage hole is considered as an internal drainage boundary and is integrated into the overall numerical model for determining the free seepage surface. In the numerical simulation procedure, all the drainage holes are treated as boun- dary conditions. Hence, the related iteration is the same as the seepage surface by the IVFM. As a result, the overall seepage field and the seepage field in the vicinity of densely-spaced drainage holes are computed simultaneously.

For reasonable representation of well drainage holes with small diameters in a numerical model, a fine discretization is required in the vicinity of such drainage holes, whereas for the overall seepage field a considerably coarser discretization is employed. Hence, special attention has to be paid to the smooth transition between domains with fine and coarse meshing. In the present context, as shown in Fig. 7, a substructure technique is implemented for generating FE meshes with densely-spaced drainage holes. The basic idea is to discretize the domain under consideration disregarding the drainage holes. Finite elements containing drainage holes serve as super-elements, which are subdivided further radially in the direction of the drainage holes into two or more layers of sub elements. The substructure technique allows eliminating nodal unknowns in the interior of the super-elements. The number of the drainage holes embedded in a super-element is determined from the size of the super-element and the spacing between two drainage holes. Figure 8. shows the planar graphs of various subdivisions in one super-element and Fig. 8 (a), (b), (c), (d), and (e) shows one hole two layers of subelements, one hole three layers, one hole four layers, two holes four layers, and three holes four layers respectively.

In this direct simulation method, the operation of each drainage hole is controlled as follows. Assume that qt and qb denote the flux at the top and the bottom nodes of a drainage hole, respectively. A leakage drainage hole is inactive, if qb is negative, i.e. the drainage hole is located in the virtual domain above the free surface. Otherwise, the leakage drainage hole is active. An overflow drainage hole is only active if it is located entirely in the actual domain below the free surface, characterized by a positive value of qt.

The surface of an active leakage drainage hole satisfies the boundary conditions prescribed by Eq. (18) and Eq. (19), and the seepage surface is determined using the method presented in section 3.3. Inactive leakage drainage holes are ignored. Generally, the boundary conditions of drainage galleries or tunnels are similar to those of leakage drainage holes. Thus, drainage galleries or tunnels are treated like leakage drainage holes.

The surface of an active overflow drainage hole satisfies the boundary condition prescribed by Eq. (20). The waterhead at the nodes of the surface is usually equal to the vertical coordinate of the top node of the hole. If the overflow drainage hole is inactive, then the whole boundary is assumed to be impermeable.

Whether a drainage hole is active or inactive is determined in the course of the iterative solution procedure for the seepage surface, described in section 3.2. At the beginning of the first iteration step all drainage holes are assumed to be in the active state. Based on the results of the current iteration step, the fluxes qb for leakage drainage holes and the fluxes qt for overflow drainage holes are computed, then the states of the drainage holes are updated accordingly. Subsequently, the boundary conditions of the drainage holes are adapted appropriately for the next iteration step. The iteration will be terminated, if the states of all drainage holes do not change in two adjacent iteration steps and if the convergence criterion (14) for the nodal waterhead vector is met.

3. RESULTS

3.1.Verification

The analytical solution for the free seepage surface in a homogeneous dam of rectangular shape is selected for verifying the IVFM. As shown in Fig. 9, the width and height of the dam were chosen as 10.0m and 12.0m, respectively. The prescribed waterheads at the upstream and downstream face are 10.0m and 2.0m, respectively. The foundation of the dam is considered as an impermeable boundary. The analytical solution of the free seepage surface is given (ZHOU et al., 1996) as:

... (21)

A regular FE-mesh, consisting of quadrilateral finite elements with dimensions of 0.5mx0.5m, is used for discretizing the dam body.

Figure 10. shows the computed location of the seepage point at the downstream face, employing an increasing number of Gauss points for the numerical integration for consistently refined FE-meshes. It confirms that the numerical solution gradually approaches the analytical solution with an increasing number of Gauss points.

Figure 11. contains a comparison of the analytical solution for the free surface, obtained from Eq. (21), with numerical solutions computed on the basis of the virtual flux method (VFM) (ZHU, 1991) and the present IVFM. The numerical solution by means of the VFM agrees well with the analytical solution, the IVFM allows even better reproduction of the analytical solution. Hence, the proposed IVFM can be employed for analyzing free surface seepage problems with a greater accuracy.

In a further step, a drainage gallery, characterized by a cross section of 1.0mx1.0m, a leakage drainage hole of 6.5 m height and a width of 0.05m and an overflow drainage hole with the size of 1.5m length and 0.05 m width are considered in the above homogenous rectangular dam. The computed free seepage surface and the waterhead contours are shown in Fig. 12. It can be seen that the free surface depresses rapidly through the leakage drainage hole and towards the downstream water level when it passes through the drainage gallery and the overflow drainage hole is filled with water.

3.2.Large-scale application in engineering practice

In this section, the IVFM, extended by considering densely spaced drainage holes as internal boundaries according to section 4, is applied for determining the monolithic seepage field in the vicinity of a hydropower station currently under construction in China.

3.2.1.Hydrogeological features

The water conservancy project is located in an L-shaped bay area where the river is deeply-incised in a V-shaped valley. The river is always in a turbulent state and the altitude, depth and width of the river are about 676~678 m, 7~11 m and 60~80 m respectively. The right bank is the concave bank of the river while the left bank is the convex bank.

The thick river bed cover of the dam foundation consisted of a floating pebble layer, pebble layer, sand - lens layer with a drifting pebble layer and an erratic boulder layer. In general, the thickness is about 40.0~60.0 m and can reach to 75.36 m in the deeplyincised channel area. i) the floating pebble layer is mainly distributed in a grade II terrace of the left bank upstream of the dam axis with a thickness of 40~50 m. The front borehole of the terrace shows that the maximum thickness is about 70.72 m and the altitude of the top surface of the layer is about 730.0~734.0 m. The basal layer directly overlies the bedrock at an altitude of about 620.0~660.0 m. ii) this pebble layer is distributed in both a I grade terrace and the bottom of riverbed. The residual thickness at the base of the riverbed is about 22.0 ~32.0 m. The altitudes of the upper and lower surfaces are about 635.0~640.0 m and 597.0~620.0 m, respectively. iii) the distribution of the drifting pebble layer is mainly on the surface of the current riverbed and the floodplains, overlapping the sand-lens layer with a drifting pebble layer with a thickness of 10.0~25.0 m.

The project area is mainly composed of pre-Sinian epimetamorphic basalt, tuff and rhyolite porphyry of the Suxiong Formation of the Lower Sinian, with the Chengjiang granite and a Quaternary loose accumulation layer. i) The bed rock of PreSinian hypometamorphic basalts is mainly exposed in the right bank valley slope and buried in the bottom of the river bed. ii) The lower Sinian Suxiong Formation rocks exhibit well developed columnar joints with columnar diameters ranging from 0.3 m to 0.8 m and the distribution area is large, mainly exposed on both sides downstream of the junction area. iii) the Chengjiang period granite intrusions are widely exposed in the left bank upstream and downstream of the junction area. The right bank represents only a small amount of the distribution in the Niri estuary. The rock mass is hard, complete and of good quality which is the main engineering utilization rock mass within the scope of the project.

3.2.2. Dam structure and power station

As shown in Fig. 13, the underground powerhouse is located on the left bank of the dam. At the upstream side of the underground powerhouse, three drainage galleries and four arrays of drainage holes, the latter with a diameter of 76.0mm each are equally spaced at 4.0m, are provided, as shown in Fig. 5. Figure 14. shows the cross section of the rock-filled dam with a maximum height of 186m. It is built on a deep overburden layer with a maximum thickness of 75m. The anti-seepage system of the dam and the foundation includes the earth core, two concrete cut-off walls and the grouting curtains. The normal water level of the reservoir and the tailwater level are 850.0m and 670.0m, respectively.

3.2.3. Numerical simulation

Figure 15. shows the three-dimensional FE mesh for this hydropower station, consisting of 74592 nodes and 67748 linear 3D finite elements. It consists of the dam body, the underground powerhouse, (approximated by a cuboid), the surrounding rock mass and 147 leakage drainage holes and 49 overflow drainage holes. Because of the large number of drainage holes, the method for considering densely spaced drainage holes in section 4 serves as an essential means for the efficient modeling of drainage holes at a similar accuracy by means of the substructure technique, however, at reduced computing time. Because of the underground powerhouse, a larger domain of the rock mass is discretized on the left bank of the dam. The permeability coefficients of the different materials are summarized in Table 1. Setting er = 10 in the convergence criterion of Eq. (14), requires 16 iteration steps to compute the free seepage surface. The computing time on a PC AMD Athlon (TM) 64 X2 Dual Core Processor 4000+, 2.11GHZ, memory 2.00G RAM, is about 3200 seconds. Compared with the direct method employed, considering drainage holes as internal boundaries, means that the substructure method is somewhat slower (Table 2).

Figure 16. shows the plan view of the distribution of the free surface in the monolithic seepage field. It nearly matches the topography and is characterized by steep gradients in the vicinity of the dam body and the underground power house, reflecting the complex boundary conditions and the effects of the seepage-control and drainage measures. The waterhead isolines, displayed in Figure 16. show a significant drop of the water head in the vicinity of the underground powerhouse. The water level of 650m at the upstream side of the powerhouse is almost equal to the bottom elevation of the lowest drainage gallery. Thus, the drainage galleries and densely-spaced drainage holes at the upstream face of the powerhouse significantly reduce the water pressure acting on the powerhouse.

The waterhead isolines for the vertical typical sections of dam and powerhouse are displayed in Figs. 17 and 18, respectively. Fig. 17 shows that: (i) the upstream water level in the rock mass is higher than the water level of the reservoir, (ii) the downstream water level in the rock mass is higher than the tailwater level, (iii) the free seepage surface drops from 850m to 675m inside the earth core, accounting for 97% of the water level differ- ence between the upstream and the downstream side of the dam, (iv) a large reduction of the waterhead is caused by the concrete walls and the grouting curtains, effectively cutting off the seepage passing through the overburden foundation, whereas smooth changes of the waterhead occur in the deep overburden.

The waterhead isolines, displayed in Fig. 18, show a significant drop of the free surface in the vicinity of the underground powerhouse. The water level of 650m at the upstream side of the powerhouse is almost equal to the bottom elevation of the lowest drainage gallery. Thus, the drainage galleries and densely-spaced drainage holes at the upstream face of the powerhouse significantly reduce the water pressure acting on the powerhouse.

4. CONCLUSIONS

Monolithic seepage problems with complex drainage systems are commonly faced in hydraulic engineering, slope engineering and underground engineering. Typically, they are non-linear problems. In this study, a new FE method was presented for monolithic seepage problems, which is characterized by analyzing the overall unconfined seepage flow by the proposed improved virtual flux method (IVFM) simultaneously with the local seepage field in the vicinity of the complex drainage system. The latter usually consists of leakage drainage holes and overflow drainage holes, which are modeled as internal boundaries.

Compared to the previously proposed virtual flux method, in the IVFM the accuracy of the computed conductivity matrices in transition elements, i.e., the finite elements which contain the free surface, is improved by increasing the number of Gauss points for the numerical integration and by introducing a regularized Heaviside function for distinguishing between the domain below and above the free surface of a particular transition element in the integration procedure. These improvements allow the use of coarser FE-meshes and in addition, adaptive remeshing can be avoided. These features are especially important for the analysis of large-scale 3D seepage problems in engineering practice.

The proposed method was verified by comparing numerical and analytical results for an academic seepage problem, for which the analytical solution is available. Subsequently, the effectiveness and robustness of the proposed method was demonstrated by the large-scale 3D numerical simulation of the monolithic seepage field in the vicinity of a hydropower station with a complex drainage system, currently under construction in China.

ACKNOWLEDGEMENT

This research was financially supported by the National Key Technologies Research & Development program (2017YFC0804600) and National Natural Science Foundation of China (Grant No. 51379220).

References

REFERENCES

ALT, H.W. (1980): Numerical solution of steady state porous flow free boundary problems.Numerische Mathematik, 36/1, 73-96. doi:10.1007/BF01395990.

AZUSA, O., TAKAHIRO, H., DENNIS, B. & JUN, S. (2018): Evaluations ofthe downward velocity of soil water movement in the unsaturated zone in a groundwater recharge area using 518O tracer: the Kumamoto region, southern Japan.- Geologia Croatica, 71/2, 65-81. doi: 10.4154/gc.2018.09.

BATHE, K.J, KHOSHGOFTA, A.R & MOHAMMAD, R. (1979): Finite element free surface seepage analysis without mesh iteration.- International Journal for Numerical & Analytical Methods in Geomechanics, 3/1, 13-22. doi:10.1002/ nag.1610030103.

BORJA, R.I. & KISHNANI, S.S. (1991): On the solution of elliptic free boundary problems via Newton's method.- Computer Methods in Applied Mechanics & Engineering, 88/3, 341-361. doi:10.1016/0045-7825(91)90094-M.

CHEN, S.H., XU, Q. & HU, J. (2004): Composite element method for seepage analysis of geo-technical structures with drainage hole array.- Journal of Hydrodynamics, 16/3, 260-266. doi:10.1016/jjfluidstructs.2003.09.001.

CHEN, S.H., XUE, L.L. & XU, G.S. (2010): Composite element method for the seepage analysis of rock masses containing fractures and drainage holes.- International Journal of Rock Mechanics and Mining Sciences, 47/5, 762-770. doi:10.1016/j. ijrmms.2010.03.011.

CHEN, Y.F., ZHOU, C.B., & ZHENG, H. (2008): A numerical solution to seepage problems with complex drainage systems.- Computers and Geotechnics, 35/3, 383-393. doi: 10.1016/j.compgeo.2007.08.005.

DESAI, C.S. & LI, G.C. (1983): A residual flow procedure and application for free surface in porous media.- Advances in Water Resources, 6/1, 27-35. doi: 10.1016/0309-1708(83) 90076-3.

FIPPS, G., SKAGGS, R.W. & NIEBER, J.L. (1986): Drains as a boundary condition in finite elements.- Water Resources, 22/11, 1613-1621. doi: 10.1029/ WR022i011p01613.

GABRIELLA, B., KISS, F.M. & LADISLAV, A.P. (2016): Hydrothermal processes related to some Triassic and Jurassic submarine basaltic complexes in northeastern Hungary, the Dinarides and Hellenides.- Geologia Croatica, 69/1, 39-64. doi:10.4154/gc.2016.04.

JANSEN, R.B., KRAMER, R.W., III, J.L. & POULOS, S.J. (1988): Earthfill Dam Design and Analysis.- Advanced Dam Engineering for Design, Construction, and Rehabilitation, 32/1, 305-308. doi: 10.1007/978-1-4613-0857-7_9.

JING, H.U., QU, W.L. & JIANG, Y.M. (2005): Study on parameter of composite element method for modeling drain hole in seepage analysis.- Water Resources & Hydropower Engineering, 6/1, 56-57.

MAO, C.X., DUAN, X.B. & LI, Z.Y. (1999): Numerical computation in seepage flow and programs application. - Hohai University Press, Nanjing, 3M-0.

NEUMAN, S.P. (1973): Saturated-unsaturated seepage by finite elements.- Journal of the Hydraulics Division Asce, 99/2, 2233-2250.

SHU, Z.Y., DENG, J.X. & LI, L.G. (2007): Application of Element Conduct matrix Regulative Method by Encrypting Gauss Point in Seepage Analysis with Free Surface.- Journal of Sichuan University, 6/1, 48-52. doi:10.3969/j. issn.1009-3087.2007.01.010.

SU, B.Y., GUO, H.X. & ZHAN, M.L. (1999): Research on the calculation of unsteady seepage by cut-off negative pressure method.- Chinese Jounal of Geotechnical Engineering, 5/2, 711-714. doi:10.1016/jjfluidstructs.2004.09.004.

WANG, X.N. & HUANG R.Q. (1997): Gauss point method for seepage analysis with free surface.- Hydrogeology and Engineering Geology, 6/3, 1M-. doi:10.16030/j. cnki.issn.1000-3665. 1997.06.001.

WANG, Y. (1998): The modified initial flow method for 3-D unconfined seepage computation.Journal of Hydraulic Engineering, 29/3, 68-73. doi: 10.3321/j. issn:0559-9350.1998.03.015

ZHANG, Y.T., CHEN, P. & WANG, L. (1988): Using initial flow method to analyze the unconfined seepage surface.- Journal of Hydraulic Engineering, 29/2, 18-26.

ZHANG, Q.F. & WU, Z.R (2005): The improved cut-off negative pressure method for unsteady seepage flow with free surface.- Journal of Geotechnical Engineering, 27/1, 48-53. doi:10.3321/j.issn:1000-4548.2005.01.007.

ZHENG, H., DAI, H.C. & LIU, D.F. (2005): Improved Bathe's algorithm for seepage problems with free surfaces.- Rock and Soil Mechanics, 26/4, 505-512. doi:10.3969/j.issn.1000-7598.2005.04.001.

ZHENG, H., LIU, D.F., LEE, C.F. & THAM, L.G. (2005): A new formulation of Signorini's type for seepage problems with free surfaces.- International Journal for Numerical Methods in Engineering, 64/1, 1-16. doi: 10.1002/nme.1345.

ZHOU, C.B., XIONG, W.L. & LIANG, Y.G. (1996): A new method for the solution unconfined seepage field.- Journal of Hydrodynamics, 10/1, 528-534.

ZHU, Y.M. (1997): The equivalent nodal flux method for the calculation of Darcy flow.Journal of Hohai University, 7/1, 105-108.

ZHU, Y.M. & SU, B.Y (1991): Procedure of virtual flux for seepage free surface with fixed mesh.- Journal of Hohai University, 5/2, 113-117.

ZHU, Y.M. & ZHANG, L.J. (1997): Solution to seepage field problem with the technique of improved drainage substructure.- Journal of Geotechnical Engineering, 19/1, 69-76.