Abstract
Displacement discontinuities method (DDM) is convenient for and efficient at dealing with discontinuity problems such as fracture and fault. However, the known method of ease iteration is not available for nonlinear joint surface problems. This paper introduces Barton-Bandis nonlinear joint deformation failure criterion, figures out the propagation model of joint through the maximum energy release rate theory of rock fracture mechanics, and establishes three-dimensional nonlinear displacement discontinuity model for rock joint. This paper gives results of the joint propagation pattern and its distribution law under tension and compression with different parameters and side pressure coefficients via compiled program.
1. Introduction
Displacement discontinuities method (DDM) is an indirect boundary elements method. By inputting relative displacement on fracture surface as unknown variables in seeking internal stress and strain of the model, DDM is convenient for and efficient at dealing with discontinuity problems. Compared to finite element method, DDM requires less variables so as to improve the calculating efficiency and accuracy and has its own advantages on joint propagation simulation.
In 1957, Rongved put forward his solution over a bounded plane area in an infinite solid. In the year 1976 [1], Crouch proposed the solution of plane elasticity problems by the displacement discontinuity method [2, 3]. The method has been developing ever since its birth; Kuriyama and Liu C. L. deduced the exact solutions of the coefficients of the triangular displacement discontinuous element in the three-dimensional condition [4, 5]. Shou developed a quadratic element of DDM and provided detailed analytical formula of the relevant calculus [6]. Wen independently developed DDM under three-dimensional elastic condition and then applied it to the research on dynamic crack researches [7]. Vijayakumar developed a three-dimensional elasticity triangular DDM and provided an adaptive integral method [8].
In terms of joint propagation, Li applied two-dimensional direct and indirect boundary element method to simulate the nonlinear joint propagation [9]. Fotoohi adopted the same method to research nonlinear crack around underground digging zones [10]. Shou adopted three-dimensional hybrid boundary elements to perform a nonlinear analysis on the weak surface around digging zone of a tunnel [11]. Shi wrote software calculating the elastic joint propagation and stress distribution in three-dimensional infinite linear elasticity area [12]. Cheng ran the numerical simulation of hydraulic fracturing for discrete fractured reservoir networks via DDM and graphic theory [13]. Kress established the influence model of crack interference on the height of crack by coupling the three-dimensional displacement discontinuity method in the hydraulic crushing simulator [14]. Cheng proposed a novel linear triangular element of a three-dimensional displacement discontinuity method [15]. Verde modeled injection in a fracture network with DDM [16]. Abdollahipour analyzed the accuracy of higher order displacement discontinuity method [17]. Zhao simulated cracks in 2D piezoelectric semiconductors by DDM [18]. Li analyzed cracks in thermo-magneto-electro-elastic media with DDM [19]. Liu proved that the equations of the DDM for solving 3D crack problems are exactly the same as the equations of the boundary element method [20]. Shen proposed 2D DDM for transversely isotropic materials [21]. Cui proposed contour integral approaches for the evaluation of stress intensity factors using DDM [22]. Xie simulated hydraulic-fracture propagation in unconventional reservoirs [23]. Chen extended DDM to solve cohesive crack propagation [24]. Many other scholars have also done a lot of meaningful research work [25–30].
In the geotechnical engineering, the initial stress state of the rock mass is a pressure state. The deformation characteristics of the discontinuities such as faults show a fairly complex nonlinearity. The simplification of the structural plane as an elastic fracture would ignore its vital character of contact. For rock discontinuity structural plane in compression shear condition, the closed condition of the structural plane should be considered (open, closed without slip, and closed slip). The known method of solving the iterative solution is only available to the linear deformation structural plane. As for the nonlinear structural plane model which conforms to the actual deformation characteristic of the fracture, it is necessary to consider not only the three closed conditions mentioned above, but also the iterative solution of the element nonlinear deformation.
This paper combines the Barton-Bandis normal, tangential deformation model and the Barton shear failure criteria to establish a nonlinear three-dimensional DDM. It invokes the maximum energy release rate theory of rock fracture mechanics to figure out the propagation pattern of a joint [6]. It takes compiled program to research the rock joint propagation process under tension and compression with different parameters and side pressure coefficients.
2. Barton-Bandis Joint Deformation Model
2.1. Normal Deformation Model
Bandis established a hyperbolic function of the natural joint normal fracture closure and effective normal stress through numerous experiments and researches:
where is the maximum closure amount of the joint.
Formulas of the initial stiffness and the maximum closer amount are as follows:
where is the joint roughness coefficient. Joint compress strength () is shear strength of joint surface. is the initial joint aperture under self-weight stress, mm. is uniaxial compressive strength. A, B, C, and D are fitting coefficients.
The normal stiffness formula could be obtained.
2.2. Tangential Deformation Model
The relationship between shear stiffness and shear stress is
where is tangential displacement.
Peak shear stiffness is
where is the maximum shear displacement , and is sample joint length.
2.3. Shear Strength
The Barton-Bandis rule was established on numerous shearing tests on rock joint [31]. It was built on more than 200 artificial sheared joint tests. In geotechnical engineering, the Barton-Bandis model is widely used in analyzing and estimating the shearing strength of a rock joint. It provides a practical method to estimate the shearing strength of a rock joint by JRC and JCS.
where is shear strength, is joint normal stress, and is basic friction angle.
When the normal stress is small and the weathering thickness of joint plane is less than 1 mm, the peak shear strength is controlled by the residual friction angle. The basic friction angle should be replaced by the residual friction angle at this time.
3. DDM Theoretical Solution
The three-dimensional triangular displacement discontinuity element is shown in Figure 1. During the loading process, the relative displacement of the rock joint occurs. As shown in (8), we define the displacement discontinuity on the discontinuous plane as displacement difference between the upper plane (seen as +) and the lower plane (seen as -).
For the triangular displacement discontinuity element, the distributions of the displacement are obtained as follows.
The distributions of the stress are obtained as follows:
where G is the shear modulus, v is Poisson’s ratio, and the kernel function I is as follows.
When the collocating point, being the center of triangular element, is not on the integration element, the integrals in (9) and (10) can be determined using Gauss Quadrature. Otherwise, when the collocating point is on the integration element, the integrals in (9) and (10) become singular since at the collocating point the integrand becomes infinite. This occurs when traction and/or displacement components on the crack are calculated and the contribution from the element that contains the collocating point is needed. This situation cannot be avoided in formation of the governing equations for DD components. We follow the procedure of Zhou [32] for this situation.
The arbitrary discontinuities are separated into N elements. Applying the basic solutions mentioned above, the local coordinates of the joint are replaced by x, y, and z in the above formula:
where , , are tangential and normal face forces of element i, respectively; , , are tangential and normal displacement discontinuities of element j. While , , , , , , , , are the stress influence coefficients, N is the number of elements.
Respectively, with the additional shear force , and the normal force reflecting the existence of real rock joint, the stress of element i can be expressed as
where and are tangential stiffness matrices obtained by element plane deformation model. is normal stiffness matrix obtained by element normal deformation model.
4. Stress Intensity Factors for Fracture
4.1. The Displacement Field and Stress Field in the Crack Front Area
For stress analysis of cracks, the most important part is near the crack front. The research of fracture mechanics is limited to the displacement field and stress field in the vicinity of the crack front. The local coordinate established at the crack front is shown in Figure 2, and the stress field formula expressed with displacement discontinuities on the upper and lower crack surface near the crack front could be obtained:
where , , and are, respectively, the positive stress in direction b and the shear stress in direction bn and bt. , , and are, respectively, the displacement in three coordinate orientations. is the shear deformation modulus. r is the distance in the local coordinate system as shown in Figure 2, and is the angle of local coordinate system as shown in Figure 2.
4.2. Stress Intensity Factor
The definition formulas of stress intensity factor are as follows:
where , , and are, respectively, the intensity factors of Mode I, Mode II, and Mode III cracks.
For the displacement discontinuity in the discrete crack front element in the three-dimensional DDM, it is expressed with , , and , and the above formulas could be further set as follows.
For practical purpose, the limits in (16) can be approximated by simply evaluating the expressions for a fixed value of r. At crack front, we could divide numerous elements whose lengths decrease geometrically, and get a series of approximate values of , , and . Through regression analysis of those approximate values, we could get the stress intensity factors in the crack front. Generally, the exponential function is adopted for fitting. Supposing that , when , , namely, the stress intensity factor of crack front. We could get the stress intensity factor of a series of points around the crack front according to (16). In order to get the value of A, B, and C, a least square fitting of the data pair should be made. Define ; in order to get the minimal value of S, respectively, seek partial derivative of S pairs of A, B, and C, to get the value of A, i.e., the stress intensity factor of crack front.
5. Joint Propagation Criteria
5.1. Joint Propagation Criterion
We adopt the maximum energy release rate (G criterion) as the criteria for 3D compound joint propagation on plane. There are two hypotheses in the criteria: (1) the joint extends along the direction in which the maximum energy release rate is generated; (2) when the energy release rate reaches the critical value in the above determined direction, the joint begins to propagate.
For elastic media, the strain energy release rate is [33] as follows.
5.2. Joint Propagation Direction
According to the energy release rate criterion, the direction of the joint propagation is the direction of the maximum value of the circumferential stress. In the normal plane of the joint frontier, the angle between the joint propagation direction and the normal n-axis, that is, the propagation angle , is as follows.
When , then the brackets take “+”; when , then brackets take the “-”; when , then it is a type I crack and .
Considering type II crack and its effect, we substitute , the effective type stress intensity factor, for in (15) [34]:
where B is g empirical coefficient and it usually takes B = 1.0.
5.3. Joint Propagation Step
According to selected joint propagation criteria, the extended step length of any crack front should be proportional to its equivalent stress intensity factor Ke. Therefore, the extended step length of crack front is
where is the maximum propagation step length for the crack front of step i-1. is the maximum of the energy release rate for all crack fronts in step i. is the maximum propagation step of the crack front of step i. is the energy release rate of the crack front at step i.
The maximum step length of the initial step is related to the initial crack area :
where E stands for the coefficient; this paper takes E = 0.01.
6. Numerical Realization
The simulation of the three-dimensional joint propagation might be distorted. While assuming that the propagation occurs in the joint frontier plane, the crack front could be treated as a set of line segments. After a propagation occurs outside this line, a new set of leading edge segments would be generated [34]. In order to show the shape of the joint, the triangles in Figure 3 are used to perform mesh generation for crack surface between the old and new crack fronts.
(a)
(b)
Joint propagation requires stepwise load. In each step, it needs to determine whether the joint reaches the propagation conditions and calculate the propagated angle and length to complete propagation of the step. If the initial load is too small, the joint would not propagate. The smaller the load step is, the more accurate the simulation would be, while the more calculation would be required. For the hyperbolic nonlinear rock joint model, iterative solution is needed for its nonlinearity. The entire calculation process requires nested iteration to complete, as shown in Figure 4.
The whole process is as follows:(a)At any step of external loading or propagation within a step of the external loading, do a DDM analysis to obtain the displacement discontinuity components (k=1,2,3). Check all elements for being open or closed according to sign of their .(b)If all elements are open, then calculate the strain energy release rate for all front elements and check for extension of the whole joint.(i)If there is any front element that extends, then add new elements to its front; go back to (a) and do a new round DDM analysis since geometry has changed; check all the elements again for their state of open or closed and repeat the process.(ii)If there is no front element that extends, then if the external loading reaches the specified value, stop the simulation; otherwise, increase the external loading; go back to (a) and do a new round DDM analysis since the external loading has changed; check all the elements again for their state of open or closed and repeat the process.(c)If there is any element that is closed, then joint stiffness takes effect.(i)Calculate the joint force for all the closed elements. The normal joint force component is computed as , where is calculated according to (2) with . The shear joint force components are computed as and , where is calculated according to (3). To check slippiness, it is necessary to calculate according to the expression given after (3).(1)If and , then no slip occurs; and are taken as the shear joint forces.(2)If and/or , then slip occurs in the direction in which and value of is taken as the shear joint force component in this direction.(ii)Add the joint force components as outlined above to the traction on all the closed elements, leave the traction on the open elements unchanged, and do a new round DDM analysis to obtain the displacement discontinuity components (k=1,2,3).(iii)Check the tolerance of for all the elements.(1)If for all the elements, then convergent solution for the external loading and geometric conditions is obtained. Propagation criterion is checked for all front elements. If there is any front element that extends, then new elements are added to their front; go back to (a); do a new round DDM analysis and repeat the process. Otherwise, if the external loading reaches the specified value, stop the simulation; if not, the external loading is increased; go back to (a); do a new round DDM analysis and repeat the process.(2)If . for any element, then solution is not convergent. Replace (k=1,2,3) with (k=1,2,3) and check all the elements for being open or closed. It is noted that it is impossible for all the elements to be open, since the external loading and geometrical conditions have not been changed and the case for all elements to be open (that is, the situation listed in (b) above) has been excluded at the present situation. Go back to (i) of (c) and repeat the process.(d)If for any element, then solution is not convergent. Formula (16) would be used to calculate , the energy release rate of the crack front element. If all the calculated are less than or equal to , (that is, fractures do not propagate), apply load and go back to (a) until fractures propagate and new crack front element appears.(e)When new crack front element appears, we need to determine whether the calculation of joint propagation should continue. If it is still required, update statistics of crack front element, and go back to (a). If not, terminate the calculation.
We use compiled software to calculate and analyze conditions in reference and have got the following comparison between result of joint propagation calculation and result in the reference, as shown in Figure 5 [35]. The calculation results show that the software can simulate engineering problems correctly.
(a) Simulation result
(b) Experiment result
7. Parameter Analysis
Although the actual joint in the rock is often a plane, its geometric boundaries would be normally irregular shapes. In order to facilitate the analysis of the universality, this paper takes the discoid joint as an example to carry out the simulation. The disc model joint angle is 45°, radius 100m, and shape as in Figure 6 in step0. The mechanical parameters of rock and joint are shown in Table 1.
7.1. Propagation Process under Compression
The three-dimensional propagation under vertical compression 100MPa is shown in Figure 6, and the whole process is extended into 11 steps. Step0 is the initial joint shape.
The joint propagation is carried out in the direction of the vertical pressure, and the propagation begins at the upper and lower ends of the joints. They, respectively, propagate upward and downward. The propagations at the upper and lower sides are the largest, and they decline to the middle, forming a wrap-like joint propagation. At the end of the propagation, variation of stress intensity factor of the joint front slows down and the joint propagation tends to be stable due to the three-dimensional wrap-like propagation shape.
7.2. Propagation Process under Tension
The three-dimensional propagation process under the vertical of 100 MPa could be seen in Figure 7. The process proceeds for 11 steps, and step0 is the initial joint shape.
The joint propagation process under tension differs from that under compression. Although both of them begin at the upper and lower ends, the upper end propagates downward and the lower end propagates upward. After three steps of propagation, the direction of joint propagation begins to incline to the horizontal direction, that is, perpendicular to the direction of tension. Like the compressing propagation, the final type presents warped shape and is eventually stabilized.
7.3. Comparison of the Propagation with Different Cohesion c
Joint cohesion c is usually small. When keeping other parameters’ values the same as those in Table 1, we simulate joint propagation under vertical compression of 100MPa with different joint c values of 0.05MPa, 1MPa, and 2MPa. The results of the simulation are shown in Figure 8.
The joint propagation peaks when c value is 0.05MPa, while it changes slightly when c value is 1MPa and 2MPa. It indicates that when c value increases to some extent, joint propagation would be limited, and the main differences of joint propagation would be reflected in propagation form. The outer edge shape of the upper and lower ends of the joint are basically the same, yet the propagation of the middle part is different. With the increase of c value, the propagation in the middle part in the vertical pressure direction is restricted, and the overall propagation of the joint tends to be stable.
7.4. Comparison of the Propagation with Different Joint Friction Angles φ
Keeping other parameters’ values the same as those in Table 1, we change the joint friction angle φ into 5°, 10°, 15°, 20°, 25°, 30°, 35°, and 40°. The results of joint propagation under vertical compression of 100MPa are shown in Figure 9.
The propagation distances change a little when φ value differs. The difference is mainly in the propagation shape. In general, with the increase of φ value, the propagation in the middle part in the direction of vertical pressure gradually decreases. The propagation form tends to be a regular pattern whose middle part is tight and two ends are loose. When the φ value is small (5° and 10°), the joint propagation pattern is irregular, and prominent propagations appear in some individual points at the upper and lower ends. The reason of that situation is that as the joint mechanical property is weak and the stress is more intense, the stress concentration coefficient of the individual points would be large.
7.5. Comparison of Different Lateral Pressure Coefficients
Using the parameters in Table 1, we simulate joint propagation under vertical compression of 100MPa with different lateral pressure coefficient λ of 0, 0.1, and 0.5 on the Y-axis. With the increase of the lateral pressure coefficient on the Y-axis, the propagation of the joint obviously tightens in the Y direction, and the horizontal propagation of the middle part is constrained, as shown in Figure 10.
When the lateral pressure coefficient on the X-axis is 0, 0.1, and 0.2, the propagation results under vertical compression of 100MPa are shown in Figure 11. With the increase of the X-side lateral pressure coefficient, the joint normal pressure increases, and the propagation under vertical compression is limited.
8. Conclusions
This paper introduces the Barton-Bandis nonlinear joint deformation failure criterion, finds out joint propagation pattern through stress intensity factor parameter in rock fracture mechanics, and establishes the three-dimensional nonlinear joint displacement discontinuity method model. We compare different propagations under tension and compression with different parameters and lateral pressure coefficients through program and get the following results.
The propagation results with vertical compression and tension shows that the direction of joint propagates is along the direction of the maximum pressure which is perpendicular to the direction of tension:(1)The propagation of the three-dimensional joint is the greatest at both ends of the pressure direction. The middle part gradually reduced into warped shape, and propagation tends to be stable.(2)The cohesion c has a great influence on the range of joint propagation. The larger the c value is, the smaller the range of joint propagation would be.(3)The effect of friction angle φ on joint propagation is mainly reflected in the propagated shape. When the value of φ is small, the propagation would be disordered, and it would be prominent in some individual points. With the increase of φ value, the propagation turns to be regular and would be in an obvious warped shape.(4)With the increase of the lateral pressure coefficient, the joint propagation range in the corresponding direction gets obviously smaller. The magnitude of the lateral pressure, especially that of joint normal pressure, plays a vital role in joint propagation.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is financially supported by National Key R&D Program of China (2017YFC0806000), National Natural Science Foundation of China (41601574), and Chongqing Basic and Frontier Research Project (cstc2015jcyjBX0118).