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

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, andproneness of cracking.On the basis of former research findings, this paper establishes aD-CRDMmethod applicable to the analysis of columnar jointed basalt, which not only integrates discrete element and equivalent finite-elementmethods, 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 methodmay 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.


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][2][3][4][5][6], this paper proposes a D-CRDM (discrete element methodcrack 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.

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 Figures 1(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.
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 Figures 1(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][8][9][10][11].

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 modulus , Poisson's ratio , cohesion , friction angle , tensile strength , 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   is usually used to describe different stress state of rock mass, which is defined  In order to define the change of deformation and strength parameters, it is supposed that each parameter satisfies the function of   the following: where  0 ,  0 ,  0 ,  0 , and  0 are Young's modulus, Poisson's ratio, cohesion, friction angle, and tensile strength of the columns of columnar jointed basalt in initial state; (  ), (  ), (  ), (  ), and (  ) are Young's modulus, Poisson's ratio, cohesion, friction angle and tensile strength for a given plastic strain state, respectively; whereas   (  ),   (  ),   (  ),   (  ), and   (  ) 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 modulus  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 , and tensile strength T along with the reduction of equivalent plastic strains.In this paper, an example with linear evolution functions is presented in Figure 2.
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-bystep updated according to (2) in iterations calculated in  the 3DEC difference dynamical iterative solution procedure (see Figure 3).
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, 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: 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  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 where   is the hydrostatic pressure.
The above equations indicate that (1) the slope of line  is only related to the friction angle ; the slope of line  decreases when  is increasing and deviates to line   , 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  on axis changes correspondently, and line     may become a parallel translation of the straight line , 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   , when the internal stress of rock mass changes, that is, when   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.

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 ( 2  3 ), where the Emeishan basalts are divided into eleven rock-flowage layers, of which the  2  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 ( 2  2 3 and  2  3 3 ).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  2  2 3 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.
The 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  .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  2  3  3 , 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.

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 Figure 6); 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.

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 mass  is established: where k is the coefficient of relaxation, V  is the velocity of longitudinal wave of rock mass at the testing position, and V 0  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 Table 3 for further details).
The 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 (Figure 8).
When the experimental cavity was cut into to the second layer, from the contour diagram of joint slipping (as shown in Figure 9) 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.In 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.0m of the surrounding rock on sidewall was obviously relaxed, which further prove the accuracy of above theoretical analysis and field testing (see Figure 10).

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 Figure 11(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.

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.

Figure 1 :
Figure 1: Failure modes of surrounding rock in the cavity.

Figure 2 :
Figure 2: Rock mass mechanical parameters changing curve diagram along with the variation of equivalent plastic strain.
Computation of elastic trial stressExamination of plasticCorrection of the stress and Updating of the material parameters of the elements yielding and computation of the plastic multiplier computation of the plastic strain

Figure 3 :
Figure 3: Procedure of the numerical computation used in CRDM.

Figure 4 :
Figure 4: Evolution characteristics of the successive yielding in CRDM.

Figure 6 :
Figure 6: Surrounding rock deformation and rock mass failure.

Figure 7 :
Figure 7: Numerical simulation results of the surrounding rock deformations.

Figure 8 :
Figure 8: Relaxation zones of joints after excavation of the first layer.

Figure 9 :
Figure 9: Relaxation zones of joints after excavation of the second layer.
Simulated results of field bearing plate

Figure 11 :
Figure 11: Contrast between field test and calculated results.

Table 1 :
Mechanical parameters of joints.

Table 2 :
Equivalent mechanical parameters of the rock.3 3 and  2 43 from the bottom to the top are considered.What is more, there is an intraformational disturbed belt of RS3311 in stratum  2  2 3

Table 3 :
Classification criteria of relaxation degrees.