# Application of D-CRDM Method in Columnar Jointed Basalts Failure Analysis

Authors: Changyu Jin; Xiating Feng; Chengxiang Yang; Dan Fang; Jiangpo Liu; Shuai Xu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Columnar jointed basalt is a type of joint rock mass formed by the combined cutting effect of original joints and aphanitic microcracks. After excavation unloading, such rock mass manifested distinct mechanical properties including discontinuity, anisotropy, and proneness of cracking. On the basis of former research findings, this paper establishes a D-CRDM method applicable to the analysis of columnar jointed basalt, which not only integrates discrete element and equivalent finite-element methods, but also takes into account the coupling effect of original joints and aphanitic microcracks. From the comparative study of field monitoring data and strain softening constitutive model calculated results, it can be found that this method may well be used for the simulation of mechanical properties of columnar jointed basalts and the determination of rock failure mechanism and failure modes, thus providing references for the selection of supporting measures for this type of rock mass.


---

## Body

## 1. Introduction

Due to the abundance of internal joints, columnar jointed basalt possesses distinct mechanical properties such as discontinuity, anisotropy and is easy to crack compared with other types of rock mass. For this type of rock mass, its surface rocks instantaneously generate tensile failure right with the beginning of excavation; then aphanitic microcracks start to develop throughout the columns along with the reduction of surface rock confining pressure and induce a secondary failure of original joints, which finally leads to a progressive failure pattern similar to “domino effect.” Currently, studies on the unloading mechanical properties of columnar jointed rock mass appear to be quite limited. Based on existing research findings on joint rock mass [1–6], this paper proposes a D-CRDM (discrete element method-crack rock mass deterioration model) analysis method by integrating equivalent continuum and discrete element and uses this method to simulate the unloading mechanical impacts of rock mass after excavation of experimental cavity. The correctness of this analytical method is verified by comparing with the monitoring results; at the same time, the unloading failure mechanism and failure modes of columnar jointed basalts are revealed as well.


## 2. Failure Modes of Columnar Jointing

Columnar jointed basalt is a type of joint rock mass formed by the combined cutting effect of original joints and aphanitic microcracks, of which the original joint is a type of geological structure generated during the magmatic condensation process. It possesses a clear column surface and a distinct column outline; its cross-section is quadrilateral or pentagonal in shape and has a mosaic structure. Due to the weak strength of original joints, original joints on the surface are prone to spall after excavation. From photos taken on the damaged sites (see Figures1(a) and 1(b)), it can be observed that such failure and stripping are mainly caused by the opening and sliding of joints, which may be described by block theories and can be analyzed by discontinuous methods such as the discrete element method.Figure 1
Failure modes of surrounding rock in the cavity.The aphanitic microcracks inside the columns are small joint planes between the two crystal grains. Such small joint planes reduce the rock integrity. Through on-site observations it can be found that large scale spall appeared in the cavern roof or arch corners and created a not very deep pit. According to the sizes of the stripping blocks, it can be initially estimated that this type of failure was induced by the extending and splitting of aphanitic microcracks. In addition, the surrounding rock on the cavern side wall generated a large amount of broken rocks, especially in the region where the cracking failures of original joints are concentrated, which is seen in the failure model shown in Figures1(c) and 1(d). Since the aphanitic microcracks are randomly developed inside the columns and it is impossible to count their number, it is more appropriate to adopt equivalent method for analysis and research in light of the research results of the literature [7–11].


## 3. D-CRDM (Discrete Element Method-Crack Rock Mass Deterioration Model)

