Based on several concepts of rock mass mechanics，this paper comprehensively reviewed and analyzed the important achievements and deficiencies in the theory and application of rock mass mechanics all around the world. In the part of rock mechanics testing machine and testing method，it mainly introduced triaxial testing machine，rigid testing machine，real triaxial testing machine，rheology testing machine，impact testing machine，high-temperature and high-pressure testing machine，multi-field coupling action testing machine，CT-Rock testing machine，in-situ test and test standards，etc. In the part of constitutive law，the whole stress-strain curve，the mechanical properties of rock under pseudo triaxial and real triaxial stresses，aging and size effect characteristics，dynamic characteristics，seepage characteristics，multi-field coupling characteristics，mechanical characteristics of structural surfaces，optoacoustic/electromagnetic/themal effects during rock deformation and failure were mainly introduced. In the theory of rock mass mechanics，the mechanical categorization of rock mass media，rock mass mechanics of matrix-fractured media，strength criteria，constitutive laws，fracture and damage mechanics，multi-field coupling model and fracture distribution model were introduced in detail. Numerical methods and software，displacement back analysis and intelligent analysis methods were introduced in the numerical calculation section. The classification，concepts and application fields of engineering rock mass mechanics and disaster rock mass mechanics were expounded，the rock mass mechanics principles of major projects such as dam engineering，tunnel engineering，mining engineering，petroleum and unconventional resource development engineering were analyzed and condensed，and the important achievements of various historical stages in the evolution of engineering rock mechanics were listed. The rock mass mechanics principles in the initiation and exacerbation of natural and engineering disasters，such as landslides，gas outbursts，rock bursts，earthquakes were combed，and the important achievements of disaster rock mass mechanics in disaster prediction and prevention technology were enumerated. Through in-depth analysis，eight unsolved problems of rock mass mechanics were put forward in the past century，including three theoretical problems of rock mass mechanics，such as media classification theory，control effect of defect levels on the deformation and failure of rock masses，and anisotropic rock mass mechanics theory and analysis methods，etc. The other five questions are related to nonlinear rock mass mechanics，namely，rock mass scale effect，time effect，unified theory of catastrophe-chaos-percolation in the process of instability and failure of rock mass system，and predicting the failure time，failure location and released energy in the instability and failure process of the intact rock sample and the rock mass system.

The size effect of rock joint shear strength is a basic research issue in the field of rock mechanics and engineering geology. Due to the lack of mature direct shear test equipment on series size specimens at home and abroad，the research process on shear strength size effect of rock joints is slow. In this paper，a shear strength size effect direct shear combined test system of rock joints is developed，and the key technologies and methods are put forward such as fast conversion of multi-scale jointed specimens，classification of load measurement and control system，equivalence check of overlapping specimens，dynamic detection in testing process and shear loss measurement of structural planes，providing a basis experimental support platform for investigating the size effect of rock joint shear strength and other key scientific issues. The preliminary application results of the developed test system are demonstrated by two-level direct shear apparatus. It is shown that the test system can solve the problem of automatic transformation of specimens in the range of 100 times of the rock joint size from 10 cm×10 cm to 100 cm×100 cm，involving the process of transporting，positioning，loading and resetting of multi-scale specimens. The test system can also solve the problem of accuracy and stability of the measurement control system in the range of 500 times of the load from 2 kN to 1 000 kN. The accuracy level of the test system reaches 0.5 through the identification of a testing organization. The shear strength test results of overlapped specimens of 30 cm×30 cm and 40 cm×40 cm in two-level direct shear apparatus are compared，which shows that the convergence of the two sizes of samples is good. For the specimens with the same structural plane sizes under different loads，the peak values of the tangential load are related to the applied normal loads. For the specimens with different structural plane sizes under the same load，however，the ratio of the tangential load peak value to the normal load has a stable trend，except for a slight fluctuation in a certain size range of the structural plane.

