An efficient three-dimensional cell-centered Finite Volume Method(ccFVM) is developed for single- phase fluid flow in fractured porous media. It works on arbitrary grids and is parallelized through OpenMP. With pressure centered in the cell，the ccFVM employs spring-in-series model for space discretization and explicit difference scheme for time discretization. It uses dynamic relaxation technique for element-by-element iteration. Numerical tests indicate that the ccFVM is as accurate as FEM，but more efficient. The parallel procedure obtains a speedup of 4.0 on CPU i7-3770 and a speedup of 4.2 on CPU i7-4770. High parallel efficiency of over 50% is achieved on both machines.

For the geological disasters of slope failure under seismic loads，the failure criterion has been a hot and difficulty issue in rock mechanics. The main objective in the study is to combine the limit analysis and wave theory to carry out a new evaluation of slope stability from the perspective of power，that is，the rate of external work is equal to the rate of energy dissipation. By means of the methods of strength reduction and overload factor，the factors of safety were adopted to evaluate the quantitative safety of the rock slope. The oblique incident P wave across a single joint in jointed slope was considered. The variations of the critical height of jointed slopes and the transmission and reflection of stress waves were investigated. Then，within the critical height of jointed slopes，the variation of the critical height was studied under the stress wave. Finally，based on the energy criterion that the rate of external work is equal to the rate of energy dissipation，the variations of FOS were investigated by means of the methods of strength reduction and overload factor. The sliding energy of the block in limit state was also calculated. The investigations revealed that the higher the cohesion of the joint，the higher the self-stable height of the jointed slope. Therefore，the higher critical heights of the jointed slope during earthquakes result in the significant improvement of the seismic performance of the jointed slope. As the peak strength of the incident P-wave increases，the reduction degree of the critical height of jointed slope also increases. As the incidence of the stress wave decreases，the variation degree of the critical height of the jointed slope also gradually increases，indicating that the failure possibility of the slope is also increasing during earthquakes. The variation of FOS is just opposed to that of the transmitted waves. The higher the rock slope is，the more unstable the block is. Therefore，the corresponding sliding energy is decreasing as the jointed slope is in limit state.

The three-dimensional morphology of a rock joint have a great impact on its shear behavior. In order to study ofrelationship between the three-dimensional morphology characteristics and peak shear strength of rock joints，fresh tensile rock joints of granite and sandstone are prepared anddirect shear tests are carried out under different constant norm loads. The morphology of each jointwas measured before the direct shear tests.The three-dimensional morphological parameters of each joint were calculated. Then，by analyzing the test results，a new peak shear strength criterion of rock joints is proposed using twothree-dimensional morphological parameters which are termed as the maximum apparent dip angle and a roughness parameter C. Direct shear test data of tensile rock joints collected from Grasselli are used to further verify the proposed criterion. A comparison among the new criterion，Grasselli’s criterion and Xia?s criterion are made from the perspective of both the rationality of the formula and the prediction accuracy for the 50 groups of joints. It turn out that the new criterion have some advantages in predicting peak shear strength of rock joints.

We carry out different angle shear tests and fatigue tests，to research salt rock compression-shear mechanical properties and crystal dislocations characteristics. Based on the damage mechanics theory，dislocation-damage relationship model is established. The main conclusions are：(1) Because that salt rock crystal dislocation slip resistance decreases，salt rock deformation resistance，peak stress and maximum deformation decrease with increasing shear deformation angle. (2) Fatigue life is associated with cyclic stress limit level. At different shear angles，if stress level is the same，the fatigue life is substantially the same. If the stress value is the same，the greater the angle，the shorter fatigue life. (3) Reverse shear influences hugely the fatigue life of rock salt，with more than 90% decrease. (4) Salt rock fatigue crack formation is related with movement of dislocations squeeze. Crack growth is the result of alternating slippage for multiple slip systems. The dislocation-damage model is established，so that the damage can essentially reveal the destruction process of salt rock. The greater the fatigue damage results in more obvious decrease of the carrying capacity.

