Numerical Simulation on the Effect of Rock Joint Roughness on the Stress Field

In view of the influence of Joint Roughness Coefficient (JRC), which is for quantitative description of the joint surface roughness, on the stress field of the rock mass, compression test and shear-compression test were simulated on models with different joint roughness. The photoelasticity technique is applied to examine the feasibility of numerical simulation. The results show that numerical simulation results are in agreement with the results of photoelastic experiments. The stress concentration area is distributed near the joint plane. Thus, the joint plane controls the shear strength of the rock. In compression test, the maximum shear stress of the model is proportional to JRC and the normal pressure. In shear-compression test, when the ratio of the axial shear to the normal pressure is small, the maximum shear stress is nonlinearly positively correlated with JRC. When the ratio of the axial shear to the normal pressure is relatively large, the relationship curve between the maximum shear stress and JRC is parabolic. When the JRC is small, as the ratio of the axial shear force to the normal pressure increases, the maximum shear stress changes abruptly, and the maximum shear stress after the mutation decreases significantly. The reason is that the upper and lower parts of the model have slipped, resulting in a redistribution of stress. In addition, when the JRC is 6 to 12, it is more likely to cause stress concentration.


Introduction
Rock discontinuities not only have impacts on the stability of geotechnical structures, but also effect the geotechnical operations like rock drilling, cutting and blasting [1]. Among discontinuities such as bedding planes, faults, joints and fractures, rock joints play the critical role in the effect on the mechanical behavior [2]. Considerable efforts have been devoted to study the effect of surface roughness on the stress field of the rock joints in the past several decades. Photoelasticity technique has been widely used to reveal the distribution of the stress and stress concentrations of geomaterials subjected to various loading conditions [3]. The combination use of photoelasticity and high-definition photography is a simple and effective method for the study of stress field, owing to that it can clearly see the full-field stress distribution [4]. Photoelastic samples under uniaxial compressive were made experiments to analyze the influence of the fractal dimension on the joint intensity and joint deformation [5]. Photoelasticity was applied to study the effect of the rock joints on the shearing property and the obtained results show that the rock joint roughness has vital impacts on the shear dilatation and contact property [6].Although photoelastic tests are significantly advantageous of stress field visualization, there 2 are some shortcomings of experiment techniques. In some tests, it is inevitable that the presence of initial stress in complex models during specimen preparation will disturb the experimental results. Furthermore, there are a large number of other issues linked with the deviations between experimental observations and real results. Consequently, as the computer technology develops rapidly in recent years, numerical simulations also seem to be a powerful tool to characterize the stress field of models during loading [7]. Finite element methods could handle well the problems that experimental and theoretical methods cannot overcome. This paper focuses on the effect study of rock joint on the stress field of the rock masses by using finite element method. Photoelasticity technique is applied to examine the feasibility of numerical simulation. Maximum shear stress is chosen to quantify the effect of joint roughness on the stress field, JRC is used to quantify the rock joint roughness. The shear strength of rock masses is much lower than their compressive strength in rock masses due to the existence of joints. Therefore, this paper concentrates on acquiring the relationship between the maximum shear stress and the joint roughness of the rock mass.

Method Verification
FEM simulation and photoelastic testing were conducted to obtain the mechanical response of the same loaded disc, and the obtained stresses of the specific points by two methods were compared to examine the validity of FEM simulation. Since the verification model is relatively simple, it can be considered that there is no interference from the initial stress in the photoelastic test. Therefore, it is effective to compare the numerical simulation results with the photoelastic test results.

