Surrounding rock mass cracking Fractal analysis method established based on particle flow discrete elements

Abstract

The invention discloses a surrounding rock mass cracking Fractal analysis method established based on particle flow discrete elements, and belongs to the technical field of rock mass fracture analysis. Constructing a surrounding rock body model according to the physical and mechanical parameters and the structural characteristic parameters of the target surrounding rock body; carrying out a uniaxial simulation test on the surrounding rock model by using particle flow discrete element software, and carrying out benchmarking on a simulation result and an indoor uniaxial compression test result of a target surrounding rock body to generate an initial numerical model; correcting the initial numerical model according to a comparison result of the multi-axis test and the multi-axis simulation test, and establishing a roadway surrounding rock data model under the influence of different stress numerical values; actual engineering load data are obtained, load application is given to the roadway surrounding rock data model for rock mass cracking simulation, crack distribution in the simulation process is obtained, expansion trend analysis is conducted on the crack distribution through a fractal dimension analysis method, and a rock mass cracking analysis result is generated. According to the invention, reliable safety early warning information can be provided for tunnel and tunnel engineering personnel.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of rock mass fracture analysis, and particularly relates to a surrounding rock mass cracking Fractal analysis method based on a particle flow discrete element. BACKGROUND

[0002] Under the background of the country's vigorous promotion of scientific and technological innovation and digital transformation, especially in the critical period of energy transformation and geological disaster prevention, it is necessary to explore and analyze the changes of underground rock mass.

[0003] The mechanical properties of surrounding rock mass mainly depend on its structural characteristics, the characteristics of structural rock mass and the properties of structural surface, and the existence of cracks will further affect its mechanical behavior and stability. In addition, the initial stress field of surrounding rock is the fundamental force of deformation and failure of surrounding rock, which directly affects the stability of surrounding rock. Although a large number of researchers have carried out a lot of research in this field, and have used different theories and methods to analyze the stress characteristics and cracking conditions of surrounding rock mass, due to the complexity of mechanical properties, the difficulty of integrating the stress state of surrounding rock mass, the difficulty of simulating the data of actual stress, and the insufficient consideration of the crack factor, etc. The strength detection system and even the safety warning system are still facing difficulties in development. SUMMARY

[0004] Therefore, the purpose of the present application is to provide a surrounding rock mass cracking Fractal analysis method based on a particle flow discrete element, which is used to numerically simulate the behavior of rock mass under multi-axial compression conditions by using the discrete element method, so as to obtain the characteristics of structure degradation and the trend of crack change of surrounding rock mass under different stress values, and to use the fractal dimension analysis method to study the crack distribution and describe the extension trend of the crack, so as to improve the accuracy and reliability of the existing surrounding rock mass strength detection and safety warning system, and to provide safety warning for tunnel and tunnel engineering personnel.

[0005] In order to achieve the above purpose, the present application provides the following technical scheme: A surrounding rock mass cracking Fractal analysis method based on a particle flow discrete element, comprising: constructing a surrounding rock mass model according to the physical and mechanical parameters and structural characteristic parameters of the target surrounding rock mass; performing a virtual uniaxial compression simulation test on the surrounding rock model by using a particle flow discrete element software, and comparing the simulation results with the indoor uniaxial compression test results of the target surrounding rock mass to generate an initial numerical model; performing a multi-axial compression test on the target surrounding rock mass, and performing a multi-axial simulation compression test on the initial numerical model, comparing the test results with the simulation test results, and correcting the initial numerical model according to the comparison results to establish a data model of surrounding rock under the influence of different stress values; The actual engineering load data is acquired, the data model of the roadway surrounding rock is loaded to perform rock mass cracking simulation, crack distribution in the simulation process is acquired, and a rock mass cracking analysis result is generated by using a fractal dimension analysis method to analyze the extension trend.

[0006] Further, by field investigation, a plurality of surrounding rock blocks meeting preset test conditions are acquired as target surrounding rock bodies for indoor tests. And the physical and mechanical parameters and structure characteristic parameters of the target surrounding rock bodies are collected. The physical and mechanical parameters include particle density, elastic modulus, Poisson's ratio, internal friction angle and cohesion. The structure characteristic parameters include particle gradation, porosity and crack distribution.

[0007] Further, the generation process of the initial numerical model includes: Step 1, ball constraint making; Step 2, generation of the particle material system; Step 3, uniformity and good contact of the model are realized by using a servo mechanism; Step 4, contact coefficient assignment; Step 5, compression process control, stress-strain curve acquisition and analysis; Step 6, if the stress-strain curve trend development in the analysis result of step 5 is consistent with the uniaxial compression test loading trend of the target surrounding rock body, it is determined that the current numerical simulation model is approximately matched with the pattern, and the initial numerical model is obtained.

[0008] Further, the multi-axial simulation compression test includes rigid biaxial simulation compression test, flexible biaxial simulation compression test and triaxial simulation compression test.

[0009] Further, the triaxial simulation compression test includes: Step A, sample preparation: generate a uniform particle assembly with a specified porosity to simulate the sample preparation process in the actual test; Step B, pre-pressing: apply an initial pressure on the sample to ensure that the internal stress of each particle is zero by using Newton's second law, so as to realize the densification of the sample, reduce the porosity inside the sample, and ensure that the sample meets the initial state of the actual material; Step C, cementing: introduce cementing force between the particles to simulate the bonding effect between the particles, so as to improve the overall strength and stability of the sample and ensure that the model can withstand the specified stress during the loading process; Step D, confining pressure: apply a uniform pressure around the sample to simulate the lateral pressure condition in the triaxial compression test; Step E, loading: apply an axial force on the sample to simulate the compression process in the triaxial compression test; Step F, data acquisition: acquire stress-strain data in the compression process to obtain simulation test results.