A finite element method(FEM) is proposed to solve the excavation problem in sand. Based on the existing research achievements of stress release and elastic-plastic unloading theory of bulk solids，the non-linear characteristics of load-displacement curve is equivalent to the gradually release process of the excavation load，during which the loading modulus is the same with the unloading one. The excavation load considering the internal friction resistance of material is deduced，after introduce the octahedral shear stress to calculate the internal friction force according to the existing research achievements，the relationship between the internal friction resistance of material and the residual stress during excavation is established according to the virtual work principle；and the FEM of stress release，which takes the excavation load into consideration，is carried out by adjusting the unloading stiffness of material during excavation on the basis of stress release properties. The derivation of node equilibrium equations indicates that the excavation load has influence on the whole computational domain and also has corresponding relation with the stress release effect；A static loading-unloading model test of sand was designed and simulated by the proposed FEM as well；the residual stress calculated by this FEM was basically identical with the results of the test，and can be also explained by the existing stress release research results. Therefore，this FEM provides the calculation of excavation problem with a numerical simulation method.

The mechanism of rockburst is intricate，which is controlled by factors such as lithology，stress and excavation conditions，and the physical simulation test for rockburst has some unique advantages. On the basis of summarizing previous studies on rockburst with physical simulation test，the shortcomings in the present study and the future direction of study were described in detail from the aspects of the test instruments，the similar materials，the loading conditions and the methods of excavation，supporting and monitoring. The following conclusions are obtained. Establishing the similarity criterion of brittleness suitable for rockburst and developing similar materials with low strength and high brittleness are basic testing requirements for the successful physical simulation of rockburst. Solving the problems of the similarity criterion and the manufacturing technology for rigid structural plane theoretically and experimentally respectively are keys for the structural rockburst research. With the required loading accuracy and stressing uniformity，designing the large-tonnage loading device and the TBM which considers the complex rock-machine interaction，and the drill and blast construction technology，are guarantees for the study of different factors influencing the rockburst. Devising diverse information monitoring systems that are matched with the typical deformation and failure characteristics of hard rock cracking is prerequisite to successfully interpret the behavior of rockburst.

In order to investigate the mechanism of grouting strengthening the fault gouge of Yonglian tunnel，a test system including a module of grouting technology，a module of injection medium and a module of information collection was designed. The orthogonal uniaxial compression tests，tests for the typical characteristics of pulp-rock interfaces and the scanning electrum micrographs(SEM) tests for fault gouge were carried out. The uniaxial compressive strength of the fault gouge was increased by 181%–2 535% after grouting. The strengthened mass was damaged through the process of the cracking of pulp veins，the breaking of the pulp-rock interfaces and the failure of the soil grout composite. The fault gouge was strengthened by the skeleton strengthening mode，and the key factor for raising the grouting effect was the grouting pressure. The cohesion of the pulp-rock interfaces was increased by 93%–274% after grouting. The influences on the strength of pulp-rock interfaced under the different grouting pressures，materials and initial dry densities were analyzed. When the cement slurry was used，the cohesion was increased by 125%–148%. The fault gouge was transformed from the floc structure into the overall densified structure after grouting. The pulp-rock interface was composed of the cementation plane，the transition area of permeability and the transitional area with micro factures. The direct and indirect strengthening modes of the pulp-rock interface are defined.

Tests on red sandstone under shear creep condition were done with both self-developed shear creep microscopic testing device and monitoring software of coal or rock micromechanics properties. Micro-cracks temporal and spatial evolution law of red sandstone under shear creep condition and the relationship between morphological properties of rupture surface and shear creep strength were analyzed. Slope slip mechanism is discussed in this paper，too. The results show that creep deformation and micro-cracks propagation were developed simultaneously under shear creep condition；the propagation direction of micro-cracks has bias in some degree compared with shear stress，and multiple bifurcations occur in the process of crack propagation；most cracks are along grain boundary，and the irregularity of crack propagation path is affected by the distribution of rock grains；micro-cracks are mostly formed by tension and shear together and the differences of stress response and inconsistent deformation between different components of rock and the interaction of micro-cracks promote the initiation and propagation of micro-cracks. Relation curve between roughness of rupture surface and shear creep strength is non-linear increasing，and the curvature increases at first then decreases.

Theory of structure stability was applied for the problems of rock mechanics. Two simple models were adopted for the analysis of reversed fault earthquake and rockburst during shaft excavation，and their analytical solutions were obtained. The equilibrium path curve for the shaft excavation links the non-dimensional displacement with the pressure of borehole wall，while the equilibrium path curve for the reversed fault earthquake links the non-dimensional displacement at far-field with the constrained force. The instability of the rock structures was described as the instability at the point of maximum value. The rockbursts and reversed fault earthquake are two different kinds of instabilities. The rockburst during the shaft excavation is the instability occurred at the maximum stress and with a displacement jumping，while the reversed fault earthquake is the instability occurred at the maximum displacement with a force jumping.

