Estimation of fracture systems parameters in rock massif by the finite element method

. Results of analytical, laboratory and mine studies of parameters of disintegration zones in structurally-heterogeneous rocks and fracture systems around of the deep mine roadways are presented. A mathematical model is proposed, which is realized with the help of the procedures of simulation modeling by finite element method. Based on the calculation of possible directions of the rock shearing (the “rock shear sites”) and breakage of bonds between the elements, the model allows determining orientation of dominant master fracture system development with taking into account natural structural defects in the rock massif. A new algorithm was created for determining master cracks and boundaries of disintegration zones with the layered rock, which differs by determining directions of rock shear sites in the elastoplastic problem considering residual strength and main structural defects of the structurally-heterogeneous rock massif. Regularities of organization of unidirectional rock shear sites and their distribution under the influence of mining operations in zones with inelastic deformations were established.


Introduction
Experience of preparatory roadway operation in the Ukrainian mines over the last ten years has shown that 30 % of them are re-supported. It is caused, first of all, by deepening and increased pace of mining operations and intensive manifestations of rock pressure. One of the ways to solve the problem of providing stability and proper performance of preparatory roadways is to elaborate new methods and facilities for their supporting with maximum use of load-bearing capacity of the rock massif, which is determined by computational or experimental methods. Theoretical methods for solving such problems give mathematically incorrect results, since strain increase occurred during the rock loosing is accompanied by stress decrease in absolute magnitude, while assignment of structural weakening coefficients for theoretical calculations is very problematic. In such conditions, quantitative analysis of stresses and strains can be correctly carried out with the help of situational simulation geomechanical models based on dependencies of continuum mechanics and the theory of limit states.
When solving complex problems of analyzing the rocks massif stress-strain state, the combined procedures of finite element method and initial stress method are used [1][2][3][4][5]. These methods have quickly occupied certain position in the practice of complex engineering calculations. However, despite the great number of works devoted to the computational methods for modeling fracturing processes, there is still no clarity in choosing mathematical model adequate to the conditions of structurally inhomogeneous rock deformation, methods of its implementation and ways of its adaptation to real conditions. Methods for analyzing geomechanical processes of such complex systems require further improvement, since there are no reliable geomechanical models and modern specialized computational programs, which take into account specific deformation of the fractured structurally-heterogeneous rocks. Despite significant progress in rock massif modeling, it is still quite uncertain of how to determine zoning and intensity of fracture systems, boundaries of zones with layered rock disintegration, spatial orientation and location of technogeneous cracks. However, it is the presence of developed fracture system that is the dominant factor for mining operations, which negatively affects stability of mine roadways and gas-content in the mining district of the mines. Therefore, to improve safety and efficiency of underground mining operations, the live task is to develop new methods, algorithms and modern specialized programs, which will determine parameters of zoning and spatial orientation of fracture system development in the rock massifs.

Methods
Methodology for organization of computational processes, methods and algorithms for determining parameters of zoning and spatial orientation of fracture system development in the rock massif; methods for constructing the structures and software models of information systems for developing basic elements for software and information technology for geomechanical process calculation by the finite element method; mathematical modeling of the rock-massif stress-strain state by the finite element method; analysis and assessment of parameters of fracture systems in the structurally-heterogeneous rock massif performed with the help of information system "GEO-RS © " developed under the guidance and with direct participation of the author [6][7][8].