D-CRDM analysis method simulates original joints by 3DEC block theory and develops constitutive deterioration model of jointed rock by using its secondary development function in order to implement coupling analysis of these two kinds of joints.After excavation, the aphanitic microcracks on the surface of surrounding rocks gradually extend and crack under the effect of unloading, resulting in the deterioration of overall mechanical properties of columnar jointed basalt. From a macro perspective, this type of deterioration can be considered as a changing process of Young’s modulusE, Poisson's ratio μ, cohesion C, friction angle ϕ, tensile strength T, and other mechanical parameters of the rock mass. Therefore, when adopting numerical simulation method to study the mechanical behaviour of columnar jointed basalt, its constitutive model should accurately reflect two fundamental features of the rock mass; that is, the mechanics parameters change along with the damage of surrounding rock both during and after the yielding course of rock mass. In rock deterioration process, each parameter change should be determined by the mechanical state of rock mass, while the mechanical index of equivalent plastic strain ε-p is usually used to describe different stress state of rock mass, which is defined
(1)ε-p=∫23(ε1pε1p+ε2pε2p+ε3pε3p).
In order to define the change of deformation and strength parameters, it is supposed that each parameter satisfies the function of ε-p the following:
(2)E(ε-p)=E0·fE(ε-p),μ(ε-p)=μ0·fμ(ε-p),C(ε-p)=C0·fC(ε-p),ϕ(ε-p)=ϕ0·fϕ(ε-p),T(ε-p)=T0·fT(ε-p),
where E0, μ0, C0, φ0, and T0 are Young’s modulus, Poisson’s ratio, cohesion, friction angle, and tensile strength of the columns of columnar jointed basalt in initial state; E(ε-p), μ(ε-p), C(ε-p), φ(ε-p), and T(ε-p) are Young’s modulus, Poisson’s ratio, cohesion, friction angle and tensile strength for a given plastic strain state, respectively; whereas fE(ε-p), fμ(ε-p), fc(ε-p), fφ(ε-p), and fT(ε-p) are the evolution functions which may be linear, piecewise, or nonlinear and can be determined in the aid of the experiences or the correlated experimental test curves.After excavation of underground cavity, the surrounding rock mass suffered different degrees of damage due to the cracking of aphanitic microcracks in surface columns and new fractures generated hence with. Theoretical analysis and engineering practice both show that the cracks and expansion of joints reduce the rock integrity, resulting in reduced resistance to deformation of rock mass, and its elastic modulusE also decreases with it. When the rock mass is within the scope of elastic range, the Poisson’s ratio is generally constant; however, Poisson’s ratio will change along the rock mass stress state alteration when the rock mass is beyond the scope of elastic range. Generally, after dilatation of rock mass, the Poisson's ratio of rock mass often increases until it reaches 0.5. Therefore, in rock degradation process, the computation of deformation parameters is a progressively variant procedure along with the change of the mechanical state of rock mass. Since the breaking of original joints continuously reduces columns confining pressure, the tiny cracks inside the columns are progressively cut through, leading to a macroscopic manifestation of large deformation as well as the stripping and falling of aphanitic microcracks. These local deterioration phenomena of rock mass reduce its strength and integrity, which corresponds to the reduction of mechanical parameters including friction angle, cohesion, and tensile strength. Therefore, the process of aphanitic microcracks breakage and new fracture development may be considered as a continuous disintegrating course of the microstructures of the columns, that is, the process of reductions of friction angle φ, cohesion C, and tensile strength0020T along with the reduction of equivalent plastic strains. In this paper, an example with linear evolution functions is presented in Figure 2.Rock mass mechanical parameters changing curve diagram along with the variation of equivalent plastic strain.
(a)
(b)
(c)
(d)
(e)To achieve the dynamical evolution of above mechanical parameters along with the variation of plastic strain, combined with Mohr-Coulomb elastoplastic constitutive model, the mechanical parameters of the materials are step-by-step updated according to (2) in iterations calculated in the 3DEC difference dynamical iterative solution procedure (see Figure 3).Figure 3
Procedure of the numerical computation used in CRDM.In the process of circulative calculation, the shear yield function uses Mohr-Coulomb yield criterion, which is expressed as a function of equivalent plastic strain; namely,(3)fs=σ1-σ3Nφ+2C(ε-p)Nφ,Nφ(ε-p)=1+sin[φ(ε-p)]1-sin[φ(ε-p)].
Taking into account the weakness of tensile strength of columnar jointed basalt, Rankine’s strength criteria of maximum tensile stress is also considered when judging whether the rock mass enters into shaping state:
(4)ft=σt-σ3.
After the cracking of aphanitic microcracks, the rock mass goes into yield condition, and since each parameter constantly changes along with the equivalent plastic strain, the yielding surface of this model varies dynamically not as that given in the common elastoplastic constitutive relations. According to the research results of Yingren et al. [12], the straight line equation of side AB in the π plane and the equation of distance between the centre of the π plane and the origin of the principle stress space coordinates system are given as
(5)x=2ccosφ+sinφ3y-σmsinφ,d=13(σ1+σ2+σ3)=3σm,
where σm is the hydrostatic pressure.The above equations indicate that(1) the slope of line AB is only related to the friction angle φ; the slope of line AB decreases when φ is increasing and deviates to line AB′, which means some changes have taken place in the shape of the yielding surface (see Figure 4(a)); (2) when the cohesion changes, the intercept of line AB on x-axis changes correspondently, and line A′B′ may become a parallel translation of the straight line AB, which means that the area of the yielding surface in the π plane will expand in a similarity form as shown in Figure 4(b). As the distance between the π plane and the origin of the principle stress space coordinates system depends only on the hydrostatic pressure σm, when the internal stress of rock mass changes, that is, when σm changes, the yielding surface will be scaled in the meridian plane along the isoclinic line direction (see Figure 4(c)). In addition to the dynamic changes of tensile strength along with the evolution of equivalent plastic strain, this constitutive model presents sufficiently well the dynamic changes of the mechanical properties of columnar jointed basalts after their yielding.Evolution characteristics of the successive yielding in CRDM.
(a)
Evolution of the yielding surface in theπ surface
(b)
Isoclinic scaling of the yielding surface in theπ surface
(c)
Isoclinic scaling of the yielding surface