[0010] Further, in the indoor compression test of the target surrounding rock mass, the method for acquiring the crack distribution after compression comprises: acquiring a theoretical crack distribution range of the target surrounding rock mass after applying a corresponding load, denoted as a prior crack area, and generating a corresponding binary reference mask;wherein the number of prior crack areas is not less than one, and the pixel value of the prior crack area in the binary reference mask is a first value, and the pixel value of the non-prior crack area is a second value; According to the pre-trained convolutional neural network model, the original image of the compressed target surrounding rock mass is recognized to obtain a plurality of crack information, and the pixel point coordinate set of the crack edge of each crack information is extracted to obtain a to-be-verified edge pixel set;wherein the number of to-be-verified edge pixel sets is the same as the number of crack information; Spatially match any to-be-verified edge pixel set with all prior crack areas, and count the number of first pixel points in each prior crack area in the to-be-verified edge pixel set after spatial matching, and calculate the ratio of each first pixel point number to the total number of pixel points in the to-be-verified edge pixel set, denoted as the edge pixel overlap degree; If all edge pixel overlaps of any to-be-verified edge pixel set are lower than the preset overlap threshold, the to-be-verified edge pixel set corresponding to the crack information is determined as invalid crack and is removed, otherwise it is denoted as valid crack; Acquire all valid cracks to generate the crack distribution result of this compression test.

[0011] Further, a surrounding rock mass crack fractal analysis method based on particle flow discrete element comprises: for the first pixel points falling in the prior crack area, according to the Euclidean distance between the first pixel points and the center of the prior crack area, a weight is assigned, the weight is positively correlated with the distance, and the ratio of the number of weighted pixel points to the total number of pixel points is calculated to obtain a new edge pixel overlap.

[0012] The beneficial effects of the present application are: The discrete element method can accurately describe the interaction and mechanical behavior between particles in the rock mass, and is particularly suitable for studying the mechanical properties of fractured rock mass and the failure mode under complex stress state; in addition, the fractal dimension is an effective tool for quantifying the complexity and extension trend of crack distribution; therefore, the present application provides a surrounding rock mass cracking fractal analysis method based on particle flow discrete element, which numerically simulates the behavior of the rock mass under multi-axial compression conditions by the discrete element method, so as to obtain the characteristics of the structure degradation of the surrounding rock mass of the roadway under different stress values and the trend of the crack change, and the fractal dimension analysis method is used to study the crack distribution and describe the extension trend of the crack, so as to improve the accuracy and reliability of the existing surrounding rock mass strength detection and safety warning system, and provide safety warning for tunnel and tunnel engineering personnel.

[0013] Other advantages, objects, and features of the present application will be apparent to those skilled in the art from the following specification and drawings, and it is intended to cover any and all adaptations or variations of the application.

[0014] The technical solutions of the present application will be described in further detail below with the help of the drawings and examples. BRIEF DESCRIPTION OF DRAWINGS

[0015] The accompanying drawings are included to provide a further understanding of the present application, and constitute a part of the specification, illustrate the present application and explain the principles of the present application, and do not constitute a limitation of the present application. In the drawings: Figure 1 The method flow chart of the surrounding rock mass cracking fractal analysis method based on particle flow discrete element in the embodiment of the present application; Figure 2 The stress distribution design drawing of the circular cross-section roadway in the surrounding rock mass cracking fractal analysis method based on particle flow discrete element in the embodiment of the present application; Figure 3 The tangential stress nephogram of the splitting failure of the surrounding rock of the circular cavern in the surrounding rock mass cracking fractal analysis method based on particle flow discrete element in the embodiment of the present application; Figure 4 The ideal and homogeneous simulation object schematic diagram obtained in the surrounding rock mass cracking fractal analysis method based on particle flow discrete element in the embodiment of the present application; Figure 5 The force decomposition schematic diagram on the contact surface in the pre-pressing step of the surrounding rock mass cracking fractal analysis method based on particle flow discrete element in the embodiment of the present application. Detailed Implementation

[0016] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0017] like Figure 1 As shown, this invention proposes a fractal analysis method for surrounding rock mass fracturing based on particle flow discrete element method, including: S101. Construct a surrounding rock mass model based on the physical and mechanical parameters and structural characteristics of the target surrounding rock mass; S102. A virtual uniaxial compression simulation test was conducted on the surrounding rock model using particle flow discrete element software, and the simulation results were compared with the indoor uniaxial compression test results of the target surrounding rock mass to generate an initial numerical model. S103. Conduct multiaxial compression tests on the target surrounding rock mass and multiaxial simulated compression tests on the initial numerical model, and compare the test results with the simulation test results. S104. Based on the comparison results, revise the initial numerical model and establish a data model of the surrounding rock of the roadway under the influence of different stress values. S105. Obtain actual engineering load data and apply loads to the roadway surrounding rock data model to simulate rock mass cracking. S106. Obtain the crack distribution during the simulation process, and use the fractal dimension analysis method to analyze its propagation trend, generating rock mass cracking analysis results; The working principle of the above technical solution is as follows: This invention is proposed based on the defects of the existing technology. The defects of the existing technology in the analysis of cracking of surrounding rock mass will be explained in detail below. First, we will take the stress distribution of the tunnel wall in the surrounding rock mass as an example. In principle, the stress distribution around a tunnel of any shape can be studied using analytical methods in homogeneous and isotropic elastic rock masses. However, since the derivation and calculation are often very cumbersome, a computer must be used. The following content uses a simple example of a circular cross-section tunnel to illustrate the principle of the analytical method of elastic theory: If a tunnel is long enough, the stress distribution and deformation will be the same on any cross-section perpendicular to its axis and away from its end. To obtain this distribution, it is sufficient to simplify the above to a single cross-section, which is to solve the so-called plane problem.

[0018] In planar problems, an orthogonal rectangular coordinate system is used: the horizontal axis is the x-axis, and the vertical axis is the y-axis, with the origin inside the cross-sectional axis. Furthermore, polar coordinates can also be used for circular tunnels. The method involves setting the origin of the coordinate system at the center of the tunnel cross-section, such as... Figure 2 As shown.

[0019] Figure 2 For any unit cell, the three unknown stress components are , and ; the three strain components are , and , all of which must satisfy the equilibrium equations: and the compatibility equations: and the two generalized Hooke's law equations: As boundary conditions, the stresses far away from the roadway are taken to be equal to the in-situ stresses at the center of the roadway before the roadway is excavated. That is: At the boundary of the roadway, : To solve the above set of equations, that is, to find the stress and strain components, one can use a method from elasticity theory. The three stress components are expressed in terms of an unknown stress function F: Substituting these expressions into the equilibrium equations, one can show that for any form of F, the following identity is satisfied: The outer boundary of the roadway influence zone is assumed to be a circle of radius , and the radial and tangential stresses and are taken to be equal to the in-situ stresses.