Development of the method for determining master crack orientation by using procedures of the finite element method
The basic procedure of the finite element method is well known: it is used for determining theoretical elastic parameters of the rock massif stress-strain state. However, factual characteristics of the rock deformation demonstrate various inelastic properties. Some rocks, such as sandstone, are more brittle, while others (mudstone, siltstone) feature certain degree of plasticity. In the process of deformation, as it is known, rocks manifest both plastic and brittle properties in different proportions. In contrast to the brittle properties, plastic properties of the rocks are of different nature, which is explained by "viscous" bonds between dislocations appeared at the rock shearing. Therefore, factual deformation processes occurring in different types of rocks most closely correspond to the elastoplastic model of the medium: this fact was confirmed by many analytical, laboratory, and field studies [3,9]. Theoretical issues of and methods for calculating the stress-strain state of medium in the area with inelastic deformations, which allow obtaining a final elastoplastic solution for the geomechanical problem under the given conditions, are disclosed in [3-6, 8, 10].
Inelastic deformation of rocks is a result of micro shears inside the rock polycrystalline structure. Therefore, the criterion of exceedance for external stresses of bonds under the condition of non-uniform compression can be described quite accurately by the Mohr-Coulomb failure criterion [3,9]: where τ is shear stress, Pa; σ n is stress normal to the rock shear site, Pa; φ is internal friction angle, degree; τ 0 is shear strength (cohesion), Pa.
The tensile strength was limited by the criterion τ 0 /5. If one of the limiting conditions is fulfilled, then the rock massif stress-strain state and inelastic deformation zones are determined by combination of finite element method and iterative procedure of the initial stress method (a set of calculations are performed with a constant system stiffness matrix [3][4][5]). The procedure of the simulation algorithm was implemented in the computer complex "GEO-RS © " [7,11].
The factual maximum σ 1 and minimum σ 3 principal stresses, which the real rock would have during the inelastic deformation, are determined from the linear strength reduction function [3]: where ε 1 is current value of the model element deformation obtained during the calculation; σ 1 lim , σ 1 res are limit and residual maximum principal stresses, respectively, Pa; ε 1 lim is limit maximum principal deformations; k is coefficient, which characterizes degree of the rock brittleness [11]; μ is the Poisson's ratio; ctgξ is parameter, which determines law of plastic flow.
Fracture is the most important property of rocks, which determines their structure and spatial heterogeneity and affects their strength and stability. Natural (primary) fracturing is interpreted as a combination of discontinuities in the continuous medium, with no or very insignificant displacements, and their size exceeds, by an order and more, interatomic distances in the crystal lattice [12]. Technogeneous fracturing occurs only during the construction and operation of underground mine workings and is characterized by shifts (often dynamic) along the contacts of master cracks, the cracks opening and change of their orientation in space. Even approximate forecasting of technogeneous fracturing is of considerable scientific and practical interest, since, unlike primary fracturing, development of technogeneous cracks depends not only on random natural mining-and-geological factors, but also on the controlled mining-and-technical parameters of the mining operations, including, first and foremost, design and power parameters of supports in the preparatory roadways, their orientation in geospace, parameters of methods for roadway driving and protecting, as well as systems designed for rock pressure control in the longwalls. Having available information about parameters of fracturing, it is possible, by means of technical facilities, to control resistance of the rock massif and loads on the roadway supports.
Cracks between the crystals and in the intercrystalline cement characterize microfracturing of the rocks, which, in its turn, decreases the rock total strength. Therefore, natural fracturing is usually considered as initial one and is taken into account when formulating a problem by applying structure weakening coefficients, which are determined by data of the weakened rock sample testing. This approach allows taking into account the overall degradation of the rock strength properties in the design scheme elements, which simulate the fractured massif. Since technogeneous cracks are formed during mining operations, they can be considered as additional ones, and their parameters depend on changes in the stress fields and deformations.
The cracks are formed only as a result of brittle rock breaking. Therefore, when isotropic rock with no pronounced initial fracturing is broken, the master fracture systems develop either in directions of or near to perpendicular to the maximum normal tensile stresses in case of tensile strain (tensile cracks) or in directions of maximum shear stresses (shear lines) in case of non-uniform compression (splitting cracks) conditions. The tensile cracks are short, have uneven rough surfaces, most often are open and widespread over the mined-out space of the longwall. The splitting cracks are smoother, are often closed and stretched for long distances in the rock massif and are oriented in a certain plane.
In this paper, an algorithm is implemented, which allows determining shear lines for an elastoplastic problem, as well as analyzing the entire area under the study with considering main structural defects existing in the real rock massif. The procedure is that as a result of solving the elastoplastic problem with the help of analysis block of the algorithm for calculating parameters of the rock massif stress-strain state, directions of vectors of the maximum principal stresses for each element of the design scheme are determined by the values of the previously calculated stress tensor components σ 1 , σ 3 and τ xy : where α is the angle of maximum principal stress direction towards to the coordinate axe x, deg; τ xy is shear stress, Pa.
At the same time, it is necessary to mention one of the major factors affecting the rock massif breakage: fracture systems with the same or similar orientation are developed with the greatest intensity and are organized into the master crack systems. Rather often, these cracks cause such abnormal phenomena as: decrease in stability and sharp movements of the roof and floor rocks, increase of coal and enclosing rock slip, layering of the rock massif edges, sudden roof caving, dynamic loads on the supports, etc.
Due to the fact that mechanism of plastic deformation is associated with the shears in the rock sites under condition of inelastic state, then, in zone with inelastic deformations, condition of limiting state [3] occurs in two sites, which are inclined to the maximum main compressive stress at angles of: As a result of calculations, an array [α] is formed with angles of maximum principal stress inclination towards to the coordinate axe x. Values of the internal friction angles are extracted from the model program class containing physical-mechanical properties, and then, by calculating by formula (4), a two-dimensional array is formed with angles of shearline vector inclination, which are defined for the center of gravity of each element (Fig. 1).
Since natural structural defects and structural inhomogeneity are typical for all rocks, the model takes into account the fact that fracture systems appeared during mining operations and rock breakage areas formed under the action of deviator loading occur mainly in the existing systems of natural cracks and structural defects in the rock mass. Such a development of cracks occurs when ratio of minimum component σ 3 of the stress tensor to the maximum σ 1 is less than 0.4, because at greater values, effect of structural defects is almost zero [3]. Consequently, with increasing σ 3 , direction of the dominant crack development will be determined mainly by the shear occurred in the monolith rather than towards direction of the rock layering or structural defect. In this connection, new parameters were introduced into the analysis block of the software algorithm: angle of direction of prevailing structural defects and angle θ (Fig. 1), which characterizes rate of coincidence between direction of structural defect and direction of maximum shear force. This coincidence is shown, with a given accuracy (in degrees), by the bold lines in visualization unit of the computational complex. A new service function for calculating directions of action of maximum principal and shear stresses was developed and introduced into the software package, which made it possible to interactively determine angles of structural defect inclination in two arbitrary directions and, thus, specify two dominant natural fracturing systems in the model. Since, as a result of calculations by finite element method, a huge data arrays should be processed, a new module for calculating and visualizing the shear lines was developed, tested as part of the research, introduced into the "GEO-RS © " and used for automatizing graphical processing of computational results.