Firstly，in order to validate the DEM software UDEC in simulating the wave propagation across a rock mass with two sets of nonlinear joints，two examples were computed with a nonlinear joint model embedded in UDEC. The results agreed well with the existing theories. Secondly，the influence of the angle and space of joints on the propagation of stress wave across two sets of intersecting joints was studied. It was found that the attenuation of stress wave across two sets of intersecting joints was greater than that across the parallel joints under the same angle and joint space. Transmission coefficient reduced with increasing the angle. The space of the intersecting and parallel joints had the similar influence on the transmission coefficient. Finally，a parametric study was conducted to investigate the effects of the angle ? between the joints and the stress wave，the angle ? between the two sets of joints，the space for each set of joints and the width of joint area L. The analytical results showed that the influences of ? and ? were different with the ranges changing. The similarities existed between the influence of space of intersecting joints and parallel joints. With the increasing of L，the attenuation decreased significantly firstly and then slowed down. Hyperbolic function was employed to fit the result and the fitting result had good agreement.

Grouting technique is widely used to control the water inrush hazards and increase the strength in tunneling and underground engineering. Therefore，the anti-washout property of grout under the flowing water condition is the key factor that determines the grouting effect. In order to evaluate the anti-washout property and the gelling performance of the fast-gelling grout，a specific test apparatus was designed to simulate the vital process of the grout being flushed by the flowing water before gelled. With this device，a series tests were carried out by controlling the velocity of water flow and the mixing proportion of grout. The grout retention ratio(GRR) index was proposed to evaluate the anti-washout property quantitatively. The test results showed that the grout of cement-sodium silicate achieved a large grout retention ratio when the water flow velocity was small，however，the grout retention ratio decreased significantly with the increasing of the velocity of water flow. The GRR of the grout of cement-sodium silicate was influenced by the mass ratio of water to cement and the volume ratio of cement to sodium silicate of the grout. The modified polyurethane grout had a larger GRR under the fast water flow condition and the gel was more compact，which meant that the modified polyurethane grout had the better anti-washout property and the gelling performance. According to the different anti-washout property and gelling performance of grout，a method of choosing grout materials reasonably to dealing with the water inrush was proposed.

The stability of high rock slope is one of the important issues during the construction of hydropower station. The stability of the left bank slope at Jinping I hydropower station was the key issue due to its high geo-stress，developed faults and fissures，large-scale excavation，deep unloading and complicated geological conditions. The characteristics of the microseismic activity of the left bank slope were introduced firstly. Analyses of the temporal and spatial distributions of the microseismic activity in the processes of excavation-induced unloading and under the excavation conditions during the different periods and in the different regions were conducted. The correlation among the initiation，propagation，extension and coalescence of micro-fractures of rock and the rock mass damage were discussed. The correlation between the microseismic activity and the construction response at the left bank slope was also obtained. The results showed that the microseismic activity induced by the unloading relaxation of the deep rock mass as results of the excavation of the foundation pit of dam，the X shafts for lamprophyre vein replacement and the drainage tunnels at the elevation of 1670 m and the consolidation grouting of weak layers inside the rock slope reflected in real time the construction conditions on site very well.

Mechanical experiments of the limestone samples from Longmen grottoes under the action of coupled chemical solutions and freezing-thawing processes were carried out to investigate the infiltrated water erosion and freezing-thawing damage on the surrounding rocks in Longmen grottoes. The results of the experiments demonstrated that the mechanical damage of limestone was increased gradually as with increasing freezing- thawing cycles under the coupled chemical solution and freezing-thawing processes and the damage was greater than that under the erosion of chemical solutions only. For example，for the cases of the freezing-thawing cycles of 90 under the coupled chemical solution and freezing-thawing processes，the strengths of the specimens soaked in the distilled water，Longmen water and NaCl solution dropped 50.73%，54.92% and 57.67% respectively compared with the specimen strength in the natural state. While under the erosion of chemical solutions without freezing-thawing，the strengths was dropped only 21.58%，22.88% and 28.72% respectively. The freezing-thawing deterioration model of limestone under investigation was found to be all due to the particle loss model. The abundance of the condensation nuclei and the pH value of water solution were found to be the important factors affecting the degree of damage to the limestone under the coupled chemical solution and freezing-thawing processes. An equation describing the erosion damage to limestone under coupled action was established.

