Abstract
The rock mass can be assumed to be homogeneous material from a macroscopic view; however, it is the heterogeneous material in mesoscopic scale and its physicomechanical properties are discontinuous in space. The failure of jointed rock mass was usually caused by the initiation, propagation, and coalescence of new wing cracks derived from primary joint. In order to further study the rock fracture instability, we need to study the expansion of rock cracks under external loads from the macro-meso perspective. This paper, based on the manifold cover concept, proposes a new discrete element numerical method, manifold particle discrete (MPD), combined with the particle contact model and the introduced concept of stress boundary. The proposed method can easily simulate the crack generation, propagation, and coalescence of jointed rock mass from the macro-meso perspective. The whole process of rock fragmentation is thereafter reproduced. By analyzing the manifold cover and sphere particle model, this paper constitutes the sphere unit cover function of three-dimensional manifold cover, establishes tetrahedron units, and obtains the equilibrium equation and compatible equation of the MPD model. For rock-like brittle material, crack propagation process can be simulated, and it also verifies the accuracy of the proposed numerical method.
1. Introduction
The rock is a complex mixture composed of various mineral crystals, cements, and pore defects. Usually, a large number of geological faults such as cracks, fracture surfaces, joints, holes, and fillings occur in them and they are randomly distributed in the rock. Under the action of the external load, the nucleation and expansion of rock internal microdefects and interaction between them determine the rock’s deformation and fracture characteristics [1]. Therefore, we need to study the expansion and coalescence of rock cracks from macro-meso and multiscale perspectives to know substantive characteristics of rock failure and instability. Numerical simulation of this issue has been one of the hot and difficult fields of rock mechanics numerical calculation [1–5].
Along with the rapid development of high-performance computers and the theory of numerical calculations, several models and software based on numerical analysis methods have been applied in the field of geotechnical engineering to simulate the mechanical response and failure modes of materials. Some numerical methods such as RFPA, DDA, meshless method, manifold method, boundary element method, and discrete element method were applied in the study on crack propagation, and good results have been achieved (Cundall [6], Liu et al. [7], Harrison et al. [8], Chen et al. [9], Li et al. [10], Tang et al. [11–14], etc.).
FEM, a mature application numerical method, has been widely used in numerical research and engineering [15]. For discontinuity problems such as large cracking, sliding, and separation issues, it is, however, still quite restricted. The FEM requires quite strict grid conditions; if mechanical problems contain quite complex boundary conditions or geometry, especially in three dimensions, meshing is a very tough job. The meshless method recently proposed tries to solve the problem, providing a promising solution for crack propagation and rock failure [16, 17]. In addition, the finite difference method which uses regular grid to process crack propagation is limited in calculating geotechnical engineering problems, especially when the homogeneous materials or complex boundary conditions are considered [18, 19].
The boundary element method was developed in the 1960s. It regards the border issue as the boundary integral equation and solves the numerical solution by dividing the unit on the boundary, thereby obtaining the field variables at any point within the region. For rock failure problem, the BEM proposes that the problem domain can be divided into a few numbers of subdomains according to the crack propagation path and the hypothetical crack interface. The main advantage of the boundary element method is to reduce the problem dimension; then it is easier to preprocess the data than to preprocess the finite element method and the finite difference method. However, the boundary element method cannot handle the inhomogeneity problems, and it is less flexible and effective than FEM on noncontinuous multimedia problems and nonlinear problems [20–23].
Once Doctor Shi [24] proposed the numerical manifold method in 1997; this method immediately attracted extensive interest and attention of international academic circles. The greatest advantage of this method is that it can uniformly simulate rock continuous and discontinuous deformation and, meanwhile, couple and expand the traditional finite element method, discontinuous deformation analysis method, and analytical method. The core manifold cover concept provides new analysis thinking for rock numerical simulation and can be widely applied to solve continuous and discontinuous problems of solid, liquid, or gas. For heterogeneous media that are not totally continuous or totally discrete, Doctor Shi created a new numerical calculation method, numerical manifold method, using the finite cover technique in manifold analysis based on his proposed discontinuous deformation analysis in the 1990s. This method uniformly solves mechanics problems of continuous and discontinuous deformation. Currently, the research and application of manifold method are mostly limited to two-dimensional models and there are relatively fewer theoretical research and programming of three-dimensional manifold method [25–30].
Discrete element method is a numerical method for analyzing the mechanical behavior and the movement of granular materials and was proposed the last century by the Cundall [6]. Because of the many advantages of discrete method and the limitations of the finite element method, discrete element method is inherited and developed continuously; after nearly three decades of efforts of countless scholars and research institutions, the discrete element method has gradually developed [1, 2]. Rock failure process of brittle materials like concrete is a process by a continuum to discontinuities; there are some shortcomings and difficulties to simulate its micromechanics using the finite element method in general. Numerical methods used for granular material analysis, such as the discrete element method, have great advantages in dynamic analysis after structural failure and therefore many researchers adopt the discrete element method together with other methods to simulate dynamic behaviors [6–8].
For complex rock, there is no good numerical simulation method that can truthfully describe the actual behavior of the rock mass and especially for rupture process of complex rock which is not effectively analyze. With the combination of micromechanics and statistical strength theories and continuum damage mechanics and damage mechanics, in recent years, some of the numerical simulation appeared to be based on meso-structure considerations in the aspect of rock and concrete fracture analysis (Cusatis et al. [31, 32]). Based on the manifold cover concept, this paper proposes a new discrete element numerical method, manifold particle discrete, to simulate the generation, propagation, and coalescence of rock crack from the macro-meso perspective, combined with the particle contact model and introduced concept of stress boundary.
2. MPD Model Unit
2.1. Internal Structure of Rock Mass
In the long course of the study, people tend to overlook the complex internal structure of brittle materials, and the material is generally regarded as the macrohomogeneous continuum for the description of the macroscopic constitutive relation of the material; that is, the material is assumed to have the same properties of each point. In fact, stress state of brittle materials (such as concrete) in the microscopic and mesoscopic level and mechanical properties reflected in the macroscale are quite different; macroscopic failure process of brittle materials is closely related to the heterogeneity of their microscopic. Therefore, for brittle material (such as rock) in a mesoscopic scale, the research of the whole process of crack propagation, penetration, and breaking of rock mass will help to further understand their macroscopic damage and their strength and deformation characteristics.
The rock is a complex mixture composed of various mineral particles, cements, and pore defects. As shown in Figure 1, take out a representative unit from the brittle materials, such as rock, and analyze rock the structure from the macro-meso-scale.
2.2. Generation of the MPD Model Unit
As shown in Figure 2, the microstructure of the rock, mineral grains will be equivalent to the ball units, and the sphere unit is assumed to be physical cover. Tetrahedral elements forming the force between mineral grains cements define mathematical coverage. Combined with the sphere grain unit and manifold unit, the center and the midpoint of each side of the unit are connected to form the MPD unit, as shown in Figure 2.
(a) The microstructure of the rock
(b) The equivalent unit of rock macro-meso structure
(c) Three-dimensional MPD unit
(d) Dimensional MPD unit
3. Balance Equation and Geometric Equation of the MPD Model
3.1. Displacement Analysis
Tetrahedron is divided into four regions from the center point, and is one of them (shown in Figure 3). As shown in Figure 3, assume that the center coordinate matrix of each sphere unit is ; each sphere unit has 6 degrees of freedom to define the displacement field of sphere particle using rigid body dynamics:
(a) Tetrahedron
(b) Interactions between particles
In the preceding formula, and are translational displacements of particle , while is the rotational displacement.
The general cover of manifold method is composed of mathematical cover and physical cover. The cover function and related weight function are defined based on the general finite cover and the overlapping portion is transformed into conventional units. Then, solve the established global equilibrium equations of general finite cover. Based on manifold cover, the tetrahedron formed by connection between ball particles is regarded as mathematical cover and its related weight function is . Then,
In the preceding formula, are 16 constants related to node coordinates of manifold units.
The displacement function of the manifold unit is the weighted average of displacement functions of 4 physical covers; that is, they are connected through weight function . Then, the overall displacement function of the manifold unit can be expressed by
In the preceding formula,
3.2. Strain Analysis
As shown in Figure 5, the connection length between any two ball particles in the MPD unit is . Formula (3) is used to define the displacement of the center of a regional surface of the tetrahedron, and the displacement can be expressed by . Taking formula (3) into consideration, the strain vector of the manifold unit is defined, is decomposed into normal and tangential parts, and is expressed by
In the preceding formula, , , and are mutual perpendicular direction vectors on the connecting line between particles, ; , , and , ; is the coordinate of point which is midpoint of a certain side .
4. Numerical Results
The MPD model computing framework described in this paper will be implemented in the open-source discrete element software YADE (as shown in Figure 4). The YADE computing part is programmed mainly by C++ and Python and is independently implementing new algorithms and interfaces are allowed. It can be applied to fast and concise simulating control, after processing, and debugging. YADE describes all force imposed on motion units with Newton’s second law and the numerical solution of system dynamic characteristics is conducted by using the time algorithm. In each time step, the speed and acceleration are constant. This system reproduces the evolution by using the explicit finite difference algorithm. The main feature of the discrete element is that it can make a large number of particles produce nonlinear interaction without requiring too much memory; so it will not be affected by convergence problems.
However, the explicit algorithm is not unconditionally stable and requires accurate evaluation of the stability of simulated numerical values. In the elastic region, the stable condition can be expressed by , where refers to the maximum inherent frequency of the computing system. The stable time step of the MPD model can be estimated by calculating the inherent frequency of each MPD model tetrahedron. This requires resolving the characteristic value problem , where is the rigidity matrix, while M is the quality matrix.
The above method is used to calculate the rock failure analysis, and its flowchart is shown in Figure 5.
Numerical Example 1 (simulation of the uniaxial tensile test for the intact rock). The cylindrical sample size is mm; the number of ball particles is 2500; the maximum diameter is 3.5 mm; the physical computing parameters of the computing model are listed in Table 1. Compared to the particles, in this paper, we impose force (tensile strength of cements) between particles to simulate the acting force of cements. Based on YADE platform, using ParaView triangular mesh, the calculation model is as shown in Figure 6.
As the sketch of samples after failure, Figure 7 shows that shear failure occurring on samples. As the comparison chart of the sample experimental and computing results, Figure 8 shows that curves change similarly and the peak strength is close, verifying the accuracy of the proposed numerical method and feasibility of rock failure analysis.
Example 2 (simulation of the uniaxial tensile test for the jointed rock masses). The number of ball particles is 2800; the maximum particle diameter is 3.2 mm and the average diameter is 2 mm. The model contains a jointed, and its inclination is 45°. The physical computing parameters of the computing model are listed in Table 1, and the calculation model is as shown in Figure 9.
Further evidence of the excellent modeling capability of MPD is provided in Figure 10 where contours of meso-scale crack opening at failure are shown, and it shows crack path of rectangular model under tension for various time steps: (a) step 135, (b) step 155, and (c) step 175. Figure 11 is damage cloud picture for sample after tensile damage. Velocity vector is shown in Figure 12; velocity is relatively large at the ends of the joint under the condition of tension.
(a) Step 135
(b) Step 155
(c) Step 175
Example 3 (indoor uniaxial compression test is simulated for the horizontal jointed rock masses). The specimen size is the same as the above example 1; the maximum particle diameter is 3.5 mm and the average diameter is 2.3 mm. Physical parameters of the calculation model are shown in Table 1. The calculation model of the horizontal jointed rock masses is as shown in Figure 13 and the calculation process of the horizontal jointed rock masses is shown in Figure 14 for various time steps: (a) step 100, (b) step 120, (c) step 160, and (d) step 250.
(a) Step 100
(b) Step 120
(c) Step 160
(d) Step 250
Example 4 (similar indoor uniaxial compression test is simulated for the same inclined jointed rock masses). The specimen size is the same as in the example 1. Physical parameters of the calculation model are shown in Table 1 above; and the calculation process is shown in Figure 9.
Figure 15 is the indoor uniaxial compression experimental result, and it is shown that the crack extended from both ends of the joint to the top and bottom of the specimen until the specimen overall rupture instability in the compression process. Simulated indoor specimens’ crack propagation under uniaxial compressive and fracture rock mass instability process is shown in Figure 16. Figure 17 is the damage diagram schematic after the the rupture instability of the model. Comparing between Figures 15 and 17, the calculation result is consistent with the experimental results. As the comparison chart of the sample experimental and computing results, Figure 18 shows that curves change similarly and the peak strength is close. This illustrates the correctness and feasibility of the proposed numerical method for crack propagation of jointed rock mass.
(a) Step 105
(b) Step 145
(c) Step 180
(d) Step 275
5. Conclusion
This paper, based on the three-dimensional manifold cover, constituted the three-dimensional manifold cover function and proposed a new discrete element numerical method based on manifold cover, manifold particle discrete, combined with the particle contact model. This method can be applied to the entire process of analyzing rock-like brittle material failure. It also verifies the accuracy and feasibility of the proposed numerical method for crack propagation and rock failure analysis through numerical examples.
MPD method inherited structural response process of discrete particles and used the modeling and analysis environment of the YADE. It can be performed to analyze stress and strain of each polyhedron unit, which reflects the heterogeneity of the material characteristics. By building three-dimensional model, the static fracture, the dynamic fracture, and instability failure whole process of the material can be analyzed from the perspective of macroscopic and microscopic. This method has good research and application prospects.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank the support of the National Basic Research Program of China (Grant no. 2010CB732002), the National Natural Science Foundation of China (Grant no. 5117909851379113), the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant no. 20120131110031), and the Program for New Century Excellent Talents in University of Ministry of Education of China (Grant no. NCET-12-2009).