Abstract
In this paper, a new approach has been developed for predicting the hydraulic and mechanical relationship of individual fractures subjected to normal stress and compression-shear stress. Considering that the closure process of rough fracture subjected to normal stress can be divided into two phases (linear behavior and nonlinear behavior), a relationship between normal stress and fracture aperture is derived through the minimum potential energy principle. Then, a formulation for calculating fracture permeability during shearing and compression processes is developed. Furthermore, a formulation for determining fracture aperture during the crack growth process is obtained, which is further implanted into the permeability model to predict the hydraulic behavior of fractured rock during fracture propagation. This new model not only considers the normal deformation of the fracture but also, and more importantly, integrates the effect of fracture propagation and shear dilation. Theoretical studies demonstrate that fracture permeability increases nonlinearly during fracture propagation. At last, experimental results and analytic results are compared to demonstrate the usefulness of the proposed models, and satisfactory agreements are obtained.
1. Introduction
The coupling between hydraulic and mechanical behavior is of considerable interest in several areas of rock engineering, such as underground tunnels [1], naturally fractured petroleum reservoirs [2], rock slope stability, and underground disposal of radioactive waste. Fractures of rock masses are key factors governing the coupling between hydraulic and mechanical behavior [3]. For engineering applications, an improved understanding of the stress-flow coupling of individual fracture, which is fundamental to the hydromechanical analysis of fractured rock, is needed. However, to accurately describe the hydraulic conductivity of fractured rock masses under complex stress is an issue that is not completely resolved. This difficulty mainly comes from that the heterogeneous nature of rock materials contains various natural fractures of different scales [4–6]. Consequently, many investigators have been focusing on the hydraulic and mechanical properties of fractures by conducting experimental tests and analytical research [7–9].
Previous studies of fluid flow through fractured rock have mainly focused on the effects of hydromechanical coupling, fracture characteristics (such as aperture distribution, surface roughness, and contact area), or a combination of both effects [10–12]. It is examined by a multitude of experiments that aperture distribution, surface roughness, and contact area have a significant influence on flow behaviors in fractures [13–16], and all of them are established with a function of normal or shear stress. Laboratory investigations on fractured rocks show that shear deformation can induce fracture opening, which is recognized as shear dilatancy [17]. After recognizing the dilation phenomenon in shear processes, shear dilatancy and related hydromechanical behavior of fractured rocks were studied extensively. Yeo et al. [18] investigated the effect of shear displacement on the aperture and permeability of rock fracture by using radial and unidirectional flow test and numerical simulation. Chen et al. [19] presented a strain-dependent formulation to address the change in hydraulic conductivity resulted from normal and shear loads. Wei and Liu [20] presented a fractal-based model for fracture deformation under shearing and compression. Although considerable progress has been made in the past four decades, models for describing the hydromechanic behavior with different conditions should be further studied [21–23].
The existing fractures in rock mass may close, open, or grow, and new cracks may be induced under external loads [8, 19, 24–27]. During the process of fracture propagation, the length and aperture of fracture, which is widely regarded as factors influencing the flow properties of the rock mass, will change. In all the above models for describing fracture hydromechanical properties, however, only the obvious factors such as aperture distribution, surface geometry, and contact area are well investigated; the variation of fracture length induced by fracture propagation is not taken into account. In practice, the effect of fracture propagation on fluid flow was not ignored, especially when the rock mass suffers relatively large external loads. Specifically, a variation of fracture length will change the average fracture aperture and flow path. However, to the best knowledge of the authors, limited literature for investigating fracture hydromechanical properties has considered the effect of propagation on fracture characteristics.
Therefore, in this paper, we extend the above work and develop a new approach to establish the hydromechanical behavior for fractured rock under shear and normal stress. A simple review of the continuously yielding fracture model proposed by Cundall and Hart [28] is first introduced, and then, an improved relationship for shear stress and displacement is built with the consideration of mechanical properties on fracture surface roughness. Based on the principle of minimum potential energy, the relationships between stress and deformation of open fracture are developed for the compression process and shearing-compression process, respectively. And further considering fracture propagation under compression-shear, a formulation for calculating fracture aperture in the crack growth region is obtained. Unlike most previous models, the new model gives particular emphasis to the effect of fracture propagation and shear dilation.
2. Materials and Methods
The fluid flow through fracture of rock masses has commonly been described by the cubic law. It is derived from the assumption that the incompressible and laminar fluid flow occurs between two smooth-parallel plates. It is often denoted as [29]. where is the flowrate, is the fracture aperture, is the water head difference through fracture surface, and is the constant related to gravitational acceleration, flow viscosity coefficient, and the shape and size of the fracture surface. Natural rock fractures, however, are featured with rough surfaces and variable aperture, which conflicts with the assumption. Many scholars have been focused on the studying of flow characteristics of rough fractures and proposed various modified cubic law. Neuzil and Tracy [29] put forward a probability density function to modify fracture width based on the assumption that width changes in the vertical direction and keeps a constant in a parallel direction. Walsh [15] built a modifier formula containing a ratio of area contact by imitating the heat conduction theory. To account for the deviation of the nature fractures from the ideally smooth plates, Chen et al. [19] presented an equivalent formulation for fracture permeability by adopting a dimensionless constant, ; it can be expressed by where is the fracture permeability. The aperture at any fracture deformation is given as [19]: where and are the initial aperture and the change in aperture due to external loads.
3. Relationship between Shear Stress and Shear Displacement
Previous studies demonstrate that rock joints show high nonlinear shear properties under direct shear stress. It is reflected in the features that normal stiffness varies with normal stress, and internal friction angle decreases with increasing shear displacement. Cundall and Hart [28] proposed a continuously yielding joint model (referred to as the CY model). This model was intended to produce nonlinear behavior and shear dilation and could be applied to analyze the internal mechanism of progressive damage of rock joints under shear stress. Based on the CY model, a new relationship between shear stress and shear displacement is built in this paper.
A demonstration of the CY model with a typical shear stress-displacement curve, normal displacement curve, and the bounding shear strength curve is shown in Figure 1. The response to shear stress in the CY model is calculated incrementally as [28] where is the shear displacement increment, is shear stiffness, and is the elastic recoverable ratio, which depends on the distance from the actual shear stress curve to bounding strength curve, (Figure 1). Theoretically, is written as [31]
The factor, , is used to reproduce the speed in approaching the bounding stress curve: a bigger indicates a lager speed to the bounding stress. In practice, is limited to 0.75 in order to avoid numerical noise. To simplify the equation, is set to be zero in this paper.
To analyze the effect of roughness in fractured rock, Oh [30] employed a formulation to substitute the bounding strength in the CY model with where , , and are normal stress, residual friction angle, and dilatancy angle, respectively.
On the basis of test data, Barton et al. [14] proposed an empirical expression of dilatancy angle with the joint roughness coefficient (JRC), which was defined as where JCS refers to the fracture surface compressive strength, subscript mob denotes that parameters change with the shear displacement during the process of shearing. Tse and Curden [32] proposed an empirical formulation to calculate JRC, and it would get the maximum value at the point of peak shear strength.
Combining Eqs. (4) through (7) yields
Combining with Eqs. (4) and (5) and using the initial condition “when , ” gives
Eq. (8) provides a new relationship between shear stress and displacement, which is highly nonlinear and is related to JRC. Previous studies focused on fracture shear properties are mainly suggested in an empirical manner and fail to explicitly describe the influence of mechanical properties on fracture surface roughness. Eq. (8) responds to our effort to consider the impact of fracture roughness during the process of shearing. Note that due to the failure of asperities, the dilation angle and JRC decrease with plastic shear displacement, which is aligned with fracture deformation.
4. Behavior of Fracture Deformation
In this section, we will discuss the hydromechanical behavior of a single fracture under different conditions, including normal stress, compression-shear stress, and fracture propagation process.
4.1. Fracture Closure under Normal Stress
As reviewed in the introduction section, for fractured rock, the relationship between normal displacement and normal stress is nonlinear. The deformation response of an open fracture subject to normal stress may include two phases: one is linear closure of the fracture, and the other is nonlinear closure due to contact areas. The observed nonlinear stress-deformation behavior of fracture could be attributed to the increasing contact areas and crushing of asperities. In this subsection, fracture closure due to the normal stress is investigated through the minimum potential energy principle.
An infinite plate with rough fractured surfaces under far-field compression stress is shown in Figure 2. With the working of normal stress, fracture deformation occurs, and new elastic energy is produced. Then, the energy of the fracture element is [33] where is the initial elastic energy stored in the fracture element. The new elastic energy, , induced by normal displacement, corresponds to the mode-I crack field. The work applied by the normal stress can be expressed as [33] where and are the work and displacement, respectively, and refers to displacement function. The region of integration, , depends on the fracture length . The effective working force, , is composed of at least two major parts, which correspond to linear closure and nonlinear closure. For the first phase, there are no contact areas in the fracture; thus, is the effective stress doing work. When normal stress increases to closure stress, , a portion of stresses are counteracted due to surfaces contacting and asperities crushing. A reduction factor is used to account for the effect of surfaces contacting and asperities crushing [34]. Then, can be determined by
Note that with normal stress increasing, the growth rate of effective working stress decreases because of the increasing contact areas and the number of contacts. Therefore, the factor should be dependent on effective stress.
The Griffith theory made it possible for the first time to like within the limits of the thermodynamic approach the process of failure of a solid with the energy transformations accompanying the formation of new fracture surfaces [35]. Griffith [36] argued that the value of deformation in elliptic fracture under compression is maximal in the middle and decreases to both ends. Therefore, the displacement function is assumed as where , , and are the half fracture length, correction coefficient of the fracture shape, and normal displacement, respectively. Substituting Eq. (13) into Eq. (11), and integrating, we obtain where is the width of the material. According to Ashby and Hallam [33], new elastic energy accumulated by is where and are the constant and fracture elasticity modulus, respectively.
The minimum potential energy principle for a conservative elastic system states that of all geometrically admissible displacement fields, that which brings the loaded structure into elastic equilibrium also minimizes the potential energy functional [37]. Minimizing the energy of element during the deformation process, hence
Applying Eqs. (10)–(15) into Eq. (16), and ignoring the changes of , yields
Hence, a theoretical model of fracture closure for rough fracture is formulated, which is totally determined by effective working stress and a set of parameters. It is to be noted that for the first phase () Eq. (17) can be rewritten as Hook’s law by defining as the fracture stiffness. where is the normal strain and is fracture stiffness.
4.2. Fracture Aperture under Compression-Shear Stress
A fracture subject to compression-shear stress will slip and dilate to facilitate the opening of fracture. As illustrated in Figure 3(a), the normal stress and shear stress of fracture surface are computed by where is the inclination measured in terms of the angle from the principal axes.
(a)
(b)
4.2.1. Fracture Aperture before Propagation
Considering a fracture to be embedded into a rock subject to compression-shear stress, a similar method and notations in the previous subsection are used to investigate the fracture deformation. Then, the energy of fracture element subject to compression-shear stress is given as where is new elastic energy induced by shear displacement and can be calculated by [33] where refers to shear displacement. One can observe from Figure 3(b) that the working applied by the shear load is positive and applied by the normal load is negative due to shear dilation. Then, the whole acting, , applied by and can be expressed as where is the shear displacement function and can be similarly given as follows
Substituting Eqs. (13) and (23) into Eq. (22) yields
The relationship between the total fracture displacement () and or can be written as where is the dilation angle. Applying Eqs. (13), (15), and (20)–(25) into Eq. (16) yields
As discussed before, the fractured rock under shear stress will dilate to facilitate the opening of fracture, which means few stresses are counteracted during the dilation process. Therefore, the function of working stress is set to . Combining Eqs. (25) and (26), the change of aperture subject to compression-shear stress is
Combining Eqs. (2), (3), and (27), the stress-dependent formulation of fracture permeability is determined by
Eq. (9) presents a relation between shear stress and shear displacement . Accordingly, a formulation of strain-dependent permeability is developed by combining Eq. (9) and (28). It is important to note that our Eq. (28) is applicable when the rock dilates. Other mechanic deformations can be similarly analyzed with this method.
4.2.2. Fracture Aperture during Propagation
Eq. (27) describes fracture aperture change during compression-shear stress without considering the effect of fracture propagation. In a general case, fractures under high stress involve aperture changes in the crack growth region. As illustrated in Figure 4(a), when stress (or) increases to the initial cracking value, stress-induced cracks are produced. Assuming that stress-induced cracks are straight lines and parallel to (Figure 4(a)), then the normal stress of wing crack, , equals to . A number of experimental results and theoretical calculations demonstrate that the initial cracking is in accordance with mode-I under compression-shear [24, 38, 39].
(a)
(b)
Similar to Eq. (22), the total work done by , , and is where is the displacement of the wing cracks. The integral region, , depends on the wing crack length . Martin [40] provided a general equation for estimating the length of the crack extension once the cracks initiate; it was written as: where , , and are the curve-fitting parameters.
In order to illuminate the fracture deformation, the system (original crack and wing crack pair) is replaced by an equivalent single straight crack [33]. Similar to Eqs. (13) and (23), the following relationships are formulated by substituting for
Assuming that the displacement of the wing crack, , is proportional to the distance from the wing crack tip, then it can be defined by where is the constant. Integrating Eq. (29) with Eqs. (31) and (32), is derived to be where . By employing the Taylor series expansion, Eq. (33) can be adequately approximated by
Note that the value in the bracket of the first item approximately equals to 1 for ; then, Eq. (34) is reduced to the following form:
Combining Eqs. (15), (16), (20), (21), and (35) yields
Therefore, the aperture and fracture permeability during the fracture propagation process can be expressed as
Eq. (38) made it possible for the first time to predict the permeability evolution of rock during fracture propagation.
5. Model Validation
In this section, we mainly focused on evaluating the usefulness of the proposed models against the experimental observations from the literature.
5.1. Validation of the Proposed Model under Normal Load
Many investigators have explored the relationship between fracture aperture and normal stress by conducting experiments. Two experimental datasets from Ref. [41, 42] are selected to validate the model proposed in Section 4.1.
The information on the datasets adopted in this paper is shown in Table 1. The values for fracture length and initial aperture are directly determined from laboratory tests. Fracture elastic modulus can be determined from the stress-strain curve in the compression test; in this paper, is set to 4 GPa for both fractured rocks. The determination of threshold stress , for intact rock, may utilize the method of fracture volume strain proposed by Martin and Chandler [43]. For separated rock, is estimated from the closure curves by identifying the turning point from linear behavior to nonlinear behavior. Moreover, an empirical function of the factor is proposed based on the experimental data in Figure 5, in which is a parameter related to contact areas and asperity degradation. The remaining parameter values ( and ) are identified by fitting Eq. (17) with the measured data in Figure 5.
(a)
(b)
The comparisons between the observed closure curves and theoretical results from Eq. (17) are depicted in Figure 5, which manifests that the mechanical deformation of the fractured rock under normal stress is well reproduced by the proposed model. It should be noted that for the two experimental datasets fitting, the values of the correlation coefficient () are above 0.95. This confirms that the proposed model is capable of reproducing the change of fracture aperture during the normal loading process.
5.2. Validation of the Proposed Model under Compression-Shear Load
Esaki et al. [44] conducted a comprehensive study regarding shear deformation and dilatancy dependence of hydraulic conductivity for rock fractures through laboratory experiments. The coupled shear-flow tests were conducted on granite samples with a size of 120 mm in length, 100 mm in width, and 80 mm in height. Note that the tested sample has been split into two parts to create an artificial fracture. Eq. (27) should be applied for calculating fracture deformation rather than Eq. (37). The values of the parameters used for matching the data of Esaki et al. [44] are listed in Table 2. The initial aperture and fracture length were directly obtained from Ref. [44]. He also provided values for JRC and JRS, which are 9 and 162 MPa, respectively. The dilation angle, , is identified as a constant for simplicity. In this case, its value can be determined from Eq. (7). The bounding strength can be calculated using Eq. (6), but the better alternative is to directly identify it from the tests. The rest of the parameter values are determined by fitting the experimental curves using relevant equations. It should be noted that the coupled shear-flow tests were conducted for one cycle of forwarding and reverse shearing, and only the forward experimental results are adopted in this paper. The results from the reverse shearing process can be also matched in the same method.
The comparisons between the aperture analytically predicted by using the proposed model given in Eq. (27) and those measured from coupled shear-flow tests are shown in Figure 6. As can be seen in Figure 6, for different normal stresses, the observed shear dilatancy behavior is fairly well captured by Eq. (27) with a single group of parameter values. Slight discrepancies between the experimental curves and analytical curves are observed when the normal stress is small. This can be explained by the fact that the dilation angle undergoes negative decay under small normal stress.
(a)
(b)
(c)
(d)
5.3. Hydraulic Behavior during Crack Propagation Process
Chen et al. [45] performed several triaxial compression tests on granite and found that the evolution of rock permeability experienced a clear decrease phase in the initial microcrack closure region and a dramatic permeability increase phase in the crack growth region. In this subsection, the observed data in the crack growth region are compared with the predicted results. However, the hydraulic behavior with multiple sets of fractures is complex and obviously distinguished from that of a single fracture due to varied aperture values and sizes in each fracture. Min et al. [46] proposed the concepts of “equivalent frequency” and “equivalent aperture” to approximate the complex fracture deformation. Using this method, the fracture permeability with multiple sets of fractures may be formulated in a simple form as follows: where the subscript eq refers to demote the equivalent concept, , and the stress-induced crack length, , can be obtained using Eq. (30).
For simplicity, two-dimensional analysis is employed here, and the associative flow rule is used for modeling the dilatancy of fracture. Assuming that the equivalent fracture is orientated at the angle of failure surface (i.e., ), the stresses ( and ) of fracture surface can be determined by Eq. (19). Since there is no model of plasticity, the dilatancy angle is assumed equal to the friction angle. The parameter values identified in the model are summarized in Table 3, in which the initial aperture, elastic modulus, and friction angle are directly obtained from Ref. [45]. The remaining parameters are obtained by fitting the measured data in Figure 7 using Eq. (39). Note that , , and in Eq. (39) are set to 1 to reduce the number of fitting parameters. The curves illustrated in Figure 7 depict that the evolution of permeability in the crack growth region is fairly well captured by the analytic model, supporting the usefulness of our approach.
6. Conclusions
The hydromechanical behavior of fractured rock is an important issue for practical applications. This paper develops the relationships between hydraulic and mechanical behavior of individual fracture subjects to normal and shear loads. Based on this study, the following conclusions can be drawn: (1)The formulas for describing the hydraulic and mechanical relationship in the fractured rock are derived through the minimum potential energy principle, including the cases of fractured rock under compression, under compression, and shearing. The proposed approach enables easier fracture deformation analysis to include the effect of shear dilation and is validated by favorable comparisons between measured and analytical results(2)On the basis of the minimum potential energy principle and cubic law, a permeability model with consideration of the fracture propagation is derived. Theoretical results show that the permeability increases nonlinearly during fracture propagation, and the proposed model is capable of predicting the permeability evolution during fracture propagation
Data Availability
The data is available on request by contacting the corresponding author. [email protected]
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (No. 52064006 and No. 52004072), the Research Fund for Talents of Guizhou University (Grant No. 201901), Specialized Research Funds of Guizhou University (Grant No. 201903), and Guizhou joint support project ([2018]2787, [2020]4Y044).