## 4. Experimental Cavity of Columnar Joints

To verify the correctness of this analysis method and to analyze the unloading mechanism of columnar joints, simulation and analysis were carried out in the experimental cavity of columnar joints located in the Baihetan hydropower station, which is situated in Ningnan county of Sichuan province and Qiaojia county of Yunnan province, the downstream of the Jinsha River. Its dam foundation and the main part of cavity groups are located in the basalts of Upper Permian Emeishan formation(P2β3), where the Emeishan basalts are divided into eleven rock-flowage layers, of which the P2β3 is a basaltic formation mainly consisted of microcrystalline aphanitic joints and oblique lamination. It is the major stratum outcropped in the project area, and columnar joints are developed in partial areas of the middle part of this stratum (P2β32 and P2β33). In order to explore the unloading properties of columnar jointed basalts and to ensure the stability of dam foundation and underground cavities on the long run, an experimental cavity was excavated in the P2β32 stratum for long-term monitoring. This experimental cavity is located in the downstream side of cave PD34 at right bank prospecting line II, with an axial direction of N45°E that is almost consistent with the strike direction of stratum. The cavity has a floor elevation of Δ727.4 m, a vertical embedded depth of about 300 m, and is approximately 160 m away from the bank slope. The in situ stress values in this region are relatively small: the first and the third principal stresses are 4.38 MPa, and 3.23 MPa respectively, with little difference between the two. With reference to the data obtained from laboratory experiments and field triaxial load/unload tests, the calculation parameters of rock mass and structure surfaces were determined as in Tables 1 and 2.Table 1
Mechanical parameters of joints.
Normal stiffness/GPa
Shear stiffness/GPa
Frictional angle/°
Cohesion/MPa
Tensile strength/MPa
52.97
12.88
24.3
0.0
0.0Table 2
Equivalent mechanical parameters of the rock.
Young's modulus/GPa
Poisson’s ratio
Friction angle/°
Cohesion/MPa
Tensile strength/MPa
Initial
Residual
Initial
Residual
Initial
Residual
Initial
Residual
Initial
Residual
32.30
15.20
0.21
0.50
35.0
25.0
1.30
0.50
1.25
0.0The experimental cavity of columnar jointed basalts, with a total length of 70 m, is divided into two experimental segments, that is, the supported and unsupported experimental segments of 30 m and 40 m, respectively. The cross-section of this experimental cavity is city-gate shaped with dimensions of 13 m × 6.5 m (height × width). According to the geological survey for monitoring section in the unsupported segment, a two-dimensional computational model for the experimental cavity is constructed in which three rock flow age strataP2β32, P2β33 and P2β34 from the bottom to the top are considered. What is more, there is an intraformational disturbed belt of RS3311 in stratum P2β32. In order to eliminate the influence of boundary constraints on calculation and to accelerate computing speed, the rock masses within the range of three times the span of cavity and twice the height of cavity were cut into columnar joints, while the rock masses outside this scope were modeled as continuum media. The columnar jointed basalts in rock flowage stratum P2β33, where the experimental cavity is located, are ash black ones whose columns are 1~5 m long, 15~25 cm in diameter, and 70°~85° in dip angle; see Figure 5 for specific model. This model adopts tetrahedral elements, with an internal column elements size of 0.15 m, an integrity rock unit of 0.5 m, and total concrete units of 15,000.Figure 5
Calculated model of computation.


