# Three-Dimensional Propagation Simulation and Parameter Analysis of Rock Joint with Displacement Discontinuity Method

Authors: Ke Li; Xinghong Jiang; Hao Ding; Xuebing Hu
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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 hissolution 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δn and effective normal stress σn through numerous experiments and researches:(1)δn=σnδnmaxKniδnmax+σnwhereδnmax is the maximum closure amount of the joint.Formulas of the initial stiffnessKni and the maximum closer amount are as follows:(2)Kni=-7.15+1.75JRC+0.02JCSaj(3)δnmax=A+BJRC+CJCSajDwhereJRCis the joint roughness coefficient. Joint compress strength (JCS) is shear strength of joint surface. aj=JRC/50.2σc/JCS-0.1 is the initial joint aperture under self-weight stress, mm. σc is uniaxial compressive strength.A,B,C, andD are fitting coefficients.The normal stiffnessKn formula could be obtained.(4)Kn=Kni1-δn/δnmax2


### 2.2. Tangential Deformation Model

The relationship between shear stiffnessKs and shear stress τp is(5)Ks=τpδswhereδs is tangential displacement.Peak shear stiffnessKsp is(6)Ksp=τpδspwhereδsp is the maximum shear displacement δsp=Ln/500JRC/Ln0.33, and Ln 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 byJRC andJCS.(7)τ=σntan⁡JRClog10⁡JCSσn+φbwhereτ is shear strength, σn is joint normal stress, and φb 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φb should be replaced by the residual friction angle φb at this time.


## 2.1. Normal Deformation Model

Bandis established a hyperbolic function of the natural joint normal fracture closureδn and effective normal stress σn through numerous experiments and researches:(1)δn=σnδnmaxKniδnmax+σnwhereδnmax is the maximum closure amount of the joint.Formulas of the initial stiffnessKni and the maximum closer amount are as follows:(2)Kni=-7.15+1.75JRC+0.02JCSaj(3)δnmax=A+BJRC+CJCSajDwhereJRCis the joint roughness coefficient. Joint compress strength (JCS) is shear strength of joint surface. aj=JRC/50.2σc/JCS-0.1 is the initial joint aperture under self-weight stress, mm. σc is uniaxial compressive strength.A,B,C, andD are fitting coefficients.The normal stiffnessKn formula could be obtained.(4)Kn=Kni1-δn/δnmax2


## 2.2. Tangential Deformation Model

The relationship between shear stiffnessKs and shear stress τp is(5)Ks=τpδswhereδs is tangential displacement.Peak shear stiffnessKsp is(6)Ksp=τpδspwhereδsp is the maximum shear displacement δsp=Ln/500JRC/Ln0.33, and Ln 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 byJRC andJCS.(7)τ=σntan⁡JRClog10⁡JCSσn+φbwhereτ is shear strength, σn is joint normal stress, and φb 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φb should be replaced by the residual friction angle φb at this time.


## 3. DDM Theoretical Solution