[0020] For the in-situ stresses, one can transform from the x, y coordinate system to the r, coordinate system. The angle between the direction of the principal stress and the direction of the in-situ stress is equal to According to the basic formulas of elasticity theory, the stresses can be expressed in the following form: To simplify the calculations, one can express , in terms of and as follows: The following can be derived The peripheral stress at the time: The outer boundary stress Composed of two components.

[0021] The first component, the radial stress, does not change with the angle and is uniformly distributed, expressed as: The second component is expressed as: In addition, the outer boundary also acts tangential stress, which can be determined by the second formula in the above description of the peripheral stress corresponding formula: It is worth mentioning that, in order to reduce the repetition, this paper omits the detailed derivation process, and aims to guide the technical personnel to better understand the technical content of the scheme through the derivation results; After the above derivation, the elastic mechanics solution of the surrounding rock mass under pressure is studied. A circular ring or cylinder is provided, the inner radius is a, the outer radius is b, the inner pressure is , the outer pressure is , and the stress component and external conditions are as follows: Stress component: Boundary conditions: Bring them into the deformation compatibility equation and the balance differential equation (i.e. partial differential equation) recorded above, and get the solution of Lame as follows: If only the inner pressure acts, then The solution is simplified as: When the outer radius of the circular ring or cylinder is infinite (i.e. b→ ), it can be regarded as an infinite thin plate with a circular hole, and the above solution is: If there is no inner pressure, i.e. The solution is simplified as: In addition to the ideal elastic theory calculation solution, scholars have also carried out numerical simulation research, such as Figure 3The tangential stress nephogram of the circular cavern surrounding rock splitting failure shows that the tangential stress concentration effect at the horizontal hole wall is remarkable, and needs to be focused on.

[0022] The above calculation can obviously show that the numerical simulation research does not show the corresponding crack occurrence and expansion mechanism, therefore, the sampling of the rock mass at the hole wall for crack analysis research proposed in the application is feasible and pioneering.

[0023] After the defects of the prior art are clear, the application proposes a surrounding rock mass fracturing Fractal analysis method based on particle flow discrete element, for solving the problems in the prior art. Specifically, the technical scheme of the application mainly uses a particle flow discrete element program (hereinafter referred to as PFC software) for corresponding digital simulation research. The application of the PFC software aims to construct a numerical model under the simplified laboratory test conditions according to the site, and combine the mechanical parameters of the surrounding rock mass obtained by the test (also can combine the existing literature) to correct the numerical model, and establish the roadway surrounding rock data model under different stress numerical effects, provide data support for the strength monitoring analysis of the surrounding rock mass, and obtain more in-depth information of the loaded body which cannot be obtained by relying on mechanical loading alone. Then according to the construction scheme designed after the site survey, the actual engineering load data is extracted, and the load is applied to the roadway surrounding rock data model for rock mass fracturing simulation. In this simulation compression process, the crack distribution is obtained, which is binarized, the effective crack contour is retained, and the mathematical software is used for dimension calculation. The calculation method of fractal dimension is relatively rich, but considering that the rock mass material has high self-similarity, the box dimension method commonly used at present is preferably used for calculation. Since the box dimension method belongs to the prior art, and the content of the application is described above, those skilled in the art can clearly understand how to use the box dimension method to implement the crack analysis process, so the application will not be described again. The crack analysis method proposed by the method is integrated into the safety early warning system to improve the accuracy and reliability of the existing surrounding rock mass strength detection and safety early warning system.

[0024] The beneficial effects of the above technical scheme are: through the above technical scheme, the behavior of the rock mass under multi-axial compression conditions is numerically simulated by using the discrete element method, so as to obtain the characteristics of the structure deterioration of the roadway surrounding rock mass and the trend of the crack change under different stress values, and the fractal dimension analysis method is used to study the crack distribution and describe the crack expansion trend, so as to improve the accuracy and reliability of the existing surrounding rock mass strength detection and safety early warning system, and provide safety warning for tunnel and tunnel engineering personnel.

[0025] In one embodiment, a plurality of surrounding rock blocks meeting preset test conditions are obtained as target surrounding rock bodies for indoor tests through field investigation; And the physical and mechanical parameters and structural characteristic parameters of the target surrounding rock bodies are collected; The physical and mechanical parameters include particle density, elastic modulus, Poisson's ratio, internal friction angle and cohesion. The structural characteristic parameters include particle gradation, porosity and fracture distribution. In addition to meeting the requirements of basic tests, the test conditions should also include shape matching, existing crack matching, and physical and mechanical parameter matching and structural characteristic parameter matching, so as to ensure that the target surrounding rock bodies used before each compression test are basically in the same state, with an error of not more than 0.005, to improve the comparison effect of subsequent simulation tests and actual test results. Similarly, when constructing the surrounding rock body model, the above test conditions should also be added for corresponding construction, so as to generate a simulation model with high matching degree. The working principle and beneficial effects of the above technical solution are that the actual data of the target surrounding rock body are obtained, which is beneficial to providing reliable basic data for subsequent numerical simulation research.