Determination of Fringe Constant
The material parameter, fringe constant, was determined by disk radial loading test. As a critical design parameter, the fringe constant of a photosensitive material is of great importance on the stress field analysis based on the stress optical method [3]. Principal stress difference and fringe order are related by fringe constant. Thus, it is possible to compare the experimental results with other results such as simulation results and theoretical solutions according to this parameter. From formula (1), it can be seen that the fringe constant is determined by the wavelength of the light source and the stress-optical coefficient of the photoelastic material. Therefore, the fringe constant can be obtained from a sample with a theoretical solution of stress as shown in formula (2).
Where f is fringe constant,  is light source wavelength, C is stress-optical coefficient.
where n is the fringe order of interest, t is the thickness of the sample, (σ 1 − 2 )is principle stress difference.
The photoelastic material used in this experiment is polycarbonate. The detail the specimen is shown in figure 1, and the experimental device is shown in figure 2. The light source is a monochromatic light with a wavelength of 532 nm. Figure 3 shows the dark field optical layout of photoelasticimeter.   Before loading, no fringe pattern is displayed on the disc. The fringes appear constantly from the upper and lower ends and move towards the side boundaries with the increase of load. The corresponding load was recorded under the circumstance of the integer fringes reached the center of the disk. The high definition camera captured the fringe pattern of the circular disc under load. Figure 4 shows the pictures of the radially loaded disk with one to four fringe orders. The final result of fringe constant f is obtained by averaging the four test values with a view to measurement error. The fringe constant of photosensitive material is determined according to formula (3) [8].
where is the diameter of the circular disc, is the fringe order at the center of the circular disc, and is the corresponding load.    Figure 5 shows the relationship between the loading at the center of the specimen and fringe order, where the slope of the broken line decreases after the point A. This phenomenon is caused by the specimen having plasticity at the end point. The fringe constant is determined as 7.986N/mm.

Photoelastic Testing
The experiment to obtain fringe constant is also the experiment to obtain the principal stress difference at point A and point B as shown in figure 6. The load when the fringe passes through point A and point B is recorded, thus the experimental data corresponding to the fringe order and the load are obtained. The specific data are shown in table 1. According to formula (3), the principal stress differences at point A and point B under a specific load are obtained.

Numerical Simulation
This model cannot be reduced to a two-dimensional model because the sizes of the model in three dimensions are close. Therefore, we set a three-dimensional model as shown in figure 7. To keep consistency with the photoelastic experiment, we set up a loading table with a relatively large elastic modulus to exert uniform load on the specimen. The elastic modulus and Poisson's ratio of the specimen are set to be the same as the experiment, which are 2.3 GPa and 0.39. The horizontal displacements of the upper and lower loading tables are constrained, only the vertical displacement of the loading tables is permitted, in order to simulate the actual loading process. Actually, principle stress difference cannot be obtained directly in the post-processor module. Consequently, data processing of σ x , and is shown as follows in order to convert into principle stress difference.  As shown in table 1, the data sets obtained by these two methods are closely matched. Most relative errors are lower than 4%, in addition to some outlier highlighted in red. In other words, the finite element method has been proved to provide a basis for reliable mechanical analysis.

Modelling Procedure
A series of numerical models were set to study the mechanical response of rock joints under uniaxial compression so as to obtain the relationship between JRC and the maximum shear stress of the rock. The roughness of the rock joint is simulated through changing the shape of the contact surface. JRC is calculated by formula (5).
where 0 = 10 , tooth amplitude, is the length of the profile curve. The detailed geometry property of the models is displayed in table 2. Rock type varies in nature. In this study, we select granite, one of the most widely existed rock, as research object. However, the density, elastic modulus and Poisson's ratio of granite are not certain due to the environmental complexity, these values are in a range. For the convenience of research, we set the density of granite as 3×10 9 t/mm 3 , the elastic modulus as 75GPa, and the Poisson's ratio as 0.2. Meanwhile, we only consider the relationship between JRC and the maximum shear stress in elastic  The accuracy of boundary conditions and contacts has a great impact on the reliability of numerical simulation. Complex boundary conditions and contacts are not easy to control in modeling, which may lead to the inaccuracy of numerical simulation [9]. In this study, the boundary conditions of the model are established to make the boundary conditions of the model as consistent as possible with the actual engineering onditions. The detailed boundary conditions are as follows:  The maximum pressure applied on top and bottom surface of the model is 10MPa.  The horizontal and vertical displacements of top surface and bottom surface of the specimen are both fixed to enable the uniaxial compression process.  The type of sample surface interaction is defined as a surface-surface contact with finite sliding.
Due to the roughness of the surface of the two specimens, the tangential behavior of contact interaction was set as punishment, and the friction coefficient was 0.4. The normal behavior of contact interaction is set to hard contact, and separation is not allowed after contact. Due to the discontinuity of specimens, the mesh of finite element model is of great significance. It seems impossible to divide complex models in a high-quality structured way, in which case unstructured grids are a better alternative [10]. Considering the complexity of the model, we combine uniform grid with adaptive grid refinement to improve the performance of simulation.   figure 9 shows, the element sizes are different for different regions. The reason is that the adaptive grid partitioning technology can subdivide the grid according to the stress level to meet the requirement of computational accuracy. In other words, stress concentration occurs in a dense grid. The simulation results show that the stress in joint surface area is higher than other areas, and the stress field distribution characteristics have not changed significantly with the increase of JRC.  Firstly, it can be determined from figure 10 that the maximum shear stress is proportional to normal pressure. In particular, different roughness, different scale constant. Secondly, the relationship between the maximum shear stress and JRC will be discussed. Figure 11 illustrate the fitting curves of maximum shear stress of the models with different JRC under uniaxial compression. One surprising finding is that the maximum shear stress is proportional to JRC at a given load, although it is well known that the maximum shear stress increases with normal pressure.