The three-dimensional triangular displacement discontinuity element is shown in Figure1. 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 -).(8)Dxx,y,0=uxx,y,0--uxx,y,0+Dyx,y,0=uyx,y,0--uzx,y,0+Dyx,y,0=uzx,y,0--uzx,y,0+Figure 1
3D displacement discontinuity element.For the triangular displacement discontinuity element, the distributions of the displacement are obtained as follows.(9)ux=Dx21-νI,z-zI,xx+Dy-zI,xy+Dz-1-2νI,x-zI,xzuy=Dx-zI,xy+Dy21-νI,z-zf,yy+Dz-1-2νI,y-zI,yzuz=Dx1-2νI,x+zI,xz+Dy1-2νI,y+zI,yz+Dz21-νI,z-zI,zzThe distributions of the stress are obtained as follows:(10)σxx=2GDx2I,xz-zI,xxx+Dy2νI,yz-zI,xxy+DzI,zz+1-2νI,yy-zI,xxzσyy=2GDx2νI,xz-zI,xyy+Dy2I,yz-zI,yyy+DzI,zz+1-2νI,xx-zI,yyzσzz=2GDx-zI,xzz+Dy-zI,yzz+DzI,zz-zI,zzzτxy=2GDx1-νI,yz-zI,xxy+Dy1-νI,xz-zI,xyy+Dz-1-2νI,xy-zI,xyzτyz=2GDx-νI,xy-zI,xyz+DyI,zz+νI,xx-zI,yyz+Dz-zI,yzzτzx=2GDxI,zz+νI,yy-zI,xxz+Dy-νI,xy-zI,xyz+Dz-zI,xzzwhere G is the shear modulus, v is Poisson’s ratio, and the kernel function I is as follows.(11)Ix,y,z=18π1-ν∬D1x-ξ2+y-η2+z2dξdη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:(12)S1i=∑j=1NA11ijD1j+∑j=1NA12ijD2j+∑j=1NA13ijD3jS2i=∑j=1NA21ijD1j+∑j=1NA22ijD2j+∑j=1NA23ijD3jS3i=∑j=1NA31ijD1j+∑j=1NA32ijD2j+∑j=1NA33ijD3jwhereS1i, S2i, S3i are tangential and normal face forces of element i, respectively; D1j, D2j, D3j are tangential and normal displacement discontinuities of element j. While A11ij, A12ij, A13ij, A21ij, A22ij, A23ij, A31ij, A32ij, A33ij are the stress influence coefficients, N is the number of elements.Respectively, with the additional shear forceK1iD1i, K2iD2i and the normal force K3iD3i reflecting the existence of real rock joint, the stress of element i can be expressed as(13)K1iD1i+S1i=∑j=1NA11ijD1j+∑j=1NA12ijD2j+∑j=1NA13ijD3jK2iD2i+S2i=∑j=1NA21ijD1j+∑j=1NA22ijD2j+∑j=1NA23ijD3jK3iD3i+S3i=∑j=1NA31ijD1j+∑j=1NA32ijD2j+∑j=1NA33ijD3jwhereK1i and K2i are tangential stiffness matrices obtained by element plane deformation model. K3i 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 Figure2, and the stress field formula expressed with displacement discontinuities on the upper and lower crack surface near the crack front could be obtained:(14)σbbr,θ=μ41-νrunθ=π-unθ=-πsin⁡θ2cos⁡θ2cos⁡3θ2+ubθ=π-ubθ=-πcos⁡θ21+sin⁡θ2sin⁡3θ2σbnr,θ=μ41-νrunθ=π-unθ=-πcos⁡θ21-sin⁡θ2sin⁡3θ2-ubθ=π-ubθ=-πsin⁡θ2cos⁡θ2cos⁡3θ2σbtr,θ=μ4rutθ=π-utθ=-πcos⁡θ2Figure 2
Local coordinate at the crack front.whereσbb, σbn, and σbt are, respectively, the positive stress in direction b and the shear stress in direction bn and bt. un, ub, and ut 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:(15)KI=limr→0⁡2πrσbbr,θθ=0KII=limr→0⁡2πrσbnr,θθ=0KIII=limr→0⁡2πrσbtr,θθ=0whereKI, KII, and KIII 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 withDb(r), Dn(r), and Dt(r), and the above formulas could be further set as follows.(16)KΙ=-μ41-νlimr→0⁡2πrDbrKΙΙ=-μ41-νlimr→0⁡2πrDnrKΙΙΙ=μ4limr→0⁡2πrDtrFor practical purpose, the limits in (16) can be approximated by simply evaluating the expressions for a fixed value ofr. At crack front, we could divide numerous elements whose lengths decrease geometrically, and get a series of approximate values of KI, KII, and KIII. 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 K⌢=A-BrC, when r=0, K⌢=A, namely, the stress intensity factor of crack front. We could get the stress intensity factor Ki 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 (ri,Ki) should be made. Define S=ΣK⌢i-Ki2=ΣA-BriC-Ki2; in order to get the minimal value of S, respectively, seek partial derivative of S pairs of A, B, and C, ∂S/∂A=∂S/∂B=∂S/∂C=0 to get the value of A, i.e., the stress intensity factor of crack front.