[0026] In one embodiment, the generation process of the initial numerical model includes: Step 1, production of the constraint ball; Step 2, generation of the particle material system; Step 3, use of the servo mechanism to make the model uniform and well-contacted; Step 4, assignment of the contact coefficient; Step 5, control of the compression process, acquisition of the stress-strain curve for analysis; Step 6, if the analysis result of step 5 shows that the stress-strain curve trend is consistent with the loading trend of the target surrounding rock body in the uniaxial compression test, it is determined that the current numerical simulation model is approximately matched with the pattern, and the initial numerical model is obtained. The working principle of the above technical solution is that: in the PFC software, the boundary of the preset range is drawn by the boundary command, and the constraint particles of the specified particle size are generated in the boundary, and only axial movement is allowed (only when uniaxial compression simulation is performed, subsequent verification and correction are performed according to the specific axis to re-limit the constraint condition); after determining the constraint particles, the surrounding rock mass model established in advance is used to restore the mesostructure of the surrounding rock mass (including particle size distribution, porosity and other parameters), to provide a test sample for subsequent compression test, including generating a particle assembly, adjusting the initial distribution of particles and the like; after the above work is completed, the servo mechanism is used to make the model uniform and well contacted, and the servo mechanism is related to the calculation core of the particle flow numerical method and is crucial to the calculation of the model. When the particle assembly is initially generated, the existence of the servo mechanism avoids a large amount of overlap of the particles, releases the strain energy stored in the particles, and avoids the particles flying out of the wall. Therefore, the servo mechanism must be used to reduce the possible overlap between the particles and reduce the contact force between the particles as much as possible. When the numerical test simulation is performed, the servo control can help control the size and balance of the load. Since the servo mechanism is commonly used in the field, the specific principle of the servo mechanism is not described herein; after the particle model is determined, the contact model is selected according to the actual situation, and the corresponding contact coefficient is set to value the particle model, and the compression process is controlled according to the pre-determined compression test process, the stress-strain curve is obtained for analysis, and when the trend of the stress-strain curve development is consistent with the mechanical loading trend of the test piece, it is considered that the numerical simulation model is approximately matched with the test sample, and on the basis of the ideal simulation object, the actual rock mass test sample (i.e., the target surrounding rock mass described above) is adjusted to match the actual rock mass test sample (i.e., the target surrounding rock mass described above); Figure 4 As shown in the figure, the position distribution of the ball is shown when the sample is formed, the small balls are uniformly distributed, and the position change of the ball on the z-axis is indicated by different colors, and on the basis of the ideal simulation object, the actual rock mass test sample (i.e., the target surrounding rock mass described above) is adjusted to match the actual rock mass test sample (i.e., the target surrounding rock mass described above); The beneficial effects of the above technical solution are that: through the above technical solution, the numerical model is completed, and reliable data support is provided for subsequent numerical simulation.

[0027] In one embodiment, the multi-axial simulation compression test includes a rigid biaxial simulation compression test, a flexible biaxial simulation compression test and a triaxial simulation compression test; The working principle and beneficial effects of the above technical solution are that: in order to ensure the accuracy and reliability of the simulation results, the numerical model needs to be continuously adjusted, therefore, after the uniaxial simulation compression test, the multi-axial simulation compression test is added to improve the accuracy of the subsequent simulation test, and to provide a reliable test model for subsequent crack analysis.

[0028] In one embodiment, the triaxial simulation compression test includes: Step A, sample preparation: generate a uniform particle assembly with a specified porosity to simulate the sample preparation process in actual tests; Step B, pre-pressing: apply an initial pressure on the sample to ensure that the internal stress in each particle is zero using Newton's second law, to achieve compaction of the sample, reduce the porosity inside the sample, and ensure that the sample meets the initial state of the actual material; Step C, cementation: introduce cementation forces between particles to simulate the bonding effect between particles, to improve the overall strength and stability of the sample, and ensure that the model can withstand the specified stress during loading; Step D, confining pressure: apply a uniform pressure around the sample to simulate the lateral pressure conditions in triaxial compression tests; Step E, loading: apply an axial force on the sample to simulate the compression process in triaxial compression tests; Step F, data acquisition: acquire stress-strain data during compression to obtain simulation test results; The working principle of the above technical solution is: the multi-axial simulation compression test proposed in the present application includes a biaxial simulation compression test and a triaxial simulation compression test. Here, the triaxial simulation compression test is taken as an example to describe the basic idea of the multi-axial simulation compression test for the person skilled in the art, to help the person skilled in the art better understand the technical content recorded in the present technical solution; In addition to data acquisition, it mainly includes sample preparation, pre-pressing, cementation, confining pressure, and loading.

[0029] Sample preparation refers to generating a relatively uniform particle assembly with a specified porosity to simulate the sample preparation process in actual tests.

[0030] Pre-pressing refers to applying an initial pressure on the sample to ensure that the internal stress in each particle is zero using Newton's second law, which helps to achieve compaction of the sample, reduce the porosity inside the sample, and thus more closely approach the initial state of the actual material. The force decomposition on the contact surface during pre-pressing is shown in Figure 5 .

[0031] Cementation refers to introducing cementation forces between particles to simulate the bonding effect between particles. Cementation can improve the overall strength and stability of the sample, allowing it to withstand greater stress during loading.

[0032] Confining pressure refers to applying a uniform pressure around the sample to simulate the lateral pressure conditions in triaxial tests. The application of confining pressure is crucial for simulating the lateral pressure conditions in triaxial tests, as it directly affects the stress-strain behavior of the sample. In the following section, the strength deformation and failure conditions of rock mass under different confining pressures will be demonstrated in detail in the stage report.

[0033] Loading refers to exerting an axial force on the sample to simulate the compression process in triaxial test. During the loading process, the sample is subjected to axial stress, thereby deforming and failing.