Cardox tube system(CO2 gas cracking technology) is conductive to reducing the side effect of vibrations in rock excavation. To quantify this reducing effect，field tests of gas cracking and explosive blasting were conducted at the foundation pit of a nuclear power plant. At the same time, the ground vibrations induced by the two methods were monitored respectively. On the basis of energy equivalence，the phase transition energy of liquid CO2 was calculated and converted to emulsion explosive mass with the same amount of energy. By introducing a PPV(peak particle velocity) ratio of the two methods，the attenuation characteristics of ground vibrations were compared and the considerable discrepancy between the PPVs was revealed，so as to provide guidance for the optimization of rock excavations and technical configurations of Cardox tube system. It is observed that, under the same equivalent charge amount，the vibrations induced by CO2 gas fracturing are much less than those by explosives，and attenuate more rapidly. Generally，the latter are at least 5 times more than the former，and the difference keeps growing with increasing the distance and decreasing the equivalent charge amount，which means that Cardox tube system is more efficient in vibration reduction in middle to far areas from the rock excavation point，and especially efficient in local and small-scale rock excavating. It is also indicated that the explosive consumption per unit volume of rock is about 5 times the equivalent value of Cardox tube system，as is the reason for the remarkable discrepancy in PPVs by the two methods.

In view of the lack of understanding of the toppling deformation of anti-dip soft-hard interbedded rock slopes，based on similarity principle and geological data，three groups of slope physical models were established with cement and gypsum as similar materials，among which one group is an anti-dip layered rock slope with hard rock and the other two groups are soft-hard interbedded rock slopes with different layer thickness ratios. By carrying out centrifugal model tests and using image measurement technology，the differences of toppling deformation and failure modes between the anti-dip soft-hard interbedded rock slopes and the anti-dip layered rock slopes are studied，and the influence of the layer thickness ratios of soft rock to hard rock on the overall toppling deformation degree and the ultimate bearing capacity of the slope is analyzed. The following conclusions are obtained by the tests：(1) The pattern of toppling deformation of the anti-dip soft-hard interbedded rock slope is different from the anti-dip layered rock slope with hard rock. The former has two fracture surfaces including the primary fracture surface and the secondary fracture surface. The secondary fracture surface is first formed and the upper rock mass has an instability. Then，the deep discontinuous bending zones coalesce with each other，and finally the main fracture surface is formed and the slope fails as a whole collapsing downward. (2) Almost no toppling occurs in the rock mass below the primary fracture surface during toppling deformation，so the primary fracture surface can be defined as the boundary line between the toppling and non-toppling rock masses. The development depth of the secondary fracture surface is shallower than that of the primary fracture surface，but the shallow rock mass above the secondary fracture surface has a greater toppling deformation degree and is more prone to failure. The secondary fracture surface is the most dangerous fracture surface during the process of slope toppling deformation. (3) Due to the weak strength of the soft rock，the fracture surface of the anti-dip soft-hard interbedded rock slope is arc-shaped，which is different from that of single lithology layered anti-dip rock slope. (4) The existence of the soft rock also has influence on the ultimate bearing capacity and the toppling deformation degree of the slope，varying with the layer thickness ratio of the soft rock to the hard rock. Compared with the single hard rock layered anti-dip rock slope，the ultimate bearing capacity of the soft-hard interbedded rock slope with a thickness ratio of 1∶1 increases and the toppling deformation degree decreases，while for the soft-hard interbedded rock slope with a thickness ratio of 2∶1，the ultimate bearing capacity decreases and the toppling deformation degree increases.