## 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 Figure2, and the stress field formula expressed with displacement discontinuities on the upper and lower crack surface near the crack front could be obtained:(14)σbbr,θ=μ41-νrunθ=π-unθ=-πsin⁡θ2cos⁡θ2cos⁡3θ2+ubθ=π-ubθ=-πcos⁡θ21+sin⁡θ2sin⁡3θ2σbnr,θ=μ41-νrunθ=π-unθ=-πcos⁡θ21-sin⁡θ2sin⁡3θ2-ubθ=π-ubθ=-πsin⁡θ2cos⁡θ2cos⁡3θ2σbtr,θ=μ4rutθ=π-utθ=-πcos⁡θ2Figure 2
Local coordinate at the crack front.whereσbb, σbn, and σbt are, respectively, the positive stress in direction b and the shear stress in direction bn and bt. un, ub, and ut 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:(15)KI=limr→0⁡2πrσbbr,θθ=0KII=limr→0⁡2πrσbnr,θθ=0KIII=limr→0⁡2πrσbtr,θθ=0whereKI, KII, and KIII 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 withDb(r), Dn(r), and Dt(r), and the above formulas could be further set as follows.(16)KΙ=-μ41-νlimr→0⁡2πrDbrKΙΙ=-μ41-νlimr→0⁡2πrDnrKΙΙΙ=μ4limr→0⁡2πrDtrFor practical purpose, the limits in (16) can be approximated by simply evaluating the expressions for a fixed value ofr. At crack front, we could divide numerous elements whose lengths decrease geometrically, and get a series of approximate values of KI, KII, and KIII. 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 K⌢=A-BrC, when r=0, K⌢=A, namely, the stress intensity factor of crack front. We could get the stress intensity factor Ki 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 (ri,Ki) should be made. Define S=ΣK⌢i-Ki2=ΣA-BriC-Ki2; in order to get the minimal value of S, respectively, seek partial derivative of S pairs of A, B, and C, ∂S/∂A=∂S/∂B=∂S/∂C=0 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 valueGc in the above determined direction, the joint begins to propagate.For elastic media, the strain energy release rate is [33] as follows.(17)G=1-ν2EKΙ+1-ν2EKΙΙ+1+νEKΙΙΙ


### 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θ0, is as follows.(18)tan⁡θ02=14KΙKΙΙ±KΙ2KΙΙ2+8WhenKΙΙ<0, then the brackets take “+”; when KΙΙ>0, then brackets take the “-”; when KΙΙ=0, then it is a type I crack and θ0=0.Considering type II crack and its effect, we substituteKΙe, the effective type stress intensity factor, for KΙ in (15) [34]:(19)KΙe=KΙ+BKΙΙΙ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(20)Δai=Δamaxi-1GiGmaxiwhereΔamaxi-1 is the maximum propagation step length for the crack front of step i-1. Gmaxi is the maximum of the energy release rate for all crack fronts in step i. Δai is the maximum propagation step of the crack front of step i. Gi is the energy release rate of the crack front at step i.The maximum step length of the initial stepΔamax0 is related to the initial crack area S0:(21)Δamax0=ES0where E stands for the coefficient; this paper takes E = 0.01.


## 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 valueGc in the above determined direction, the joint begins to propagate.For elastic media, the strain energy release rate is [33] as follows.(17)G=1-ν2EKΙ+1-ν2EKΙΙ+1+νEKΙΙΙ