[0034] Since the technical solution is mainly implemented in PFC software, part of the code information of PFC software implementing the simulation process is provided here to help those skilled in the art better understand the triaxial simulation process. First, the first to tenth lines of the sample are as follows: 1 new 2 def chicun_par 3 4 sample_rad=0.2 5 sample_hight=sample_rad*4 6 7 Keli_rdmin=0.006 8 Keli_rdmax=0.009 9 End 10 @chicun_par The functions of different lines (also known as Line, Line N, and Line N) are as follows: Line 1 represents creating a new command file; Line 2 defines a function for subsequent Line 10 to call and execute all commands in the process; Line 4 sets the sample radius to 0.2; Line 5 sets the sample height to 4 times the radius, i.e. 0.8; Line 7 sets the minimum radius of the particle to 0.006; Line 8 sets the maximum radius of the particle to 0.009; The main code information of the 14th to 38th lines of the sample is as follows: 14 domain extent [-sample_rad*1.5] [sample-rad*1.5][-sample_rad*1.5][sample_rad*1.5]... &15 [-sample_hight*0.5*1.5] [sample_hight*0.5*1.5] 17 [n=1.4] 18 Wall generate cylinder base 0 0 [-sample_hight*0.5*n] axis 0 01... &19 | | | | height {sample_hight*n] radius [sample_rad]capfalse false 21 wal1 generate plane position 0 0 [sample hight*0.5] dip 0 ddir 0 22 wal1 generate plane position 0 0 [-sample hight*0.5] dip 0 ddir 0 24 bal1 distribute radius [keli_rdmin][keli_rdmax]porosity 0.28... &25 range cylinder end1 0 0 [sample_hight*0.5-keli_rdmin]... &26 end2 0 0 [-sample_hight*0.5+keli_rdmin] radius [sample_rad-keli_rdmin] 27 28 cmat default model 1inear method deform emod 10e8 kratio 1.5property fric 0.5 30 bal1 attribute density 2.7e3 damp 0.7 31 cyc1e 2000 ca1m 50 33 solve arat 1e-4 35 bal1 delete range cylinder end1 0 0 [sample_hight*0.5]... &36 End2 0 0 [-sample_hight*0.5] radius [sample_rad]not 38 Save sample wherein Line 14 to Line 15 set the range of the simulation domain according to the sample radius and height to ensure that the simulation domain is large enough to accommodate the subsequently generated particles and wall; Line 17 sets the variable n for adjusting the size of the wall; Line 18 to Line 19 generates a cylindrical wall as the bottom boundary of the simulation according to the sample radius and height, with the center of the cylinder as the origin; Line 21 to Line 22 generates two plane walls as the side boundaries of the simulation; Line 24 to Line 26 distributes particles within the specified cylindrical region, satisfying that the particles are between the maximum radius and the minimum radius, and the porosity is 0.28; Line 28 sets the default contact model (cmat) as the linear deformation model, the elastic modulus (emod) as 10e8, the Poisson's ratio (kratio) as 1.5, and the friction coefficient (fric) as 0.5; Line 30 sets the density of all particles as 2700 kg / m³ and the damping coefficient as 0.7; Line 31 performs 2000 iterations of calculation, and after each iteration, 50 balance calculations are performed to stabilize the system; Line 33 performs solution calculation with an allowed relative error of 0.0001; Line 35 to Line 36 deletes particles located outside the cylindrical region to ensure that all particles are within the simulation domain; Line 38 saves the current simulation state as "sample".

[0035] The main code information of the first line to the 22nd line about pre-pressing is as follows: 1 restore sample 2 def wp_wall 3 wp_up=wall.find(2) 4 wp_down=wall.find(3) 5 wp_rr=wall.find(1) 7 loop foreach vt wall.vertexlist(wp_rr) 8 | vert_in_ce=vt 9 endloop 10 | 11 end 12 @wp_wall 14 [tzz=-1e6] 15 [trr=-1e6] 16 [sevro_fac=0.5] 18 [do_zservo=ture] 19 [do-rservo=true] 21 [sevro_frep=500] 22 [timestepNow=global.step-1] Line 1 restores the previously saved simulation state "sample"; Line 2 to Line 11 define functions for determining the various vertices of the wall; Line 12 calls the functions and executes all the commands defined therein; Line 14 sets the stress in the z direction to -1e6; Line 15 sets the stress in the radial (r direction) to -1e6; Line 16 sets the coefficient for servo control to 0.5; Line 18 enables servo control in the z direction; Line 19 enables servo control in the radial direction; Line 21 sets the frequency for servo control to 500Hz; Line 22 sets the current time step to the global step minus 1.

[0036] The main code information for the pre-press from Line 23 to Line 52 is as follows: 23 def servo_walls 24 computer_wallstress 26 if timestepNow<global.step then 27 get_gain(sevro_fac) 28 timestepNow+=sevro_freq 29 endif 31 if do_zservo=true then 32 z_vel=gz*(wszz-tzz) 33 wall.vel.z(wp_up)=-z_vel 34 wall.vel.z(wp_down)=z_vel 35 endif 37 if do_rservo=true then 38 r_vel_mag=(-1)*gr*(wsrr-trr) 39 loop foreach vt wall.vertexlist(wp_rr) 40 mag=math.sqrt(wall.vertex.pos.x(vt)^2+wall.vertex.pos.y(vt)^2) 41 fang_normal_x=wall.vertex.pos.x(vt) / mag 42 fang_normal_y=wall.vertex.pos.y(vt) / mag 44 r_vel=vector(fang_normal_x,fang_normal_y,0)*r_vel_mag 46 wall.vertex.vel(vt)=r_vel 47 endloop 50 endif 52 end where Line 23 to Line 24 defines a function to calculate the stress state of the wall; Line 26 to Line 28 checks whether the current time step is less than the global time step, if so, the servo control gain is obtained, and the current time step is updated; Line 31-Line 35 if the z-direction servo control program is enabled, the velocity in the z-direction is calculated based on the gravitational acceleration (gz), the target z-direction stress (wszz) and the current z-direction stress (tzz), and the z-direction velocity of the upper wall and the lower wall is reset; Line 37 to Line 38 if the radial servo control is enabled, the magnitude of the radial velocity is calculated based on the gravitational acceleration (gr), the target radial stress (wsrr) and the current radial stress (trr); Line 39 iterates through all the vertices of the right wall; Line 40 calculates the distance of each vertex to the origin; Line 41 to Line 42 calculates the normal component of the vertex; Line 44 calculates the radial velocity based on the normal component and r_vel_mag; Line 46 sets the velocity of the vertex.

[0037] The main code information of Line 55 to Line 77 about pre-pressing is as follows: 55 def computer_chicun 56 x_pos=wall.vertex.pos.x(vert_in_ce) 57 y_pos=wall.vertex.pos.y(vert_in_ce) 58 wlr=math.sqrt(x_pos^2+y_pos^2) 59 wlz=wall.pos.z(wp_up)-wall.pos.z(wp_down) 60 end 61 62 def computer_wallStress 63 computer_chicun 64 ding_yuanmianji = math.pi * wlr ** 2 65 wszz = (wall.force.contact.z(wp_down) - wall.force.contact.z(wp_up)) * 0.5 / ding_yuanmianji 67 ce_mianji = 2 * math.pi * wlr * wlz 68 wsrr = 0 69 loop foreach ft wall.facetlist(wp_rr) 70 ft_fangxiang = wall.facet.normal(ft) 71 loop foreach ct wall.facet.contactmap(ft) 72 force_in_facet = contact.force.global(ct) 73 wsrr += -(math.dot((force_in_facet, ft_fangxiang))) / ce_mianji 74 endloop 75 endloop 77 end Line 55 defines a function for calculating the size of the wall; Line 56 gets the x-coordinate of the specified vertex; Line 57 gets the y-coordinate of the specified vertex; Line 58 calculates the distance from the origin to the vertex, i.e. the radius of the wall; Line 59 calculates the height of the wall, i.e. the difference in z coordinates of the upper and lower walls; Line 62 defines a function for calculating the stress of the wall; Line 63 calls the function to get the size of the wall; Line 64 calculates the area of the circle where the vertex is located; Line 65 calculates the average stress in the z direction based on the difference in contact force between the upper and lower walls and the area of the circle; Line 67 calculates the lateral area of the cylinder; Line 68 initializes the radial stress to 0; Line 69 iterates through all the faces of the right wall; Line 70 gets the normal vector of the face; Line 71 iterates through all the contact points on the face; Line 72 gets the global force of the contact point; Line 73 accumulates and updates wsrr.