### 4.1. Analysis of the Displacement Field of Surrounding Rock

The excavation process of the experimental cavity is analyzed by utilizing the simulation method established in this paper. The results reveal that during the excavation of the first layer, deformation of surrounding rock at surface was relatively small; its numerical simulation result is 2.35 mm. During the excavation of the second layer, deformation at each monitoring point increased rapidly (as shown in Figure6); the numerical simulation value and monitoring datum reached 16.68 mm and 17.80 mm, respectively. These numerical simulation results reveal that the principle cause for the large deformation of surrounding rock is the existence of normal and shear deformations of original joints in rock mass, which are mainly concentrated within 0~4 m of the surrounding rock surface. In addition, the deformation of surrounding rock is also influenced by the dip angles of columnar joints and manifests distinct anisotropic properties; that is, the deformation value of the cavity roof is larger than that of the cavity floor, and the deformation value at the downstream side of cavity is larger than the displacement value of cavity sidewall at the upstream side (as shown in Figure 7). From the comparative study of in situ monitoring results and calculated results of strain-softening constitutive model, it can be found that D-CRDM may accurately simulate the displacement field alteration of surrounding rock during excavation.Figure 6
Surrounding rock deformation and rock mass failure.Numerical simulation results of the surrounding rock deformations.
(a)
Normal deformation of original joints
(b)
Tangential deformation of joints
(c)
Deformation of surrounding rocks


### 4.2. Relaxation of Joints

The failure of jointed rock mass is often due to the expansion and cracking of joints, while the development of these fissures will lead to the reduction of wave velocity. According to the reduction degree of wave velocity and by modeling on the rock mass integrity evaluation system, a relaxation coefficient of rock massk is established:
(6)k=vp0-vpvp0,
where k is the coefficient of relaxation, vp is the velocity of longitudinal wave of rock mass at the testing position, and vp0 is the eigen wave velocity of nonrelaxed rock mass.In order to distinguish the relaxation degree of rock mass in a more detailed manner, the relaxation degree of columnar joint is subdivided into highly, weakly, and slightly relaxed ones (see Table3 for further details).Table 3
Classification criteria of relaxation degrees.
k
k
≤
0
0
<
k
≤
20
%
20
%
<
k
≤
40
%
k
>
40
%
v
p/(m/s)
v
p
≥
5000
5000
>
v
p
≥
4000
4000
>
v
p
≥
3000
v
p
<
3000
Relaxation evaluations
Nonrelaxed
Slightly relaxed
Weakly relaxed
Strongly relaxedThe calculated results obtained by D-CRDM method indicate that after excavation of the first layer of the experimental cavity, the open original joints of upstream and downstream cavity sidewalls were concentrated within 0.5 m of the cavity surface layer. Due to the effect of dip angles of original joints, the joints development area at arch roof basically diffuses along the direction of joint dip angle, with an opening depth almost twice the height of cavity excavation, that is, around 10 m. This is generally consistent with the distribution regularity of relaxation zones of rock mass obtained by in situ acoustic wave tests (Figure8).Relaxation zones of joints after excavation of the first layer.
(a)
Original joints
(b)
Aphanitic joints
(c)
Relaxation zone of columnar jointsWhen the experimental cavity was cut into to the second layer, from the contour diagram of joint slipping (as shown in Figure9) it can be found clearly that the opening zones of original joints on the cavity roof and floor were significantly enlarged, and their heights were also almost equivalent to the cavity height, that is, about 10 m. In addition, a cracking zone of certain depth also appeared on the cavity sidewall. By combining the normal and tangential displacement contour diagrams of original joints, it reveals that the failure and relaxation degree of original joints on cavity sidewall is more intense than those on the roof and floor, with a failure depth of 2~3 m. Due to the same effect of original joints, the joint relaxation zones still develop along the joint dip angle direction, which fully complies with actual observations.Relaxation zones of joints after excavation of the second layer.
(a)
Original joints
(b)
Aphanitic joints
(c)
Relaxation zone of columnar jointsIn order to verify the analysis results obtained from wave velocity tests and numerical simulations, a three-dimensional scanning was carried out by using a borehole camera technique on the surrounding rocks inside the borehole that was drilled along the downstream sidewall of experimental cavity. The results showed that the rock mass within 0.0~4.0 m of the surrounding rock on sidewall was obviously relaxed, which further prove the accuracy of above theoretical analysis and field testing (see Figure10).Figure 10
Borehole camera information of columnar joints.