## 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θ0, is as follows.(18)tan⁡θ02=14KΙKΙΙ±KΙ2KΙΙ2+8WhenKΙΙ<0, then the brackets take “+”; when KΙΙ>0, then brackets take the “-”; when KΙΙ=0, then it is a type I crack and θ0=0.Considering type II crack and its effect, we substituteKΙe, the effective type stress intensity factor, for KΙ in (15) [34]:(19)KΙe=KΙ+BKΙΙΙ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(20)Δai=Δamaxi-1GiGmaxiwhereΔamaxi-1 is the maximum propagation step length for the crack front of step i-1. Gmaxi is the maximum of the energy release rate for all crack fronts in step i. Δai is the maximum propagation step of the crack front of step i. Gi is the energy release rate of the crack front at step i.The maximum step length of the initial stepΔamax0 is related to the initial crack area S0:(21)Δamax0=ES0where 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.Figure 3
The new and old joint front edge surface mesh.
(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 Figure4.Figure 4
Flowchart of joint propagation.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 componentsDki (k=1,2,3). Check all elements for being open or closed according to sign of their D3i.(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 forceK·D for all the closed elements. The normal joint force component is computed as Kn·D3i, where Kn is calculated according to (2) with δn=D3i. The shear joint force components are computed as Ks·D1i and Ks·D2i, where Ks is calculated according to (3). To check slippiness, it is necessary to calculate τp according to the expression given after (3).
(1)
IfKs·D1i<τp and Ks·D2i<τp, then no slip occurs; Ks·D1i and Ks·D2i are taken as the shear joint forces.(2)
IfKs·D1i>τp and/or Ks·D2i>τp, then slip occurs in the direction in which Ks·D>τp and value of τp 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 componentsDki+1 (k=1,2,3).(iii)
Check the tolerance ofDki+1-Dki for all the elements.
(1)
IfDki+1-Dki<tolerance 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)
IfDki+1-Dki>toleranceDki+1-Dki>tolerance. for any element, then solution is not convergent. Replace Dki (k=1,2,3) with Dki+1 (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)
IfDki+1-Dki>tolerance for any element, then solution is not convergent. Formula (16) would be used to calculate G, the energy release rate of the crack front element. If all the calculated Gs are less than or equal to Gc, (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 Figure5 [35]. The calculation results show that the software can simulate engineering problems correctly.Figure 5
Trajectory of joint propagation under compression.
(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 Figure6 in step0. The mechanical parameters of rock and joint are shown in Table 1.Table 1
Parameters of the joint propagation model.
Rock Joint  Elastic modulus (GPa) Poisson’s ratio Shear coefficient (m) Shear coefficient (n) c (MPa) φ (°) Maximum closure (m) Initial normal stiffness (MPa)  4 0.2 0.01 0.2 0.05 38.6 1 15Figure 6
3D nonlinear joint propagation process under compression.


### 7.1. Propagation Process under Compression

The three-dimensional propagation under vertical compression 100MPa is shown in Figure6, 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 Figure7. The process proceeds for 11 steps, and step0 is the initial joint shape.Figure 7
3D nonlinear joint propagation process under tension.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 Cohesionc

Joint cohesion c is usually small. When keeping other parameters’ values the same as those in Table1, 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.Figure 8
Comparison of joint propagation with different cohesion under compression.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 Table1, 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.Figure 9
Comparison of joint propagation with different friction under compression.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 Table1, 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.Figure 10
Comparison of joint propagation under compression with different lateral pressure coefficient in axial-Y.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.Figure 11
Comparison of joint propagation with different lateral pressure coefficient in axial-X under compression.


## 7.1. Propagation Process under Compression

The three-dimensional propagation under vertical compression 100MPa is shown in Figure6, 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 Figure7. The process proceeds for 11 steps, and step0 is the initial joint shape.Figure 7
3D nonlinear joint propagation process under tension.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 Cohesionc

Joint cohesion c is usually small. When keeping other parameters’ values the same as those in Table1, 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.Figure 8
Comparison of joint propagation with different cohesion under compression.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 Table1, 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.Figure 9
Comparison of joint propagation with different friction under compression.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 Table1, 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.Figure 10
Comparison of joint propagation under compression with different lateral pressure coefficient in axial-Y.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.Figure 11
Comparison of joint propagation with different lateral pressure coefficient in axial-X under compression.


## 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.



---

*Source: 3164817-2019-07-16.xml*