In order to study the change of the wave velocity and the law of energy evolution of red sandstone under triaxial loading，compression and shear wave velocity tests of red sandstone in the full stress-strain process under different confining pressures were carried out by using the Rock Top multi-field coupling tester and the JSR-DP300 ultrasonic testing system. The results show that the compression and shear wave velocities increase as a power function of the confining pressure，accompanied by increasingly nonlinear characteristics，and that the wave velocity data points are distributed in a long and narrow interval. The confining pressure only has a gain effect on the macroscopic strength of rock，and the microscopic difference is still very obvious. In the process of triaxial compression，the compressional wave velocity can be divided into three stages including continuous growth，oscillation and rapid attenuation. Under the condition of high confining pressures，the phase characteristics of the shear wave velocity are not obvious，and the curve of the shear wave velocity is nearly linear and relatively smooth. The greater the confining pressure is，the weaker the gain of the wave velocity is in the axial loading process. The attenuations of the compression and shear wave velocities are synchronous，and the compression wave velocity is more sensitive to the state of stress difference. Under different confining pressures，the overall shapes of the energy density curves are consistent，and the attenuation interval of the wave velocity corresponds to the main growth interval of the dissipated energy. Before the rock damage stress，there is a good linear correlation between the compression and shear wave velocities，and the inflection point of the compression and shear wave velocities can be used as a mark to judge the damage stress. With increasing the axial strain，the dynamic elastic modulus under different confining pressures generally increases at first，then keeps relatively stable，and finally rapidly decreases. The dynamic Poisson's ratio shows an overall growth trend. After entering the IV stage，the curve fluctuates violently and the dynamic Poisson¢s ratio increases rapidly. The stage characteristics of the two curves are not obvious.

Precise localization of microseismic events plays an important role in rockburst monitoring and early warning，and although the positioning error of the sectional wave velocity model is substantially lower than that of the single wave velocity model，the positioning error is still large(28.51 m of an example in this paper). Therefore，this paper proposes a new microseismic localization method based on Snell's law and Cuckoo search algorithm to improve the localization accuracy in layered rock masses. The results show that the Cuckoo search algorithm can combine long and short step lengths to search for spatially optimal solutions using the Lévy flight search strategy，presents good robustness in the process of seismic source inversion by simulating the specialized breeding parasitic behavior of Cuckoos，and can achieve reliable localization of microseismic events. It is also indicated that the combination of the Cuckoo search algorithm and the Snell¢s law can effectively overcome the error caused by simplifying the elastic wave propagation path between two points to a straight line，and can significantly reduce the positioning error compared to the sectional wave velocity model(down to 0.12 m in this paper)

In order to study the rock strength size effect under dynamic impact，nine groups of granite samples with different length-diameter ratios(0.5，0.6，0.7，0.8，0.9，1.0，1.2，1.5 and 2.0) were prepared，and dynamic impact compression tests were carried out under different impact speeds using f 50 mm split Hopkinson bar (SPHB). The mechanical properties of granite specimens with different length-diameter ratios were studied，and the absorption energy of the broken specimens was analyzed. The results show that，under three impact velocity compression conditions，the maximum differences of the compressive strength of the samples with different length-diameter ratios are 10.88%，3.95% and 8.17%，respectively，and that，under the same impact velocity，the compressive strength does not change obviously with the change of the length-diameter ratio. The length-diameter ratio of the rock specimens is negatively correlated with the peak strain of the rock specimens，and the fitting data curve shows that there is a certain nonlinear relationship between them. The length-diameter ratio of the specimens is positively correlated with the elastic modulus，and the fitting data curve shows that there is a good linear relationship between them. The failure shapes of the rock specimens with different length-diameter ratios are different while the failure forms are basically the same. The magnitude of the incident energy is related to the impact velocity but independent of the length-diameter ratio of the rock specimens. The fitting curves between the reflected energy and the transmitted energy with the length-diameter ratio are contrary. The absorbed energy of rock specimen breakage does not change significantly with the length-diameter ratio，but the absorbed energy per unit length of rock specimen breakage decreases with increasing the length.