Modelling Procedure
To be consistency with the former simulation, there are not any changes in the geometry and material properties compared with former models. Considering the actual situation and the convergence of the model, the boundary conditions of the model under shear and compression are more complicated. The detailed boundary conditions are as follows:  The maximum pressure exerted on the top surface of the specimen is 10MPa.  Horizontal uniform loads of 20MPa were applied to the left and right sides of the specimen to obtain the mechanical response of rock joints under shear loads.  The horizontal displacement of the top surface of the specimen is fixed to exert the axial compression.  All the freedoms of the bottom surface of the specimen is constrained.  The type of surface interaction between two specimens is also defined as surface-surface contact with finite sliding. The tangential behavior of contact interaction is set as penalty and the friction coefficient is 0.4 due to the roughness of two sample surfaces. The normal behavior of contact interaction is set to "hard" contact, and separation is not allowed after contact. For the sake of improving the precision of the simulation results, adaptive mesh refinement methods are also applied in these models.

Results and Discussion
The simulation results are similar to the uniaxial compression testing, the stress in joint surface area is higher than other areas, and the stress field distribution characteristics have not changed significantly with the increase of JRC.To begin with, shear-compression ratio n is defined as follows.   Next, cases where is more than 0.6 will be discussed. There are two maximal points for all curves and the maximal points are located at the position where JRC equal to 6 and 12. When JRC is less than 6, the maximum shear stress increases with the increase of JRC, and when JRC is greater than 12, the maximum shear stress decreases with the increase of JRC. It can be seen that when JRC is greater than 6 and less than 12, the joints are more prone to stress concentration. The tendency of the curves shown in figure 14 differ from the curves shown in Figure 13 when JRC is greater than 15. In uniaxial compression testing, we get a conclusion that the maximum shear stress is proportional to JRC under a given load. Thus, when n is quite small and JRC is greater than 15, the maximum shear stress is nonlinear proportional to JRC.  According to figure 15, we can determine that the maximum shear stress is proportional to shear pressure when JRC is greater than 3. The mutation of the maximum shear stress when JRC equals 3 and shear pressure equals 1.4 MPa is shown in figure 16. The reason is that there is a relative slipping between the top part and the bottom part of the model. Whether occurred slipping is determined by JRC and . There is a threshold that slipping will happen when JRC is lower than the threshold and is big enough. Conversely, slipping will not happen when JRC is bigger than the threshold. The phenomenon is referred as self-locking.

Conclusion
The following conclusions were drawn:  The stress in joint surface area is higher than other areas, and the stress field distribution characteristics have not changed significantly with the increase of JRC. The shear strength of the joint controls the shear strength of the rock.  Under uniaxial compression, the maximum shear stress of the model is proportional to JRC and the normal pressure, and the greater JRC, the greater the proportional constant.  Under shear-compression, when n is relatively small, the maximum shear stress is positively correlated with the JRC nonlinearity. With the increase of JRC, the positive correlation gradually weakens. When n is relatively large, the relationship curve between the maximum shear stress and JRC is parabolic. When JRC is small, as n increases, the maximum shear stress changes abruptly, and the maximum shear stress after the abrupt change decreases significantly.  Under the same loading conditions, when JRC is greater than 6 and less than 12, the joints are more prone to stress concentration.