### 4.3. Field Bearing Plate Test

In order to further explore the compression features of columnar joints, a bearing plate test was carried out in the areas where the columnar joints are densely distributed and where the aphanitic microcracks are developed. After the circulative loading and unloading tests of multiple experimental points, it is discovered that most stress-displacement curves are of concave type; that is, under the effect of circulative loading and unloading step-by-step, the Young's modulus of rock mass gradually increases; refer to Figure11(a) for specific stress routes. Using the formerly proposed method of D-CRDM, a simulation analysis was made to the experiment. As the original joints gradually close during the step-by-step pressuring process of bearing plate, together with the decline of equivalent plastic strain of columns and the gradual increase of Young's modulus, the numerical simulation results basically match with the experimental results, thus further verify that this model complies with the reloading mechanical properties of columnar joints after unloading.Contrast between field test and calculated results.
(a)
Typical testing curves of field bearing plate
(b)
Simulated results of field bearing plate


## 4.1. Analysis of the Displacement Field of Surrounding Rock

The excavation process of the experimental cavity is analyzed by utilizing the simulation method established in this paper. The results reveal that during the excavation of the first layer, deformation of surrounding rock at surface was relatively small; its numerical simulation result is 2.35 mm. During the excavation of the second layer, deformation at each monitoring point increased rapidly (as shown in Figure6); the numerical simulation value and monitoring datum reached 16.68 mm and 17.80 mm, respectively. These numerical simulation results reveal that the principle cause for the large deformation of surrounding rock is the existence of normal and shear deformations of original joints in rock mass, which are mainly concentrated within 0~4 m of the surrounding rock surface. In addition, the deformation of surrounding rock is also influenced by the dip angles of columnar joints and manifests distinct anisotropic properties; that is, the deformation value of the cavity roof is larger than that of the cavity floor, and the deformation value at the downstream side of cavity is larger than the displacement value of cavity sidewall at the upstream side (as shown in Figure 7). From the comparative study of in situ monitoring results and calculated results of strain-softening constitutive model, it can be found that D-CRDM may accurately simulate the displacement field alteration of surrounding rock during excavation.Figure 6
Surrounding rock deformation and rock mass failure.Numerical simulation results of the surrounding rock deformations.
(a)
Normal deformation of original joints
(b)
Tangential deformation of joints
(c)
Deformation of surrounding rocks


## 4.2. Relaxation of Joints