The roughness of joints in rock mass is very important for the calculation of the shear strength. However，the determination of JRC(joint roughness coefficient) based on the 10 standard curves established by Barton is often subjective and uncertain，so it is important theoretical significance and practical application values to do further study about this issue. In this paper，all the standard joint roughness curves by Barton are quantized. The length of the datum line L0 is selected as 10 cm and the contour length L of each curve is measured. The functional relationship between the section length ratio of the contour length to the datum line length L/L0 and JRC can be fitted. Based on the ellipse definition and the sine function method，the curves are redistributed according to their contour lengths，and the maximum undulation Ry1 and Ry2 under two extreme distribution states are obtained. Making the maximum undulations correspond to the upper and lower limits of JRC，then，the functional relationship between Ry and JRC can be fitted. Finally，by bivariate fitting of the above two functional relationships，the functional relationship between JRC with L/L0 and Ry is obtained with a fit goodness of 99.24%. Taking the discontinuities of the slopes in Heihe reservoir area and Zhen¢an County of Xi¢an city as the cases，comparison between the developed approach with the classical JRC calculation method is performed. It is found that the calculation results are similar while the change trend of the JRC value calculated by the proposed function is more stable than the classical method，verifying the feasibility of the developed method.

Reservoir resources often occur in anisotropic jointed rock mass and fracture mode of anisotropic reservoirs is relatively complex，which make it difficult to carry out efficient exploitation only based on engineering experiences. To better understand the fracture characteristics of anisotropic reservoirs，to realize the control and prediction of fracture expansion morphology，and hence，to collect reservoirs more scientifically，this paper uses the discrete element method(DEM) to study fracture characteristics of the anisotropic sandstone by semi-circular bending test under asymmetric loading. By adjusting the asymmetric degree of the bottom bearing， the control of rock mass fracture from pure I mode to pure II mode can be achieved. The main results are as follows：(1) The peak load of the semi-circular bending samples shows obvious anisotropy，showing distribution forms of “W” and “increasing-stable-increasing” with increasing the joint angle，and a trend of “precipitous drop-stable” with increasing the asymmetric coefficient. (2) SCB samples with different notch angles and joint distributions show different failure modes mainly including offset type，bedding type，cut-bedding type and composite bedding type. With increasing the asymmetric coefficient，the stress field in the samples changes，making fracture path of the samples deviate along the movement direction of the bearing. (3) The fracture toughness of the samples obtained by solving the geometric factor，which is obtained based on the calculation of J-integral and by which 3 kinds of asymmetric loading ways of pure II mode，pure I mode and compound mode are designed，shows obvious anisotropy with the change of the joint angle and presents a trend of “precipitous drop and stable” with increasing the asymmetric coefficient. (4) Impacted by the joints，the initial crack angles of pure I mode and pure II mode are respectively distributed in the range of 0°–36.9°and 47.68°–76.12°，which are not exactly the same as the theoretical values.

In order to research the supporting mechanical characteristics of advanced hydraulic support group and bolt (cable) coupling support for surrounding rocks，the Mindlin-Reissner plate theory was employed to describe the mechanical characteristics of the surrounding rock，and the energy functional was used to describe the supporting force of the equivalent support group and the bolt(cable) action. A mechanical model of the advanced hydraulic support group-anchor coupling support of surrounding rocks was established，in which the penalty function was introduced to constraint the model boundary strategies and the system control equations were obtained based on Hamilton principle. According to the advanced support parameters of the 15106 working face in Wenjiazhuang coal mine of Pingshu company of Shanxi Yangquan Group，the deformations as well as the bending moment and stress distributions of the roof under three strategies including the anchor support，the equal strength support and the non-equal strength support were obtained by solving the model with Ritz method. Comparisons of the supporting effects among three supporting strategies show that，compared with the anchor support，the maximum displacement，the bending moment and the stress of the roof under advanced hydraulic support group anchor coupling support are reduced by 36.87%，35.85% and 47.81%，respectively，indicating that the coupling support can effectively maintain the stability of the roof，and that，compared with the equal strong coupling support，the maximum displacement，the bending moment and the stress of the roof under the non-equal strength coupling support condition are reduced by 10.29%，12.2% and 29.53%，respectively，proving that the non-equal strength coupling support was effective to maintain the stability of the surrounding rock. Consequently，a “continuous gradient” non-equal strength support strategy of advanced hydraulic support group anchor coupling support was put forward. Laboratory tests by using the three-dimensional similar simulation experiment platform designed by the Key Laboratory of Mine Subsidence Disaster Prevention and Control of Department of Education of Liaoning Province were carried out to verify the effectiveness of the “continuous gradient” non-equal strong coupling support strategy，and the experimental results were compared with the theoretical results. It is revealed that the maximum displacement error and the maximum stress error of the equal strength support and the non-equal strength coupling support were respectively 7.2 mm and 0.48 MPa and 5.9 mm and 0.19 MPa，showing the accuracy of the mechanical model and the effectiveness of the mechanical characteristics of the coupling support considering the experiment complexity. Finally，the field application experiment of the “continuous gradient” non-equal strength support strategy was carried out in the 15106 working face of Wenjiazhuang coal mine to verify the effectiveness of the research results. The research results provide a theoretical support for the safety and stability maintenance and the intelligent support control of deep mine surrounding rocks.