The artificial boundaries in numerical analysis may introduce fictitious scattered waves，which inevitably affect the simulation results. The artificial boundary problem in the discontinuous deformation analysis (DDA) method had been investigated in detail. Firstly，the contributions to the DDA equations for nonreflecting boundaries(including the viscous boundary and the viscoelastic boundary) were deduced based on the Newmark method. Secondly，to accurately simulate the motion due to external source wave such as earthquakes，a free-field boundary in DDA was introduced. Lastly，the unified static-dynamic boundary was implemented in DDA to ensure the consistency of the transformation of the boundaries. The analysis of examples with the improved DDA program showed that each artificial boundary theory has been correctly implemented. The viscous boundary was found to be highly effective for the absorption of the reflection wave at the artificial boundaries，and the viscoelastic boundary adequately simulated the elastic recovery of the infinite domain. The DDA model with the free-field boundaries at both sides is better than that imposed by the viscous boundary. The DDA with the unified boundary can completely simulate the static and dynamic calculation process.

The weak layer in a slope always controls slope deformation and failure mode. Under an earthquake，the complex interaction mechanism between the weak interlayer and the seismic waves hinders the understanding of the dynamic responses of this type of slope. Based on the existing knowledge，shaking table tests were conducted for four slope models containing a weak interlayer. The effects of the weak interlayer were uncovered by comparing the responses of the models with and without a weak interlayer. Both the horizontal(PHA) and vertical(PHA) acceleration responses were amplified along the altitude for all four slope models，especially in the locations with the elevation higher than the interlayer. In the locations with the elevation lower than the interlayer，the PHA amplification factor remained within 1.5. In comparison with the homogeneous model，the existence of an interlayer resulted the two component responses in the upper part of slope strongly depended on the features of the interlayer and the input intensity. When the input intensity was weak(≤0.3 g)，the responses for the models with a thin or a thick interlayer were stronger than in the homogenous model. With the increasing of the input intensity，the thick interlayer generated an effect of seismic isolation，which caused the PHA and PVA responses at the slope crest decreasing by 50% and 70% respectively. The horizontal displacements on the slope surface were amplified and reached the maximum value at the crest. The initial deformation occurred at the top of for the homogeneous model slope. However，for the four models with a weak interlayer，the deformation was firstly observed in the position above the interlayer and below the crest. The location and the degree of the slope failure were closely related to the thickness and the dipping angle of the weak interlayer. The model with a horizontal interlayer caused the failure location higher than the models with an anti-dip interlayer and the model with a thin interlayer caused the failure more severe than the one with a thick interlayer.

The nonlinear equations of gas seepage in coal considering adsorption effect was established based on the relationship between the permeability of coal and the quantity of adsorbed gas from the earlier experimental results. The dimension and the physical meaning of the parameters in the derived equations are clearer. Darcy?s law is a special case of the derived equations. An analytical solution of the gas seepage equation was derived. This analytical solution reflected the migration of the gas in coal seams faithfully. There was a good agreement between the results from the new model and the experiment results of gas flow in coal sample.

The variations of flow，chemical fields and aperture in a fracture could be described by the model of coupled fluid flow and chemical dissolution，which were affected by the dissolution rate and the initial distribution of aperture. Two factors，the surface reaction and the mass transfer，were considered to control the dissolution process. A rough fracture was generated numerically with the fractal theory and a model of coupled fluid flow and chemical dissolution was thus established and numerically solved. The results showed that the distribution of aperture，flow velocity and concentration of  appeared as inhomogeneous wormholes. The flow through the fracture increased during the whole period of simulation，while the growth rate of flow increased sharply at a certain moment，called the breakthrough time. The factor controlling the chemical reaction varied during the dissolution，was the surface reaction initially and then was turned to the mass transfer on the upstream and extended to downstream gradually.

The nonlinear elastic nature of surrounding soil of pile was simulated with a hyperbolic model，and the time effect of the settlement and the distribution of pile shaft resistance and the axial force of a single pile under static loading were investigated. The pile shaft was divided into finite segments from the head to the toe. The settlement of any pile segment at any moment was deduced with the program of wave equation analysis based on the decomposition of traveling wave. The parametric study was conducted to analyze the influence of pile-soil parameters on the settlement behavior of a single pile. The theoretical calculated s-lgt and q-s curves were compared with the testing results. The results demonstrated that the solutions simulated well the settlement behavior of a single pile under static loading.