The failure of jointed rock mass is often due to the expansion and cracking of joints, while the development of these fissures will lead to the reduction of wave velocity. According to the reduction degree of wave velocity and by modeling on the rock mass integrity evaluation system, a relaxation coefficient of rock massk is established:
(6)k=vp0-vpvp0,
where k is the coefficient of relaxation, vp is the velocity of longitudinal wave of rock mass at the testing position, and vp0 is the eigen wave velocity of nonrelaxed rock mass.In order to distinguish the relaxation degree of rock mass in a more detailed manner, the relaxation degree of columnar joint is subdivided into highly, weakly, and slightly relaxed ones (see Table3 for further details).Table 3
Classification criteria of relaxation degrees.
k
k
≤
0
0
<
k
≤
20
%
20
%
<
k
≤
40
%
k
>
40
%
v
p/(m/s)
v
p
≥
5000
5000
>
v
p
≥
4000
4000
>
v
p
≥
3000
v
p
<
3000
Relaxation evaluations
Nonrelaxed
Slightly relaxed
Weakly relaxed
Strongly relaxedThe calculated results obtained by D-CRDM method indicate that after excavation of the first layer of the experimental cavity, the open original joints of upstream and downstream cavity sidewalls were concentrated within 0.5 m of the cavity surface layer. Due to the effect of dip angles of original joints, the joints development area at arch roof basically diffuses along the direction of joint dip angle, with an opening depth almost twice the height of cavity excavation, that is, around 10 m. This is generally consistent with the distribution regularity of relaxation zones of rock mass obtained by in situ acoustic wave tests (Figure8).Relaxation zones of joints after excavation of the first layer.
(a)
Original joints
(b)
Aphanitic joints
(c)
Relaxation zone of columnar jointsWhen the experimental cavity was cut into to the second layer, from the contour diagram of joint slipping (as shown in Figure9) it can be found clearly that the opening zones of original joints on the cavity roof and floor were significantly enlarged, and their heights were also almost equivalent to the cavity height, that is, about 10 m. In addition, a cracking zone of certain depth also appeared on the cavity sidewall. By combining the normal and tangential displacement contour diagrams of original joints, it reveals that the failure and relaxation degree of original joints on cavity sidewall is more intense than those on the roof and floor, with a failure depth of 2~3 m. Due to the same effect of original joints, the joint relaxation zones still develop along the joint dip angle direction, which fully complies with actual observations.Relaxation zones of joints after excavation of the second layer.
(a)
Original joints
(b)
Aphanitic joints
(c)
Relaxation zone of columnar jointsIn order to verify the analysis results obtained from wave velocity tests and numerical simulations, a three-dimensional scanning was carried out by using a borehole camera technique on the surrounding rocks inside the borehole that was drilled along the downstream sidewall of experimental cavity. The results showed that the rock mass within 0.0~4.0 m of the surrounding rock on sidewall was obviously relaxed, which further prove the accuracy of above theoretical analysis and field testing (see Figure10).Figure 10
Borehole camera information of columnar joints.


## 4.3. Field Bearing Plate Test

In order to further explore the compression features of columnar joints, a bearing plate test was carried out in the areas where the columnar joints are densely distributed and where the aphanitic microcracks are developed. After the circulative loading and unloading tests of multiple experimental points, it is discovered that most stress-displacement curves are of concave type; that is, under the effect of circulative loading and unloading step-by-step, the Young's modulus of rock mass gradually increases; refer to Figure11(a) for specific stress routes. Using the formerly proposed method of D-CRDM, a simulation analysis was made to the experiment. As the original joints gradually close during the step-by-step pressuring process of bearing plate, together with the decline of equivalent plastic strain of columns and the gradual increase of Young's modulus, the numerical simulation results basically match with the experimental results, thus further verify that this model complies with the reloading mechanical properties of columnar joints after unloading.Contrast between field test and calculated results.
(a)
Typical testing curves of field bearing plate
(b)
Simulated results of field bearing plate


## 5. Conclusions

By utilizing the D-CRDM analysis method established in this paper, mechanical simulations of both the loading and the reloading after unloading mechanical properties of columnar joints maybe carried out, and multiple field test results have proved the correctness of this method. Meanwhile, the research results of a variety of means also indicate that the failure and instability of columnar jointed basalts are mainly caused by the open failure of original joints; therefore, reinforcing the support of original joints is fundamental to the stability insurance of columnar jointed basalts.



---

*Source: 848324-2013-09-24.xml*