During injection heat mining，new crack passages are developed between water cooling parts and natural cooling parts of hot dry rock(HDR)，and crack roughness has direct effects on seepage heat transfer of fluids. In this study，taking the granite in Daqing area of northern Songliao basin as the research object，Combined cooling of different initial temperatures(100 ℃–600 ℃) and different cooling ratios(0% immersion cooling，25% immersion cooling，50% immersion cooling，75% immersion cooling and 100% immersion cooling) were performed on the specimens(partly placed in water and partly kept in air)，and Brazilian split tests on the heat-treated specimens were performed. The cleavage surface roughness of the samples was measured by the customized roughness profiler to clarify the variations of physical and mechanical properties(apparent temperature，upshift of immersion surface，morphology and tensile strength) and cleavage surface roughness features(height difference parameters：distributions of peak and trough of height difference and mean square error of height；texture parameters：roughness coefficient and rough angle) of granite under the specific conditions. Meanwhile，the correlation of height difference parameters and texture parameters was investigated and the tensile strength-roughness fitting curve was developed. The results demonstrate that the reduction of the apparent temperature of the rock is proportional to the initial temperature and the reduction rate of the immersion cooling part exceeds that of the natural cooling part. During immersion water cooling，the immersion surface shifts upward and the upshift is proportional to the temperature. In extreme cases，the apparent temperature may be below 400 ℃ when the initial temperature is 500 ℃–600 ℃. The distribution of the height difference peak appears to the water cooling part as the temperature increases. With an initial temperature of 600 ℃ and an immersion cooling height ratio of 3/4，the tensile strength loss is maximized(84.88%，14.84%) and crack height difference parameters and texture parameters reach their extremes. The fitting coefficients of the tensile strength and roughness features at different temperatures are relatively high. However，the fitting coefficients at different water cooling ratios are relatively low but still remain in a reasonable range. Therefore，the determination of the rock cleavage surface roughness based on the tensile strength is highly feasible. This study provides references to the understanding of crack propagations in HDR in case of thermal cracking.

Aiming to optimize the retrieval and analysis process of non-destructive testing(NDT) of GFRP anchorage system in earthen sites，the Hilbert-Huang transform(HHT)method was used in this study to analyze the test signals of various intact and defective anchoring models. In the process of analysis，the noise in the detected signal is firstly eliminated by Butterworth filter method，and then the signal after denoising is analyzed by Fourier and HHT methods successively. It is found that this method can clearly identify the defect and the reflection signal of the bottom of the rod，so as to reverse the defect position and the anchorage length more accurately. It is also shown that the pulp defect length，reversed according to the difference of the reflection times of complete and defect models，is accordance with the given conditions. The results provide a systematic analysis scheme for the inversion of the anchorage length of soil GFRP and the location and size of the pulp defect.

