Abstract
Water flow and hydromechanical coupling process in fractured rocks is more different from that in general porous media because of heterogeneous spatial fractures and possible fracture-dominated flow; a saturated-unsaturated hydromechanical coupling model using a discontinuous deformation analysis (DDA) similar to FEM and DEM was employed to analyze water movement in saturated-unsaturated deformed rocks, in which the Van-Genuchten model differently treated the rock and fractures permeable properties to describe the constitutive relationships. The calibrating results for the dam foundation indicated the validation and feasibility of the proposed model and are also in good agreement with the calculations based on DEM still demonstrating its superiority. And then, the rainfall infiltration in a reservoir rock slope was detailedly investigated to describe the water pressure on the fault surface and inside the rocks, displacement, and stress distribution under hydromechanical coupling conditions and uncoupling conditions. It was observed that greater rainfall intensity and longer rainfall time resulted in lower stability of the rock slope, and larger difference was very obvious between the hydromechanical coupling condition and uncoupling condition, demonstrating that rainfall intensity, rainfall time, and hydromechanical coupling effect had great influence on the saturated-unsaturated water flow behavior and mechanical response of the fractured rock slopes.
1. Introduction
Fractured rock masses including numerous discontinuities with various attitudes and different scales are complicated geological media that have undergone a long period of geological evolution. And the corresponding prediction and description of deformation, stress, and groundwater flow in the fractured rocks are of interest in many geotechnical areas such as slope stability [1, 2] and seepage control for hydraulic projects [3, 4]. However, different from the general porous media, the permeability of intact rock is very low compared with that in fractures [5]; water flows predominantly through a single fracture or fracture networks [6]; it is more complex to describe the seepage behavior in such rocks. From the achievements about water flow in fractured rocks, numerical methods are often used to model seepage behaviors through fractured rocks, including the equivalent continuum approaches such as FDM, FEM [7], and BEM, the dual-continuum approach introduced by Barenblatt (1960) commonly used in mimicking water flow and mass transport in fractured porous medium [8, 9], and discrete fracture network (DFN) [10–13]. However, the continuum model is to view the fractured rock masses as an equivalent porous medium in which the heterogeneity of rock masses is averaged based on the concept of representative element volume (REV). Thus, it is desirable to evaluate spatial and random distribution of the fractured rock flow [12], and then a sound numerical method of intact rock and fractures should be proposed to model water movements in the deformed rocks for performance assessments and design optimization in a lot of rock engineering such as rock slopes and dam foundation.
And above researches are mainly about rock steady and nonsteady seepage focusing on water flow in the fully saturated mode; the unsaturated flow behavior in such fractured rock masses is very often encountered in geotechnical problems such as atomized rain or reservoir water level variation on the dam and foundation, embankment, and rainfall infiltrating into the rock slope. Whereas few investigations about saturated-unsaturated hydromechanical coupling model in the fractured rocks have been implemented, so the unsaturated flow behavior especially reflecting the hydromechanical coupling effect in fractured rocks cannot be ignored and should be deeply investigated considering the rock deformation. A discontinuous deformation analysis (DDA) [14] may be employed for time discretization and spatial approximation to solve the rock deformation. And in the DDA method, the displacements and strains of the blocks are determined by solving a system of simultaneous equilibrium equations which are formulated by minimizing the total potential energy of the block system, similar to the conventional FEM. In addition, the DDA method allows the blocks to have various modes of motion, such as sliding, rotation, and opening, so that large displacements can be simulated, which is similar to the DEM [15, 16]. Meanwhile, with continuous contributions to its improvements and applications, the DDA was well developed in terms of both theory and computing codes [17–21] to analyze the mechanical response of the rocks with random fractures or faults. And then, in order to determine the effective saturation and relative permeability [22, 23], the Van-Genuchten model was written in saturated-unsaturated seepage codes differently treating the rock and fractures permeable properties describing the constitutive relationships during the unsaturated flow processes. Thereby, combined with the DDA and seepage codes of fractured rocks, a saturated-unsaturated hydromechanical coupled model using DDA has been proposed.
For the sake of better applying the proposed coupled model to the corresponding practical engineering, the validity of the developed method was verified by the consistency of numerical results and computing results by the other method DEM in the analysis of a dam foundation. And also, a fractured rock slope considering rainfall infiltration as a case study and the water pressure distribution on the fault surface and inside the rocks was demonstrated with rainfall time increasing, and also the displacement, principle stress distribution, and slope stability were described considering the hydromechanical coupling and uncoupling process, indicating that rainfall intensity, rainfall time, and hydromechanical coupling effect under unsaturated conditions have great influence on the movement of seepage flow through fractured rock slopes and corresponding stability.
2. Unsaturated Hydromechanical Coupling Model of Fractured Rock Mass Using DDA
In view of the complexity of fractured network and hydraulic parameters determination, some traditional analysis models are limited for saturated-unsaturated fluid flow of fractured rocks. For better describing the hydraulic and mechanical characteristic of the rock masses, an unsaturated hydromechanical coupling model based on DDA was brought forward to simulate the water pressure, fracture deformation, and stress distribution to analyze the seepage field and stress field variation.
2.1. Unsaturated Flow Model of Fractured Rock Masses
It is known that the entire area of the fractured rocks can be divided into saturated domains and unsaturated domains. The movement of partially saturated flow in an unsaturated domain appears because of the capillary effect. The governing equation of unsaturated water flow in porous media can be written as where is water density; S is water resource; is porosity; is saturation, ; is water content, ; is Darcy velocity.
According to Darcy’s law, where is permeability tensor related to water content ; is head; is the saturated permeability tensor of fractured rocks .
Substituting (2) into (1) gives the unsaturated flow expressions in the fractured rocks as follows: where is the relative permeability, in unsaturated field, and in saturated domains; is pressure head and ; is water capacity and in saturated domains; is special parameter, in unsaturated domains, and in saturated filed; is the storage coefficient, for the unsaturated rock, and is constant for saturated rocks.
To solve (3) in the saturated-unsaturated flow analysis, initial and boundary conditions are needed as follows.
Initial condition is expressed as
The boundary conditions generally including the water head, flux, and seepage-face boundaries are described by where and are normal flow rate; is direction cosine of external normal; is initial time; is head boundary; is flow rate boundary; is saturated overflow boundary; is unsaturated overflow boundary.
2.2. The Theory of DDA Method
Discontinuous deformation analysis (DDA) is a numerical method that can be used to simulate the discontinuous deformation behavior of jointed rocks [14]. The method not only can model all kinds of blocks cut by the fractures, but also deals with the contact [24] and restraint between all the blocks. And the method determined by minimizing the total potential energy [19] of all blocks is similar to FEM and allowing all the blocks to have various mechanical motions is similar to DEM.
2.2.1. Displacement Functions of Deformable Blocks
Different from FEM selecting the node displacements as the degrees of freedom, DDA selects the rigid body displacements and the strains of block elements as the basic unknown variables to solve for solutions. DDA conducts calculations based upon time steps in which the large displacement and deformation are the accumulating results of the small ones obtained from time-stepwise calculations. Suppose that an arbitrarily shaped polygonal block has uniform stress and strain, the unknown variables including parallel movement , rotation , normal strain, and tangential strain of the block can be written as where are, respectively, shape matrix and displacement matrix.
2.2.2. Motion Equation of Block System
Suppose blocks in the block system, the potential energy of a block, includes the elastic strain energy , initial stress potential energy , point force potential energy , linear force potential energy , body force potential energy , inertial energy , and displacement constraint energy can be written as
According to the principle of minimum potential energy, derive (7) to obtain the stiffness matrix expression and load matrix of a single block,where and
Thus, the equations of motion of the blocks can be given in the matrix form,
2.2.3. Treatment of Water Load
Water flow behavior in the fractured rocks is significant for building the unsaturated hydromechanical coupling model; water pressure should be correctly determined on the blocks. Shown in Figure 1 representing the water pressure on a block, suppose and at points , respectively, and the water pressure at any point can be described by (10), where
Thus, the energy of th block considering the water pressure as the linear force can be expressed by
Then, derive (11) to get the right load item ,
Assuming the constant constraint element, the load item can be written as
To improve the computing precision, the higher-order nonconstant constraint element is introduced to determine the load items with the following two ways.
Assuming the regular rectangular element of a block with four points under global coordinate; two transformation matrices are written as (14) and (15). And Substituting (15) into (12) gives the corresponding load item as (16)
In general, not all the block elements are regular so that complex nonregular element is needed to be deeply developed. The displacement function will be expressed by isoparametric element function; the load item can be written by where is weight function. Suppose global coordinate of and local coordinate of and the relationship of in can be described using the shape function , written as where , , are the vertex coordinates of a block element and can be expressed by (20). And substituting (21) into (19) gives the expressions of weight function , written as (22) and shape function as (23); the load item can be expressed by (24). The expression of the water load on the block system using DDA can realize the hydromechanical coupling effect coming from the water flow inside the block and the block deformation where and
2.3. Unsaturated Hydromechanical Coupling Model of Fractured Rock Masses Using DDA
The unsaturated coupled model using DDA has been proposed to solve the water pressure variation and mechanical response in saturated-unsaturated rocks. Note that the aperture variation of rock fractures resulting in obvious water pressure variation and rock deformation, corresponding real-time coupling process cannot be ignored.
However, the real fractures are coarse and nonparallel, Cubic Law cannot be directly applied in this study. For the sake of getting new permeability tensor at any time, fractures deformation including normal deformation and shear deformation [24] should be completely dealt based on DDA calculation. As mentioned in previous studies about water flow in the fractured rocks, natural fractures cannot completely close under high normal stress with enough mechanical aperture. Assuming , are the initial fractures aperture for each side of the fracture, , representing the fracture aperture after rock deformed are written by (25) considering fracture normal deformations and (open is positive) and shear deformations and computed by DDA. And then, the corresponding flow rate can be expressed by (26); thus the equivalent hydraulic aperture can be calculated by (27) to get new permeability tensor at any time where is the dilation angle; ; is the flow rate; is the flow rate corresponding to average aperture ; and .
Therefore, the proposed coupled model can consider the block normal deformation and shear dilation when water is flowing through the rock fractures and can also be applied to determine the unsaturated hydromechanical coupling behavior in the fractured rocks. Shown in Figure 2 is the flow chart of the proposed coupled model using DDA, recording the whole calculating process.
3. Calibration of the Proposed Coupled Model Based on DDA
For the sake of calibrating the proposed hydromechanical coupling model, shown in Figure 3 is a model of concrete gravity dam with rock foundation that measures 50 m in height, 6 m in width at dam top, and 30 m in width at dam bottom, with the upstream face vertical and the downstream one having a slope of 1 : 0.682. The rock foundation is 105 m in length and 50 m in height. The calculation domain takes the dam heel as the original of the coordinates with -axis directing toward the downstream and -axis pointing upwards. The mechanical parameters of dam are the density 2400 kg/m3, Young’s module 19.6 GPa, and Poisson’s ratio 0.167. And corresponding mechanical parameters of rock foundation are the density 2650 kg/m3, Young’s module 8.0 GPa, Poisson’s ratio 0.24, internal friction angle 45°, and cohesion 2 MPa. Suppose a group of orthogonal fractures in the rock foundation with aperture of 1 mm, the normal stiffness, and tangential stiffness of 20 GPa/m. The boundary conditions include the horizontal displacement restraint for the upstream and downstream boundaries of the rock foundation and vertical displacement restraint for the bottom boundary of the rock foundation, and water level acting on the upstream face of the dam body is 50 m and on the downstream face is 0 m.
And then, both of the numerical models based on DDA and DEM methods were employed to make hydromechanical coupling calculations considering the horizontal water pressure, uplift pressure, and water flow force in rock foundation. For the sake of comparison, Figures 3, 4, and 5 give the calculating results of DDA and DEM models, respectively. The water head (unit: m), displacement (unit: m) and maximum principal stress (unit: Pa) based on proposed DDA model are shown in Figures 4(a), 5(a), and 6(a), indicating that the fractures in the upstream are tending to be tensile causing the fractures permeability increasing; however, the fractures in the downstream are tending to be compressed causing the fractures permeability decreasing. Therefore, the curves of water head in upstream move toward the upstream, and the curves of water head in downstream move toward the downstream; the variations are consistent with real seepage characteristics in the rocks. And it can be seen from the computing results shown in Figures 4(b), 5(b), and 6(b) based on DEM that the variations are nearly identical as determined by DDA. The above comparison of two methods demonstrated the feasibility and reliability of the proposed coupled model using DDA and can also be applied to analyze the mechanical characteristics of rock masses. In addition, the proposed model considering the heterogeneous deformation of the fracture deformation with different permeability and the better contact criterion and convergence criterion [19] will get more accurate computing results close to the real performance, so the model is superior to DEM to solve the practical hydromechanical problems in geotechnical engineering.
(a) Proposed hydromechanical model by DDA
(b) Hydromechanical model by DEM
(a) Proposed hydromechanical model by DDA
(b) Hydromechanical model by DEM
(a) Proposed hydromechanical model by DDA
(b) Hydromechanical model by DEM
4. Analysis of a Fractured Rock Slope Considering Rainfall Infiltration
As noted in previous investigations of unsaturated fluid flow in slope engineering [1, 2, 25], rainfall infiltration will increase the water pressure inside the rocks and fractures and weaken the rock strength parameters; the landslide often occurs in the rain season. And the analysis of the fractured rock slope has been employed to demonstrate the capacity and robustness of the proposed coupled model for describing the rainfall infiltrating process influencing the slope water pressure, displacement, and stress distribution. Different from the saturated steady seepage or unsaturated in porous media, rainfall infiltrating into the rocks will flow through the fractures, directly loading on the fractures surface causing the real-time variation of the hydraulic aperture; the permeability of the rock and fractures will change accordingly, which also implies the unsaturated hydromechanical coupling process. As shown in Figure 7, a rock slope with a height of 180.0 m, length of 20.0 m, and width of 1.0 m is considered. There is a set of obvious faults with initial hydraulic aperture of 1 mm developed in the rock slope. A series of computing results about water pressure inside the rock and on the fault surface, displacement, stress, and stability coefficient under different loading conditions have been analyzed in detail based on the proposed DDA method.
4.1. Case Introduction
As an initial condition, the rock slope is treated as a water head boundary assuming an external water level 110.0 m on the upstream face and 32.0 m on the downstream face, and the bottom boundaries were impermeable. The rainfall infiltration boundaries are imposed on the slope surface. And the initial saturated permeability coefficient of the rocks is assumed to be 1.39 × 10−5 m/s. For the unsaturated coupled model, water retention curves are very important to express the unsaturated rock characteristics; particularly the Van-Genuchten model (V-G) is often used to describe the constitutive relations such as - and -. Table 1 gives the constitutive relations based on V-G model [26]. And Table 2 gives the mechanical parameters of the rock and the fault.
4.2. Analysis for the Rock Slope Considering Rainfall Infiltration
Suppose that the rainfall lasts for 80 h with rainfall intensity of 8 mm/h and 20 mm/h, and analyze the water pressure (unit: m) at different time. The distribution of water pressure at different time is shown in Figures 8 and 9, respectively. Negative water pressure indicates unsaturated area; positive water pressure demonstrates fully saturated area, which is in accordance with the rainfall filtration mechanism of the rock slope. It can be observed from these figures that, due to the rainfall infiltration, water content of the rock slope is increasing and the matric suction is reducing gradually, and the negative pressure region continues to reduce progressively; the expansion of the transient area leads to large increase of water pressure inside the rock slope, which can drive the wetting front quickly into a deeper location. In particular, with rainfall infiltrating inside the slope, there will appear greater positive transient water pressure and larger range of transient saturated zone on the rock slope surface. The results also show that transient saturated zone will quickly be larger in shorter period of time with greater rainfall intensity.
(a) Rain falling for 20 h
(b) Rain falling for 40 h
(c) Rain falling for 60 h
(d) Rain falling for 80 h
(a) Rain falling for 20 h
(b) Rain falling for 40 h
(c) Rain falling for 60 h
(d) Rain falling for 80 h
In addition, Figure 10 presents the water pressure (unit: m) distribution on the fault surface at different rainfall time, indicating the flow behaviors of rainfall infiltration inside the fault. It can be observed from Figure 10(a) with rainfall intensity of 8 mm/h that water pressure on the fault surface varies slightly, whereas it varies more greatly with rainfall intensity of 20 mm/h shown in Figure 10(b) considering the same rainfall infiltrating time. Meanwhile, as for the larger range of the unsaturated zone inside the rock slope, greater rainfall intensity will cause more and more water infiltration into the unsaturated zone to become the saturated zone, and the groundwater level will ascend quickly so as to increase the water pressure on the fault surface. Therefore, the water pressure on the sliding surface in the unsaturated zone will increase when rain falling time lasts with same rainfall intensity. And it can be seen from the two figures that, in the same period of rainfall infiltrating time, greater rainfall intensity will cause larger water content of the rocks in a shorter period resulting in greater unit weight increasing, and the water pressure head will rise more obviously so that the water pressure in the transient unsaturated zone increases greater and more quickly, which demonstrates that the rock slope will be sliding more easily considering rainfall infiltration of greater rainfall intensity.
(a) Rainfall intensity of 8 mm/h
(b) Rainfall intensity of 20 mm/h
From the above calculating results in agreement with the water flow behaviors in the saturated-unsaturated rocks, it can be noticed that the proposed coupled model is also applicable to analyze the unsaturated hydromechanical coupling characteristics under complex conditions.
4.3. Comparison under Unsaturated Hydromechanical Coupling and Uncoupling Conditions
Based on the deeper calibration of the proposed coupled model in Section 4.2, further researches will be advanced to clearly compare the mechanical characteristic of the rock slope under unsaturated hydromechanical coupling conditions and uncoupling conditions. Figures 11 and 12 show the water pressure head distribution under coupling conditions and uncoupling conditions at time T = 60 h and T = 80 h with the rainfall intensity of 8 mm/h and 20 mm/h. It can be seen from the above figures that the water pressure on the fault surface will be larger in the transient saturated zone considering the coupling process. And the comparison under two conditions demonstrates that coupling process causes water pressure head increase more obviously with longer rainfall infiltrating time and greater rainfall intensity, which is more unfavorable for the slope stability.
(a) Rain falling for 60 h
(b) Rain falling for 80 h
(a) Rain falling for 60 h
(b) Rain falling for 80 h
Furthermore, shown in Figures 13 and 14 are the maximum principle stress contours (unit: MPa) and displacement contours (unit: m) under coupling conditions and uncoupling condition considering the rainfall intensity of 20 mm/h at time = 20 h, = 60 h, and = 80 h. It can be noted from the comparison of the figures that the stress and displacement of the rock slope are larger under coupling condition than that under uncoupling conditions, which also indicates that the unsaturated hydromechanical effect will enlarge failure probability of the slope and should not be ignored especially under rainfall infiltration conditions.
(a) Rain falling for 20 h
(b) Rain falling for 60 h
(c) Rain falling for 80 h
(a) Rain falling for 20 h
(b) Rain falling for 60 h
(c) Rain falling for 80 h
For the sake of explaining the coupling effect influencing on the rock slope stability, Table 3 gives the stability coefficient of the rock slope under different conditions. The results indicate that the safety coefficients under hydromechanical coupling conditions are smaller than that under uncoupling conditions, which is in agreement with the water pressure, stress, and displacement variation. It can be concluded that the rainfall intensity, rainfall time, and hydromechanical coupling effect have significant influence on the slope stability.
Therefore, the hydromechanical coupling process under rainfall infiltration for a rock slope will cause the transient saturated zone to be larger, and the transient water pressure will accordingly increase the downslide force and water pressure causing the effective normal stress on the sliding surface reducing too much, which is more unfavorable for the rock slope stability. The comparison between uncoupling conditions and coupling conditions proves the applicability of the unsaturated hydromechanical coupled model based DDA to analyze the mechanical response of related practical projects.
5. Conclusions
DDA similar to FEM and DEM has attracted increasingly wide attention in the geotechnical engineering; it can be applied to analyze the rock deformation and fractures aperture variation. And for the hydromechanical coupling issues of the fractured rocks, the fractures aperture variation changing the permeability will seriously influence the water pressure distribution; DDA method is more favorable to solve the fractures deformation under complex conditions. Therefore, an unsaturated coupled model based on DDA has been proposed for the analysis of saturated-unsaturated water flow behavior and mechanical response of the fractured rocks such as dam foundation and rainfall infiltration in rock slopes, and a systematic numerical code was developed to analyze the hydromechanical process, indicating the feasibility and applicability of the proposed coupled model. On the basis of this work, the following conclusions can be drawn:
For better modelling the deformation of the rocks and fractures, the higher-order isoparametric element is selected in DDA procedure to improve the numerical precision. And on the basis of Cubic Law, the fracture aperture under complex conditions can be derived considering the fracture norm deformation and shear deformation to calculate the permeability at any time. And also selecting the V-G model for describing the constitutive relations and appropriate evaluation of the capacity term, the unsaturated coupled model has been proposed to describe the hydromechanical response of the fractured rocks. Meanwhile, the effectiveness of the model applying in a dam foundation has been calibrated by the satisfactory consistency of the results by DDA and computing results by DEM.
The proposed model has been applied to describe the water variation at various time under a rainfall infiltration process in a rock slope and the unsaturated flow behavior with respect to the variation of water pressure inside the rocks and on the fault surface is well revealed. It can also be noted that the matric suction of the rock slope is smaller and smaller when rainfall infiltrating lasts, and greater rainfall intensity causes safety coefficient to reduce quickly, which is in good agreement indicating the feasibility of the proposed model.
Further research has been advanced by comparison of the pressure head, principle stress, and displacement under coupling conditions and uncoupling conditions, demonstrating that the above corresponding values are increasing and the safety coefficients are reducing at the same time under coupling conditions, and also indicating that the rainfall intensity, rainfall time, and hydromechanical coupling effect under unsaturated conditions have great influence on the stability of the rock slope. It also explains that it is more dangerous under rainfall infiltration in a rock slope considering the hydromechanical coupling process. Therefore, the proposed unsaturated coupled model based on DDA can better describe the hydromechanical response and the stability of the rock slope.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The research is financially supported by the National Natural Science Foundation of China nos. 51779021 and 51478065 and Fund of National Engineering and Research Center for Mountainous Highways no. GSGZJ-2015-09.