The formation and the influencing factors of seepage erosion in soil and the relationship between the damage to the soil structural strength and the erosion ratio were systematically analyzed. The relationship of the constituents of soil grains，the constitutive equation and governing equations of seepage erosion in the earth structure were established. A relationship of structural strength damage of soil due to seepage erosion was presented. The partial differential equations for the coupled multiple fields based on the governing equations were solved with the finite element method. A core dam was analyzed. The results of the numerical analysis validated the governing equations and the relationship between the damage elastic modulus and erosion ratio.

The ground motion in earthquake is a complex and multi-dimensional movement. The ground soil experiences the horizontal shearing in two directions with the shear directions and the magnitudes of the loading both coupling. The effect of cyclic shear direction and shear strain amplitude ? on the shear modulus G and the pore pressure u were therefore investigated to the saturated soft clay during the shear process and compression in the reconsolidation process with the VDDCSS system. The shear modulus G was found to be different in different directions due to the soil anisotropy. The shear test with small ? made G increased slightly in other directions. The changing of shear direction accelerated the declining of G(1.51 times of unidirectional attenuation of G) and promoted the development of u in the shear process(1.67 times of the increment of the unidirectional pore pressure u). The shear tests with ? greater than a value ?s had no such phenomenon. The changing of shear direction influenced greatly the amplitude of pore pressure uamp. When 0.1%＜?＜2.0%，the reconsolidation settlement ?Hrc/H increased with the increasing of ?. The changing of shear direction promoted the ?Hrc/H at the small ? and reached a value which was reached at the large ? in uni-directional shear；when ?≥2.0%，the ?Hrc/H under 2 kinds of shear-mode tended to be the same.

The bearing capacity of the square footings on the reinforced sand slope was studied with the laboratory modelling tests. The effects of the location and the number of the reinforcement layers in the slope on the bearing capacity and the failure mechanism were investigated. The experimental results indicated that there was an optimum location of reinforcement with the highest bearing capacity under the different numbers of reinforcement layers. The optimum depth of the first reinforcement layer decreased，whereas the efficiency of the sand- reinforcement system increased with the increasing of the number of reinforcement layers. According to the characteristics of the load-settlement curves，the deformation and the reinforcement rupture，the failure mechanisms of the reinforced sand slope were classified into three types，the failure above the top layer of reinforcement，the failure between the reinforced layers and the failure below the reinforced zone. Meanwhile，the failure patterns were inferred.

The classic theory of earth pressure can only be used to calculate the limit earth pressure on the retaining wall under the mode of translation. With the Mohr circle expressed with the principal stress difference and on the basis of the radial stress-strain relation obtained from the stress path triaxial tests，the nonlinear relationships of  the mobilized internal friction angle of the soil and the friction angle between the wall and the backfill soils vesus the displacement were established under non-limit state. The method of effective area ratio of wall displacement was put forward to quantify the relationship to rigid retaining wall under the deformation mode of rotation. Based on the method above and combined with the method of horizontal slices and the modified Coulumb?s formula，the formula of earth pressure distribution，the resultant force and its application point under non-limit state which considered the influence of wall deformation were deduced. The result of an example showed that the difference between the theoretical calculation value quantified by the method of effective area ratio of displacement and the measured value was small. The proposed formula well reflects the change of earth pressure with the wall displacement and is an effective extension of Coulumb?s theoretical formula.

The end-bearing piles do not meet normally the deformation compatibility between piles and soils. The deformation adjustors were thus designed to make the deformation of piles and soils compatible. The deformation adjustors were placed on top of the piles so that the end-bearing piles and soil layers together form a composite foundation. The foundation soil carried the most of the superstructure load and the purpose of the settlement control of the building was achieved with the end-bearing piles carrying the rest of the load. This method was applied in practice where the bearing capacity of natural foundation could meet the requirements but the settlement of the building was too large. Using the method presented in this paper，the foundation soil bore 55% of the superstructure load and the settlement was controlled well. A series of monitoring data showed that the scheme of the end-bearing composite pile foundation based on settlement control was reasonable and effective.