[0038] The main code information of Line 79 to Line 110 about pre-pressing is as follows: 79 def get_gain(fac) 80 gz=0 81 gr=0 82 zonggangZ=0 83 zonggangR=0 84 loop foreach ct wall.contactmap(wp_up) 85 zonggangZ+=contact.prop(ct,"kn") 86 endloop 88 loop foreach ct wall.contactmap(wp_down) 89 zonggangZ+=contact.prop(ct,"kn") 90 endloop 92 loop foreach ct wall.contactmap(wp_rr) 93 zonggangR+=contact.prop(ct,"kn") 94 endloop 96 gz=fac*math.pi*wlr*wlr / (zonggangZ*global.timestep) 97 gr=fac*2*math.pi*wlr*wlz / (zonggangR*global.timestep) 98 end 99 @computer_wallStress 100 @get_gain(@sevro_fac) 103 set fish callback -1.0 @servo_walls 105 history id 1 @wszz 106 history id 2 @wsrr 108 cycle 1 109 solve arat 1e-4 110 save yuya Line 79 defines a procedure called get_gain; Line 80 to Line 83 initializes the variables gz, gr, zonggangZ, zonggangR to 0; Line 84 to Line 94 iterates through all the contact points of the wall wp_up, wp_down, wp_rr, and accumulates the normal stiffness kn of each contact point to zonggangZ and zonggangR; Line 96 to Line 97 calculates the gains gz and gr in the z direction and radial direction based on the normal stiffness, wall size, time step, and parameter fac; Line 99 calls the procedure computer_wallstress to calculate the stress of the wall; Line 100 calls the get_gain procedure to obtain the gain for servo control, using the previously defined servo_fac as a parameter; Line 103 sets a callback function that calls the servo_walls procedure every 1.0 time unit of simulation, used to dynamically adjust the movement of the wall; Line 105 to Line 106 defines two history records to record the average stress wszzz and wsrr in the z direction and radial direction, respectively; Line 108 executes a simulation loop for one time step; Line 109 calls the solve command to set the relative error tolerance to 0.0001 and perform iterative solving; Line 110 saves the current simulation state as "yuya".

[0039] The main code information about the glueing process is as follows: 1 restore yuya 3 [pb_coh=10e6] 4 [ten_coh=2.7] 6 cmat default type ball-facet model linear method deformability emod10e8 kratio 1.5 9 cmat default type ball-ball model linearpbond method deformability... &10 emod 12e8 kratio 1.5 pb_deformability emod 54e8 kratio 1.5... &11 property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh]pb_fa 50 fric 0 12 cmat apply 13 clean 15 cycle 1 16 solve arat 1e-2 18 save yujiaojie 20 contact method bond gap [keli_rdmin*0.1] 22 cycle 1 23 solve arat 1e-2 25 save jiajiaojie Wherein, Line 1 restores the simulation state from a file named "yuya"; Line 3 sets the cohesion of the material (pb_coh) to 10e6; Line 4 sets the internal cohesion of the material (ten_coh) to 2.7; Line 6 sets the default contact model (cmat) to a ball-surface contact model, linear deformation model, with an elastic modulus (emod) of 10e8 and a Poisson's ratio (kratio) of 1.5; Line 9 to Line 11 sets the default contact model (cmat) to a ball-ball contact model, linear bond model, deformation model, with an elastic modulus (emod) of 12e8, a Poisson's ratio (kratio) of 1.5, a bond deformation elastic modulus of 54e8, a Poisson's ratio of 1.5, and sets the bond cohesion (pb_coh), bond tensile strength (pb_ten), bond friction angle (pb_fa), and friction coefficient (fric); Line 12 applies the previously defined contact model; Line 13 clears some temporary data or state in the simulation; Line 15 to Line 18 performs a simulation cycle of one time step, calls the solve command, sets the relative error tolerance to 0.01, performs iterative solving, and saves the current simulation state as "yujiaojie"; Line 20 sets the bond gap between contact points to 0.1 times the minimum radius of the particles; Line 22 to Line 25 again performs a simulation cycle of one time step, calls the solve command, sets the relative error tolerance to 0.01, and performs iterative solving to save the current simulation state as "jiajiaojie".

[0040] The main code information about this segment of adding confining pressure is as follows: 1 ;restore jiajiaojie 2 [tzz=-6e6] 3 [trr=-6e6] 5 cycle 1 6 solve arat 1e-2 8 sava weiya Line 1 restores the saved simulation state named "jiajiaojie"; Line 2 to Line 3 sets the stress target values in the z-direction and the radial direction to -6e6, i.e., the confining pressure of the rock mass is 6Mpa, which is an important parameter in rock mechanics research. The confining pressure parameter will be adjusted later to study the influence of different confining pressures on the cracking characteristics of rock mass under biaxial conditions; Line 5 performs a simulation cycle for one time step; Line 6 calls the solver for iterative solution, setting the relative error tolerance to 0.01; Line 8 saves the current simulation state as "weiya".