Loess unloading collapse widely exists in test pits and the loess between piles. On the premise of ensuring the safety of loess engineering，it is important to boost the economy and rationality of engineering projects. To achieve efficient calculation of loess unloading collapse，the inverse analysis of loess collapse tests under constant loads is carried out，and the relationship between the collapse coefficient and the modulus reduction coefficient under constant loads is established according to that the deformation modulus of loess decreases after immersion. Furthermore，the inverse analysis of loess unloading collapse tests is performed，and the relationship between the unloading collapse coefficient and the modulus reduction coefficient is proposed. A modulus reduction method for calculating the loess unloading collapse is obtained with the help of the relationship between the collapse coefficient and the unloading collapse coefficient，and verified through evaluating the relationship between the unloading collapse coefficient and the modulus reduction coefficient by setting the unloading amount equal to zero to simulate the constant loading case. The research results provide an efficient way for analyzing loess unloading collapse and improve the collapse evaluation system.

In order to analyze the particle breakage of gravel-geogrid interfaces under cyclic shearing，cyclic direct shear tests of gravel geogrid interfaces were performed under different cyclic shear displacement amplitudes，normal stresses，shear frequencies and cycle numbers. The particle distribution change before and after cyclic shearing was analyzed，and the particle breakage was quantified by the relative particle breakage ratio. The test results show that the peak shear stress of each hysteresis loop increases gradually as the number of cycles increases，and fluctuates due to particle breakage. The relative particle breakage ratio increases with increasing the amplitude of the cyclic shearing displacement，obeying a hyperbolic function relationship. A logarithmic function can be used to describe the relationships between the normal stress and the shear frequency with the relative breakage ratio. The maximum relative breakage ratios are 13.54% and 11.22% for different normal stresses and shear frequencies. The number of cycles has a significant influence on the relative particle breakage ratio，and the larger the cycle number，the larger the relative breakage ratio. Moreover，under different test conditions，the relative particle breakage ratio has a non-linear positive correlation with the input plastic work，and the relationship between them can be described by a hyperbolic function.

Aiming at the calculation of the electroosmotic consolidation model with a large spacing of same nature electrodes，the concept of punctiform electrode unit is proposed according to the electrode arrangement form，and a two-dimensional electroosmotic consolidation theory of the punctiform electrode unit is established based on the punctiform electrode unit model. The governing equation is discretized by the Galerkin method，and the finite element forms of the governing equations of the electric potential field and the pore water pressure field are developed. A two-dimensional electroosmotic consolidation calculation numerical module PyEcFem is developed by Python language. The numerical simulation of the two-dimensional effect of symmetric units and the electroosmotic consolidation calculation of asymmetric units are performed，and the experimental verification is carried out. Results show that the numbers of the drainage boundaries and the number of the anodes respectively play a major role in the drainage rate in the early and later stages of electroosmosis，and that both the drainage effect and the negative pore pressure of asymmetric units are higher than those of symmetrical units. For the symmetric units，the electrode spacing ratio can describe the potential drop amplitude on the same polarity boundary，and the smaller the same polarity electrode spacing，the smaller the two-dimensional effect coefficient，obeying a negative exponential relationship. For the electroosmotic consolidation of the asymmetric units，the limitation of the drainage capacity of the units will cause the accumulation of the pore water pressure in the cathode area at the initial stage of electroosmosis，generating local positive excess pore water pressure，which will affect the change rate and the extreme value of the negative pore water pressure.

Breakage transition matrix can effectively describe the evolution of particle size distributions(PSDs) of granular soils by establishing the relationship between the percentage vectors of each size of particles before and after breakage. The existing breakage transition matrix neglects the breakage difference of different size groups of particles，and fails to consider the fact that the daughter particles fall into the original groups after surface grinding occurs. Therefore，this paper proposes a new method for calculating the breakage transition matrix. Firstly，a Hill distribution function，which can describe the distribution of the daughter particles，was proposed based on the point load test of single particle. Such a breakaed law of single particle was then introduced to DEM simulations. A series of one-dimensional compression tests were carried out based on the replacement method and the Apollo filling mode. Finally，the breakage transition matrix was obtained by tracing the breakage paths of particles with different sizes. This research overcomes the disadvantages of the traditional replacement method in numerical simulations that the choices of the number and distribution of the daughter particles are subjective，and provides an alternative approach for studying the evolution of particle breakage of granular soils.