Verification of the model and estimation of parameters of master fracture systems
Comparison of the calculated areas of shear along the radial shear surfaces with the known solution [13]; (Fig. 2, a) was carried out on the example of a punch impacting on a coherent weighable region under the conditions of maximum critical load. The area with dimensions of 40×20 m was investigated. The following rock properties were assumed: elastic modulus 4.42·10 4 MPa, Poisson's ratio 0.3, volume weight 2.6 t/m 3 , shear strength 9 MPa, internal friction angle 33. The punch was presented as a distributed load applied to the upper edge of the design scheme. Critical load for these conditions was determined by gradually increasing values of the forces simulating the punch action up to the value at the moment when one of the elements transited to the inelastic state. This stress-strain state was fixed and shown in Figure 2, b. A map of the rock shear sites, which were determined for each element of the design scheme, was aligned with the diagram of maximum principal stress. It is known that three conditional zones are formed under the action of the punch (Fig. 2, a). Immediately under the punch, a compression zone is created in the form of an elastic "core" (zone 1), which "pushes" the rock in different directions. At the boundaries of this zone, the conditions of limiting state are created and breakage sites occur, which is confirmed by a numerical experiment with a higher load applied on the punch. Zone 2 is characterized by the shearing force increase in different directions from the compression core. Theoretically, zone 2 (the region of shearing along the radial shear surfaces) has only conditional parameters (type, shape) as its real geometric dimensions and location of the boundary line are not unambiguously determined. Nevertheless, as a result of calculation by the finite element method, the boundary line of the zone 2 is unambiguously determined as the border between the ordered and disordered shear lines (the border is shown by a small dotted line). As it is seen in Figure 2, b, its depth is ≈6.0 m. In contrast to zone 2, shear lines in zone 3 are located randomly, and deformations are characterized by displacements of some rocks in the massif towards the open space (the border is shown by the dotted lines). The model allowed determining real geometric parameters of this zone (its depth, width and configuration). As can be seen from the comparison between the Figures 2, a and  2, b, type of zones, their location and configuration, shear lines and boundary lines between the sectors with different types of deformations correspond to theoretical concepts of the rocks stress-strain state under conditions of maximum load.
The quantitative and qualitative assessment of the proposed method was studied in the model on the example of the controlled breakage of the sandstone sample and by its comparison with the results of press tests under uniaxial compression. The physical model of the sample with dimensions of 60×60 mm was digitized in the form of triangular grid with 3200 finite elements and 1682 nodes. That is, the stress value was obtained for each area of 0.75 mm. Sandstone properties were as follows: elastic modulus 4.42·10 4 MPa, Poisson's ratio 0.28, volume weight 2.55 t/m 3 , shear strength 25 MPa, internal friction angle 37. The strength characteristics established in the model corresponded to the ultimate strength for uniaxial compression of 100 MPa. The load was applied vertically in several stages, each of 2 MPa. The upper sector of the design scheme was set movable in vertical direction, and the lower one was fixed in two directions. That is, the idealized experimental conditions were simulated without the sample's upper and lower edges shearing in horizontal direction.
The sample stress-strain state was determined both for elastic and inelastic deformation stages. It is established that zones with unidirectional rock shear sites (are shown in Figure 3 Analysis of rock shearing directions in the inelastic deformation zone shows that under the given conditions of sample fixing, its edges are not deformed from the top and from the bottom simultaneously. First, its bottom part is layered, and later -its top part; this fact was confirmed by uniaxial press tests used for studying nature of the breakage (Fig. 3, d). The cone-shaped blocks in the edges, top and bottom remained undisturbed.
As it is shown in Figure 3, d, the total deformation of the real sample also consists of deformations caused by the crack development due to the action of shear stresses along the shear surfaces located at some angle to the direction of the applied loads. It is obvious that directions of dominant fracture, which is formed in zones with inelastic deformations of the real rock massif, are caused by micro and macro defects of the rock (size of mineral grains, fracturing, layering), and depend on the type of stress state.
Simulation results and results of laboratory experiment were compared in terms of forms of the sample breakage and values of longitudinal and transverse deformations. The correlation analysis established a good association by the values of the data sets: correlation coefficient for longitudinal strain was 0.94, and for transverse strain -0.96. It is established that differences in deformations between the simulated results and results of laboratory experiments are mainly explained by deliberate idealization of the rock microstructure.
In the process of the model verification, efficiency of the developed method of mathematical modeling was checked on the example of numerical experiments on studying the punch impact on the coherent weighable area under critical load and the controlled breakage of sandstone sample as an element of the massif under the effect of simplified uniaxial compression. During the model verification, more than 150 computational experiments were analyzed and processed, which showed that shear lines were in harmony with theoretical concepts and principles of the rock mechanics. At the second stage, we will consider result of calculation of a more complex multicomponent model on the example of the rock massif sector containing air roadway and the mined-out space.
In Figure 4, diagrams of orientation of maximum principal stresses and fracturing systems are shown, which are obtained by calculation of the multicomponent simulation geomechanical model for conditions of the 18 th western longwall, seam m 3 , horizon 1340 m. Reliability of the results of numerical analysis and their compliance with real deformation processes in the rock massif were provided through the schematization and introduction into the geomechanical model of specific features inherent to the rock massif under study and its individual parameters determined by experimental studies. Relationship between the values obtained by the mine observations and data of numerical modeling was established on the basis of mine research data on development of cracks and their location relatively to the roadway, layering of the massif, displacements and convergence of the roadway roof and floor, and control of boundary of inelastic deformation zone.
Analysis of the rock shearing directions shows that the sites are ordered in zones with inelastic deformations, where layering blocks are formed due to the dislocation shearing. Before some rocks are broken, the concentrated shear forces are re-oriented contributing to occurrence of directed master cracks. Systems of cracks in the main roof are developed in close to horizontal and close to vertical directions. Above the mined-out space, one of the fracture systems is ordered. As a result of the roof deformation, technogeneous cracks form horizontal master cracks and layering planes, which are parallel to the plane of the minedout coal seam. In zone of stope, the shear lines are parallel to the exposed surface, and horizontal layering is formed in the stope floor and roof as a result of the shear. Measurements of the rock layering in three deep benchmark stations of the roadway roof confirmed the nature of location of the rock layering zones.
Most of the rock pressure manifestations are accompanied by gas release into the mine roadways. At the same time, in terms of impact on the massif destruction and gas release, quantity and direction of the dominant cracks are main parameters of the fracture system. For example, the most dangerous are cracks developed in parallel to the exposed surface of the roadway. No less important is development of crack dip towards the drift driving direction. Knowing orientation of the systems of endogenous, tectonic and technogeneous cracks in the geospace, it is possible, by technical means, to increase resistance of the rock massif to arbitrary deformations and to control loads on the roadway supports. is zone of elastic deformation; is directions of vectors of maximum principal stresses; is orientation of two shear sites in each element of the design scheme; is shear lines, which coincide with the direction of the rock massif layering (± 3); is direction of the stope driving.
Thus, it is for the first time when a mathematical model is proposed, which allows determining orientation of the dominant master fracture system development with taking into account natural structural defects in real rock massif on the basis of calculation of possible directions of rock shear sites and breakage of bonds between the model elements with the help of simulation modeling by finite element methods. Despite a certain idealization, which is typical for any modeling process, the designed technology of computer analysis features essential advantages as it opens up new broad possibilities for estimating the stress-strain state of the fractured structurally-heterogeneous rock massif. In particular, new opportunities are created for solving a number of topical scientific and applied problems for determining such parameters as: zonality, development and intensity of fracture systems in the rock massif; boundaries of zones with layered rock disintegration; orientation of dominant master fracture systems. The model allows studying complex geomechanical processes in more details and correctly determining technical and technological parameters for technology of mining operations.
Introduction into production of recommendations on complex monitoring of the rock pressure manifestations and improving stability of preparatory roadways developed on the basis of validated force characteristics of the "support-rock massif" system with considering tectonic activity, inelastic deformations and technogeneous fracturing of structurally-heterogeneous rocks under the impact of winning operations allows getting the factual economic effect and proving efficiency and effectiveness of the obtained scientific results.

Conclusions
On the basis of the findings, the following scientific results are obtained: -a mathematical model is created, which allows, on the basis of calculation of possible directions of the rock shear sites and breakage of bonds between the model element with the help of procedures of simulation modeling by the finite element method, to determine orientation of dominant master fracture system development with taking into account natural structural defects in the rock massif; -a new object-oriented algorithm is developed for determining master fracture systems and boundaries of zones with layered rock disintegration; the algorithm differs by determining directions of the rock shear sites in the elastoplastic problem with taking into account residual strength and main structural defects of the structurally-heterogeneous rock massif; -regularities of organization of unidirectional rock shear sites with contouring of temporarily stable structures in zones with inelastic deformation and their distribution under the impact of mining operations in zones with inelastic deformations are established.
Practical importance of the work lies in the fact that, on the basis of the conducted research, a new technology was created for computer modeling of fracture systems and, in particular, a new version of the "GEO-RS © " software complex was developed for determining the stress-strain state of the rock massif and fracture systems, which was tested, in laboratory and field conditions, by way of comparative analysis of analytical inelastic simulation solutions in the mode of gradual deformation and real destruction of rocks.
The conducted studies disclose new possibilities for analyzing the stress-strain state of the rock massif allowing to predict the most probable directions of crack development, and, therefore, to study complex geomechanical processes in more details.