[0041] The main code information about the loading segment is as follows: 1 ;restore weiya 2 [strainRatio=10e-1] 3 [do_zservo=false] 5 ball attribute displacement multiply 0 6 wall attribute zvelocity [-wlz*strainRatio*0.5] range id 2 7 wall attribute zvelocity [wlz*strainRatio*0.5] range id 3 9 [IZ0=wlz] 10 [IR0=wlr] 11 def jiance 12 wezz=(wlz-IZ0) / IZ0 13 werr=(wlr-IR0) / IR0 14 wvol=2*werr+wezz 15 end 17 set fish callback -1.1 @jiance 20 history delete 22 history id 1 @wszz 23 history id 2 @wsrr 24 history id 3 @wezz 25 history id 4 @werr 26 history id 5 @wvol 27 call fracture.p3fis 28 @track_init 30 [stop_me=0] 31 def stop_me 32 if wezz<-20e-3 then 33 stop_me=1 34 endif 35 end 36 solve fishhalt @stop_me 37 save result Line 1 restores the saved simulation state named "weiya" before Line 1; Line 2 sets the strain ratio to 0.1; Line 3 sets a flag variable do_zservo to false, which is used to control whether to enable z-direction servo control; Line 5 multiplies the displacement attribute of all spheres by 0, i.e. resets the displacement of the spheres; Line 6 sets the velocity of the wall numbered 2, simulating the movement of the left wall; Line 7 sets the velocity of the wall numbered 3, simulating the movement of the right wall; Line 9 to Line 10 sets the variable IZ0 to the initial height of the wall wlz, and sets the variable IR0 to the initial radius of the wall wlr; Line 11 to Line 15 defines a procedure named jiance, which is used to calculate the height and radius change rate of the wall, as well as the volume change rate; Line 17 sets a callback function, which calls the jiance procedure every 1.1 time units; Line 20 deletes all historical records; Line 21 to Line 26 creates five historical records, which record the z-direction stress wszzz, the radial stress wsrr, the height change rate wezz, the radius change rate wrrr, and the volume change rate wvol, respectively; Line 27 calls the fracture.p3fis module, which is used to simulate the fracture behavior of the material; Line 28 calls the procedure named track_init, which is used to initialize the tracking variables; Line 30 initializes the variable stop_me to 0; Line 31 to Line 34 defines a procedure named stop_me, which sets stop_me to 1 if the height change rate wezz is less than -0.002; Line 36 performs the solution, i.e. stops the solution when stop_me is 1; Line 36 saves the simulation results as "result".

[0042] The beneficial effects of the above technical solution are: through the above technical solution, the multi-axial simulation compression test is completed in the PFC software, reliable data support is provided for subsequent model correction, and the accuracy of subsequent crack propagation analysis using crack distribution is further improved.

[0043] In one embodiment, when the target surrounding rock mass is subjected to indoor compression test, the method for obtaining crack distribution after compression comprises: The theoretical crack distribution range of the target surrounding rock mass after the corresponding load is applied is denoted as a priori crack area; preferably, the theoretical crack distribution range is determined through historical literature materials, and the experimental method in the historical literature materials is preferably a CT scanning test. The CT scanning test is performed on a surrounding rock mass sample (which should be consistent with the target surrounding rock mass in physical parameters) to extract a surrounding rock mass sample with the same physical and mechanical characteristics and structural characteristic parameters. The two-dimensional or three-dimensional coordinate range of the real crack after the compression test is taken as the a priori crack area; wherein the number of a priori crack areas is not less than one, that is, one or more; A corresponding binary reference mask is generated according to the a priori crack area, wherein the pixel value of the a priori crack area in the binary reference mask is a first value, and the pixel value of the non-a priori crack area is a second value; specifically including: performing grayscale and threshold segmentation processing on the image of the a priori crack area, setting the pixel value in the a priori crack area to 1 (the first value), and setting the pixel value in the non-a priori crack area to 0 (the second value); The original image of the compressed target surrounding rock mass is recognized according to the pre-trained convolutional neural network model; wherein the training method of the convolutional neural network model is a conventional training scheme, which belongs to the common knowledge of those skilled in the art, therefore, the present application will not be described in detail, in addition, the original image for recognition also needs to be preprocessed, that is, to perform Gaussian filter denoising, adaptive histogram equalization to enhance contrast, and then to perform binary segmentation through Otsu algorithm and other processing; A plurality of crack information is recognized by an edge detection algorithm, a pixel point coordinate set of the crack edge of each crack information is extracted, and a to-be-verified edge pixel set is obtained; wherein the number of to-be-verified edge pixel sets is the same as the number of crack information; Any to-be-verified edge pixel set is spatially matched with all a priori crack areas, and the number of first pixel points in the to-be-verified edge pixel set after spatial matching in each a priori crack area is counted; that is, the coordinates (x, y) of the to-be-verified edge pixel are mapped with the reference mask of the a priori crack area, and it is judged whether the pixel value corresponding to the pixel point (x, y) in the reference mask is 1, if yes, it is recorded as a matching pixel point; The ratio of each first pixel point number to the total number of pixel points in the to-be-verified edge pixel set is calculated, and is denoted as an edge pixel coincidence degree; If all edge pixel coincidence degrees of any to-be-verified edge pixel set are lower than a preset coincidence threshold, it is determined that the crack information corresponding to the to-be-verified edge pixel set is invalid crack and is removed, otherwise, it is recorded as valid crack; All valid cracks are obtained to generate the crack distribution result of this compression test; Further, the retained effective cracks can be further profile-optimized, residual noise points are removed through morphological operation (dilation, erosion), and the crack profile is simplified through polygon fitting, and a corrected crack distribution image is output for generating a crack distribution result; The technical scheme has the beneficial effects that the technical scheme can effectively solve the problem of low recognition accuracy caused by complex crack generation in compression tests, can focus on the area where cracks may theoretically be generated by introducing prior crack area constraints, can reduce the interference of non-crack area noise on the recognition result, and can reduce the probability of false recognition from the source; meanwhile, with the help of edge pixel overlap checking, the initially recognized cracks can be accurately checked one by one, the reliability of the recognition result of each crack can be quantitatively evaluated by counting the overlap ratio of the edge pixels of the crack to be checked and the prior crack area, and the false-recognized cracks can be accurately removed, so that each real crack information can be accurately recognized, and the interference of irrelevant noise or pseudo-cracks on subsequent analysis can be avoided; the subsequent comparison accuracy of the test results and the simulation results is improved, and the simulation accuracy of the simulation model for cracking simulation is improved.

[0044] In one embodiment, the surrounding rock mass cracking Fractal analysis method based on the particle flow discrete element further includes: for the first pixel points falling in the prior crack area, a weight is assigned to the first pixel points according to the Euclidean distance of the first pixel points from the center of the prior crack area, the weight is positively correlated with the distance, and a ratio calculation is performed on the number of weighted pixel points and the total number of pixel points to obtain a new edge pixel overlap degree. The technical scheme has the beneficial effects that the pixel points in the prior crack area are weighted according to the Euclidean distance from the center, which can match the distribution characteristics of real cracks, that is, the core is dense and the edge is dispersed, can reduce the misjudgment caused by edge interference, and can enhance the recognition anti-interference capability; meanwhile, the real cracks near the edge can be accurately distinguished from the pseudo-cracks, the error removal accuracy can be optimized, the space characteristics of the simulation model are more matched, and the subsequent comparison consistency is improved.

[0045] Finally, it should be pointed out that the above preferred embodiments are only used to illustrate the technical solutions of the present application and not to limit the present application, although the present application has been described in detail through the above preferred embodiments, those skilled in the art should understand that various changes can be made in form and details without departing from the scope defined by the claims of the present application.

Claims

1. A fractal analysis method for rock mass cracking based on particle flow discrete element method, characterized in that: a rock mass model is constructed according to the physical and mechanical parameters and structural characteristic parameters of the target rock mass; a virtual uniaxial compression simulation test is carried out on the rock mass model by using the particle flow discrete element software, and the simulation results are compared with the uniaxial compression test results of the target rock mass to generate an initial numerical model; a multi-axial compression test is carried out on the target rock mass, and a multi-axial simulation compression test is carried out on the initial numerical model; the test results are compared with the simulation results, and the initial numerical model is corrected according to the comparison results to establish a data model of the surrounding rock of the roadway under different stress numerical effects; actual engineering load data are obtained, and the data model of the surrounding rock of the roadway is loaded to simulate rock mass cracking, crack distribution in the simulation process is obtained, and an extension trend analysis is carried out on the crack distribution by using a fractal dimension analysis method to generate rock mass cracking analysis results. 2.The fractal analysis method for rock mass cracking based on particle flow discrete element method according to claim 1, characterized in that: a plurality of rock blocks meeting preset test conditions are obtained through field investigation as target rock masses for indoor tests; and the physical and mechanical parameters and structural characteristic parameters of the target rock masses are collected; wherein the physical and mechanical parameters include particle density, elastic modulus, Poisson's ratio, internal friction angle and cohesion; and the structural characteristic parameters include particle size distribution, porosity and fracture distribution. The generation process of the initial numerical model includes: step 1, ball constraint preparation; step 2, generation of particle material system; step 3, uniform and good contact of the model by using servo mechanism; step 4, contact coefficient assignment; step 5, control of the compression process, and analysis of the stress-strain curve; and step 6, if the stress-strain curve trend in the analysis results of step 5 is consistent with the uniaxial compression test loading trend of the target rock mass, it is determined that the current numerical simulation model is approximately matched with the pattern, and the initial numerical model is obtained. The multi-axial simulation compression test includes rigid biaxial simulation compression test, flexible biaxial simulation compression test and triaxial simulation compression test. The triaxial simulation compression test includes: step A, sample preparation: generating a uniform particle assembly with a specified porosity to simulate the sample preparation process in the actual test; step B, pre-pressing: applying an initial pressure on the sample to ensure that the internal stress of each particle is zero by using Newton's second law, so as to realize the densification of the sample, reduce the porosity of the sample, and ensure that the sample meets the initial state of the actual material; step C, add cementation: introducing cementation force between particles to simulate the bonding effect between particles, so as to improve the overall strength and stability of the sample, and ensure that the model can withstand the specified stress during the loading process; step D, add confining pressure: applying a uniform pressure around the sample to simulate the lateral pressure condition in the triaxial compression test; step E, loading: applying an axial force on the sample to simulate the compression process in the triaxial compression test; and step F, data acquisition: obtaining the stress-strain data in the compression process to obtain the simulation test results. ​ ​ ​ ​ ​ 3. The Fractal analysis method for surrounding rock mass cracking established based on particle flow discrete element according to claim 1, characterized in that, ​ ​ ​ ​ ​ ​ ​ 4. The Fractal analysis method for surrounding rock mass cracking established based on particle flow discrete element according to claim 1, characterized in that: ​ 5. The Fractal analysis method for surrounding rock mass cracking established based on particle flow discrete element according to claim 4, characterized in that, ​ ​ ​ ​ ​ ​ ​ 6. The Fractal analysis method for surrounding rock mass cracking established based on particle flow discrete element according to claim 1, characterized in that, In the indoor compression test of the target surrounding rock mass, the method for obtaining the crack distribution after compression comprises: obtaining the theoretical crack distribution range of the target surrounding rock mass after applying the corresponding load, denoted as the prior crack area, and generating the corresponding binary reference mask; wherein the number of prior crack areas is not less than one, the pixel value of the prior crack area in the binary reference mask is the first numerical value, and the pixel value of the non-prior crack area is the second numerical value; According to the pre-trained convolutional neural network model, the original image of the compressed target surrounding rock mass is identified to obtain a plurality of crack information, and the pixel point coordinate set of the crack edge of each crack information is extracted to obtain the to-be-verified edge pixel set; wherein the number of to-be-verified edge pixel sets is the same as the number of crack information. Spatially match any to-be-verified edge pixel set with all prior crack areas, and count the number of first pixel points in the to-be-verified edge pixel set after spatial matching in each prior crack area, and calculate the ratio of each first pixel point number to the total number of pixel points in the to-be-verified edge pixel set to obtain the edge pixel overlap degree. If all edge pixel overlaps of any to-be-verified edge pixel set are lower than the preset overlap threshold, it is determined that the crack information corresponding to the to-be-verified edge pixel set is invalid crack and is removed, otherwise it is denoted as valid crack. Obtain all valid cracks to generate the crack distribution result of this compression test.

7. The Fractal analysis method for surrounding rock mass cracking established based on particle flow discrete element according to claim 6, characterized in that, Further comprising: for the first pixel points falling in the prior crack area, according to the Euclidean distance between the first pixel points and the center of the prior crack area, a weight is assigned, the weight is positively correlated with the distance, and the ratio of the number of weighted pixel points to the total number of pixel points is calculated to obtain a new edge pixel overlap degree.

Images