Discrete element method-based simulation method for inherent microfractures in rock

Abstract

The invention discloses a simulation method for inherent microfractures in rock based on a discrete element method, which comprises the following steps of: constructing a discrete element numerical model sample, performing static equilibrium calculation on the sample by adopting a linear model, applying fixed confining pressure to an initial sample, and performing servo calculation; a Flat-joint contact model is added, and contact bonding mesoscopic parameters are set for static balance calculation; modifying bonding contact into non-bonding contact, adding a contact gap and performing static balance calculation; carrying out a loading test; processing a test result to obtain mechanical parameters of the sample, comparing the mechanical parameters of the numerical model test with the mechanical parameters of the rock sample obtained by the indoor test, if the two mechanical parameters are matched, outputting the result, and terminating the simulation process, otherwise, modifying the contact bonding mesoscopic parameters. And reloading the test until the mechanical parameters of the numerical test are matched with the mechanical parameters of the indoor test. According to the method, the inherent micro-fracture in the rock is simulated, and the defects of a conventional model in a discrete element method are overcome.

Description

Technical Field

[0001] The present invention relates to the field of numerical simulation methods in rock mechanics, and in particular to the simulation of the non-linear mechanical behavior of inherent microcracks in rock materials during the process of being compressed to closure, and more particularly to a simulation method based on the discrete element method for simulating the presence of inherent microcracks in rocks. Background Art

[0002] Affected by internal mineral composition, pores, microcracks and other microscopic defects, rock masses often exhibit characteristics of diversity, inhomogeneity and anisotropy. Researchers generally believe that the initial non-linear stage of the stress-strain curve observed in compression tests is closely related to the opening and closing behavior of microcracks and other defects inside the rock.

[0003] Currently, the microdefects in rock samples are usually examined through theoretical analysis and indoor or field tests. Theoretical analysis believes that compared with spherical pores, the deformation of microcracks is more sensitive to external forces, significantly promoting the sensitivity of the initial non-linear stress-strain curve.

[0004] Field investigations have shown that in the surrounding rock with stress redistribution, the modulus changes caused by joint opening and closing, and the crack interaction and failure mechanism are more complex. These studies provide an important theoretical basis and reference for a deeper understanding of the influence of microcracks on the mechanical properties of rocks. However, due to the concealment and randomness of microscopic defects, it is still difficult to directly observe the opening and closing behavior of microcracks under complex stress conditions, which limits our comprehensive understanding of their role in the macroscopic mechanical response of rocks and rock masses.

[0005] Numerical simulation provides a new way for in-depth study of rock and rock mass mechanics. Because there is no need to consider complex empirical constitutive relations, the discrete element method has been widely used in simulating and analyzing the rock failure mechanism.

[0006] Taking the particle flow software PFC as an example, this method uses explicit time-domain integral equations of motion to simulate and reproduce the processes such as crack initiation, connection and penetration during rock failure. In PFC, the parallel bond model is widely used to simulate the strength characteristics, microcrack distribution characteristics and fracture processes of intact rock masses, jointed rocks and rock joints. Previous studies have shown that in indoor tests and large-scale experiments, the parallel bond model can reproduce many characteristics of rocks, such as the elastic modulus, Poisson's ratio and uniaxial compressive strength of rocks. However, there are still some challenges, such as too low friction angle, too small ratio of uniaxial compressive strength to tensile strength and a linear strength envelope. Especially for the simulation of the microcrack closure stage of rocks, the parallel bond model has not been well solved yet.

[0007] How to solve the simulation of the discrete element method for the mechanical properties of rocks, including the problems of too low friction angle, too small ratio of uniaxial compressive strength to tensile strength, linear strength envelope, and the simulation of the nonlinear mechanical behavior caused by the gradual closure of microcracks under compression. This is of guiding value for expanding the application prospects of discrete element simulation in the field of rocks and for the stability analysis and disaster impact assessment of related rock mass engineering. Summary of the Invention

[0008] Object of the Invention: Aiming at the deficiencies of the existing discrete element simulation method in reflecting the closure process of inherent microcracks in rocks and their nonlinear mechanical behavior, the present invention proposes a simulation method for rocks with inherent microcracks based on the discrete element method. Based on the discrete element method, a numerical model of rock containing microcracks is constructed and a loading test is carried out. By setting a certain proportion of non-bonded contacts and introducing an initial contact gap in the Flat-joint contact model, the nonlinear deformation process of rocks in the initial compression stage is realistically reproduced, providing a more reliable numerical simulation method for rock mass stability analysis and rock damage mechanism research, realizing the simulation of inherent microcracks in rocks, and making up for the defects of the conventional model in the discrete element method.

[0009] Technical Solution: The simulation method for rocks with inherent microcracks based on the discrete element method of the present invention includes the following steps:

[0010] Step (1), set up rigid walls perpendicular to each other according to the size of the indoor test specimen to form a boundary constraint area (width and height are the same as the indoor specimen) matching the specimen size, and construct an initial numerical model specimen under static equilibrium. The process is as follows:

[0011] Step (1.1), generate a small number of overlapping discrete particles in this boundary constraint area. For each pair of contacting discrete particles generated, add a linear contact model, and set the normal contact stiffness k n and the tangential contact stiffness k s and the contact damping parameter to ensure the numerical stability of the contact force between particles during the iterative calculation process.

[0012] Step (1.2), perform iterative calculation and obtain the initial numerical model specimen under static equilibrium. The process is as follows:

[0013] Step (1.2.1), after fixing the position of the rigid wall, start the iterative calculation. At each iterative step, calculate the normal force F n and the tangential force F s :

[0014]

[0015] In the formula, is the normal force of the particle at time t + Δt, k nNormal contact stiffness of particle, g s is the normal overlap of particles, is the tangential force of the particle at time t, k s Tangential contact stiffness of particle, Δδ s is the incremental tangential displacement calculated for the particle.

[0016] Step (1.2.2), after calculating each iteration step, calculate the normalized unbalanced force F of the initial numerical model specimen unbal,norm , and then evaluate the current balance degree of the system:

[0017]

[0018] In the formula, F unbal,norm is the normalized unbalanced force of the calculation system, F unbal,i is the sum of all external forces and contact forces on particle i, F contact,i is all contact forces on particle i.

[0019] Step (1.2.3), when the normalized unbalanced force F of the initial numerical model specimen unbal,norm is less than the set threshold (such as 10 -5 ), that is, the system reaches static equilibrium, stop the iterative calculation and output the discrete element numerical model specimen at this time, and record the obtained specimen as the initial numerical model specimen; if not satisfied, return to step (1.2.1) to continue the iterative calculation.

[0020] Step (2), apply the target stress to the rigid wall of the obtained initial numerical model specimen, and adjust the moving speed v of the rigid wall according to the servo gain coefficient b , so that the wall stress reaches the target stress, and establish the actual in-situ stress condition;

[0021] v b = k servo (σ target - σ current ) (4)

[0022] v b is the moving speed of the rigid wall, k servo is the servo gain coefficient, σ target is the target stress on the wall, σ current is the stress on the wall.

[0023] Step (3), in the initial numerical model specimen after servo loading, modify all particle-particle linear contact models to Flat-joint contact models, set the initial mesoscopic parameters and perform static equilibrium calculations:

[0024] Step (3.1): Restart the iterative calculation. Calculate the contact forces between particles using Equations (1) and (2), and obtain the normalized unbalanced force using Equation (3). To ensure that the stress state of the specimen at this time matches the confining pressure environment set in Step (2), servo control is performed in the iterative process by combining Equation (4): If the actual boundary stress deviates from the target value, adjust the wall velocity v b , and finely adjust the boundary stress. Repeat until there is no deviation between the boundary stress and the target stress, and the system unbalanced force is lower than the threshold. Then, the specimen reaches a new static equilibrium under the Flat-joint contact model. At this time, the "specimen after servo completion" is considered to have a realistic stress field and a more reasonable mesoscopic contact relationship.

[0025] Step (4): Introduce the inherent microcracks of the rock. The process is as follows:

[0026] Step (4.1): Under the static equilibrium state, traverse the contacts between particles in the specimen and determine whether they are bonded contacts. Count the number of bonded contacts Num all . Generate a value between 0 and 1 as the non-bonded contact ratio φ in the Flat-joint contact model, and calculate and convert it to the total number of non-bonded contacts Num unbond :

[0027] Num unbond = φ.Num all (5)

[0028] In the formula: φ is the non-bonded contact ratio, Num all is the total number of bonded contacts between particles in the numerical specimen model, and Num unbond is the number of contacts set as non-bonded contacts between particles in the numerical specimen model.

[0029] Step (4.2): Select Num unbond contacts in the specimen, and change the Num unbond contacts from bonded contacts to non-bonded contacts to simulate the distribution of natural microcracks inside the rock material, and set the normal bonding strength and tangential bonding strength to zero.

[0030] Step (4.3): To further reflect the pore characteristics inside the rock, set the initial contact gap g0 in the non-bonded contacts:

[0031]

[0032] In the formula: g0 is the initial contact gap, L is the specimen length, n is the average number of particles along the specimen loading axis, γ0 is the porosity ratio, and φ is the non-bonded contact ratio.

[0033] Step (4.4), after completing the non-bonded conversion and introducing the contact gap, start the iteration again. Calculate the inter-particle contact force using Eqs. (1) and (2), and use Eq. (3) to judge whether the unbalanced force is less than the threshold value. Since introducing microcracks will affect the local force chain distribution, use Eq. (4) to adjust the boundary stress to avoid the confining pressure deviation caused by the contact change. When the numerical model is balanced again, it enters the stable state after the microcrack introduction; at this time, the specimen is in the initial mechanical environment of the real rock before loading.

[0034] Step (5), adopt a constant speed v load Control the movement of the upper and lower rigid walls to simulate the uniaxial compression test, biaxial compression and Brazilian splitting test, and monitor and calculate the axial strain, axial stress and lateral strain in real time.

[0035] Step (6), after the loading is completed, output the recorded stress-strain curve, crack distribution, and numerical test results of the failure mode, and perform data processing to extract the obtained mechanical parameters. Compare the mechanical parameters with the mechanical parameters measured in the laboratory test. If Eq. (10) is satisfied, output the results and terminate the simulation:

[0036]

[0037] where: X num is the mechanical parameter of the numerical model specimen, X exp is the mechanical parameter measured in the laboratory test, and δ is the preset error threshold (usually set according to the test accuracy and material heterogeneity, such as 5% - 10%).

[0038] Step (7), if Eq. (10) is not satisfied, there is a deviation between the two, return to Step (3) to modify the mesoscopic parameters of the Flat-joint contact model, and re-perform the loading test until the mechanical parameters obtained from the numerical test match the mechanical parameters obtained from the laboratory test.

[0039] The process of Step (2) is as follows:

[0040] Step (2.1), for the obtained initial numerical model specimen, apply the target stress σ target matching the confining pressure condition of the laboratory rock sample on the upper, lower, left and right rigid walls to simulate the magnitude of the confining pressure suffered by the laboratory rock sample in the natural state. To achieve the gradual loading of the confining pressure and stabilize it at the target stress, select the servo gain coefficient k servo , and calculate the moving speed v of the rigid wall according to the following formula b for servo control:

[0041] v b = k servo (σ target - σ current ) (4)

[0042] In the formula: v b is the current moving speed of the rigid wall, k servo is the servo gain coefficient, σ target is the target stress on the wall, σ current is the stress currently borne by the wall.

[0043] In step (2.2), when the stresses on the upper, lower, left, and right rigid walls all reach the target stress and tend to be stable, it is regarded as the completion of servo loading. The internal stress state of the initial numerical model specimen is regarded as the same as the in-situ stress conditions of the in-door specimen (or actual rock mass), which is convenient for adding the contact model and introducing microcracks in the next step.

[0044] In step (3), the mesoscopic parameters specified in the Flat-joint model are: the Young's modulus of the parallel bond the ratio of the normal stiffness to the shear stiffness of the parallel bond the particle friction coefficient μ, the tensile strength of the flat joint the cohesion of the flat joint

[0045] The process of step (5) is as follows:

[0046] In step (5.1), by recording the position changes of the upper and lower rigid walls, the axial strain ε of the numerical specimen during the loading process is obtained using the following formula y :

[0047]

[0048] In the formula, ε y is the axial strain of the numerical specimen, h0 is the initial height of the numerical specimen, and h is the height of the numerical specimen at the current moment.

[0049] In step (5.2), by recording the sum of the forces on the upper and lower rigid walls, the axial stress σ of the numerical specimen is obtained using the following formula y :

[0050]

[0051] In the formula, σ y is the axial stress of the numerical specimen, F is half of the sum of the forces on the upper and lower loading plates, and A is the area of the loading plate.

[0052] In step (5.3), by setting a measurement circle in the specimen to track the changes of particles in the lateral direction; by detecting the evolution of the diameter of the measurement circle over time, the lateral strain ε of the numerical specimen is obtained using the following formula x :

[0053]

[0054] where ε x is the lateral strain of the numerical specimen, D0 is the initial diameter of the measuring circle, and D is the diameter of the measuring circle at the current moment.

[0055] In step (6), the extracted mechanical parameters include elastic modulus, Poisson's ratio, uniaxial compressive strength, tensile strength, cohesion, internal friction angle, porosity ratio, and the ratio of the initial elastic modulus to the final elastic modulus.

[0056] In step (1), the constructed initial numerical model specimens include uniaxial compression specimens, biaxial compression specimens, and Brazilian splitting specimens.

[0057] In step (2.1), within each calculation step, the wall velocity is automatically corrected according to the difference between the actual stress and the target stress of the current wall, so that the confining pressure value within the specimen gradually approaches the set target value. When the forces on the upper, lower, left, and right rigid walls all reach and remain within the target stress range, the servo calculation is terminated.

[0058] In step (4.1), φ is the proportion of non-bonded contacts used to simulate natural microcracks in the numerical model. If φ is large, it indicates a higher density of initial microcracks in the rock; if φ is small, it indicates that the rock is more intact.

[0059] In step (4.2), by selecting the contacts at Num unbond and changing their bonded contacts to non-bonded contacts, the distribution state of natural microcracks inside the rock is reproduced macroscopically, and the zonal or weighted random method is used to realize the spatial distribution characteristics of microcracks.

[0060] In step (4.3), an initial contact gap g0 is set in the non-bonded contacts, so that these non-bonded contacts have an actual spatial gap at the initial stage of loading, thereby simulating the gradual closing process of microcracks during the compression stage.

[0061] In step (5), longitudinal loading is carried out while retaining the confining pressure environment to simulate the indoor triaxial compression test; by removing the left and right rigid walls, the indoor uniaxial compression test is simulated with the specimen being laterally free.

[0062] In step (5), during the entire process of the loading test, the axial strain, axial stress, and lateral strain are monitored and recorded in real time at each calculation time step.

[0063] In step (6), the extraction of each mechanical parameter is obtained through the following steps:

[0064] In step (6.1), for the macroscopic elastic modulus of the specimen, in the linear stage of the stress-strain curve of the uniaxial compressive test, stress values σ2 and σ1, and the corresponding strain values ε2 and ε1 are selected, and are calculated using the following formula:

[0065]

[0066] Step (6.2), for the Poisson's ratio of the specimen, during the elastic stage of the uniaxial compressive test, record the axial strain ε l and the lateral strain ε a , and calculate and obtain it using the following formula:

[0067]

[0068] Step (6.3), the compressive strength of the specimen is the peak value of the stress-strain curve of the uniaxial compression test.

[0069] Step (6.4), the tensile strength of the specimen is obtained through the ultimate load in the Brazilian splitting test.

[0070] Step (6.5), for the cohesion and internal friction angle of the specimen, based on the results of the biaxial test, draw a Mohr circle according to the failure point, and extract them using the slope and intercept of the envelope line.

[0071] Step (6.6), the void ratio of the specimen and the ratio of the initial elastic modulus to the final elastic modulus are obtained by the intersection point of the linear elastic extension line and the strain axis in the stress-strain curve of the uniaxial compression test.

[0072] Working principle: According to the dimensions of the specimens in the laboratory test, the present invention constructs a discrete element numerical model specimen, uses a linear model to perform static equilibrium calculations on the specimen to obtain the initial specimen of the numerical model specimen; applies a fixed confining pressure to the obtained initial specimen and performs servo calculations to simulate the stress state of the specimen under natural conditions; installs a Flat-joint contact model on the servoed numerical model specimen, and preliminarily sets the mesoscopic parameters of contact bonding, and performs static equilibrium calculations; modifies some bonded contacts in the specimen after the static equilibrium calculation to non-bonded contacts, and adds contact gaps to the non-bonded contacts to simulate the inherent microcracks in the rock, and performs static equilibrium calculations; conducts a loading test on the static equilibrium model specimen, and monitors and records the test results during the loading process; processes the recorded numerical test results to obtain the mechanical parameters of the numerical model specimen, and compares the mechanical parameters of the numerical model test with the mechanical parameters of the rock specimen obtained from the laboratory test to determine whether they match; if the mechanical parameters of the two match, output the numerical test results and terminate the simulation process, if they do not match, modify the mesoscopic parameters of contact bonding and re-conduct the loading test until the mechanical parameters obtained from the numerical test match the mechanical parameters obtained from the laboratory test.

[0073] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0074] (1) By introducing partial non - bonded contacts to replace full bonding, the present invention overcomes the limitations of the traditional parallel bonding model in simulating the friction angle, strength ratio, and non - linear stress - strain curve.

[0075] (2) By setting the non - bonded contact ratio and the initial contact gap, the numerical model can truly reflect the distribution of natural micro - cracks inside the rock and their closure process during the low confining pressure or initial loading stage, thus more accurately reproducing the initial non - linear segment of the rock stress - strain curve.

[0076] (3) In terms of model parameter correction, the present invention adopts a strategy of item - by - item comparison and progressive correction, and rapidly iterates under various test modes such as uniaxial, biaxial, and Brazilian splitting, improving the matching degree between the simulation results and the laboratory test results.

[0077] (4) The present invention provides an efficient and accurate numerical method for simulating the initial micro - cracks in rocks, which not only expands the application of the discrete element method in rock mechanics research and engineering, but also provides numerical support for the study of rock failure mechanisms and stability analysis. Brief Description of the Drawings

[0078] Figure 1 is a flow chart of the simulation method for the inherent micro - cracks in rocks based on the discrete element method of the present invention;

[0079] Figure 2 is a diagram of the discrete element numerical specimen established in the embodiment of the present invention;

[0080] Figure 3 is a comparison diagram of the uniaxial compression stress - strain curves of the laboratory test and the numerical test in the embodiment of the present invention. Detailed Embodiments

[0081] As Figure 1 shown, the simulation method for the inherent micro - cracks in rocks based on the discrete element method of the present invention includes the following steps:

[0082] Step (1), set up mutually perpendicular rigid walls according to the size of the laboratory test specimen to form a boundary constraint region (with the same width and height as the laboratory specimen), and construct an initial numerical model specimen under static equilibrium; as follows:

[0083] Step (1.1), generate a certain number of discrete particles that slightly overlap with each other within the boundary constraint region. For each pair of contacting discrete particles generated, add a linear contact model, and set the normal contact stiffness k n and the tangential contact stiffness k s , and at the same time define the contact damping parameter to ensure the numerical stability of the contact force between particles during the iterative calculation process.

[0084] Step (1.2), iteratively calculate and obtain the initial numerical model specimen under static equilibrium:

[0085] Step (1.2.1), after fixing the position of the rigid wall, start iterative calculation. In each iteration step, calculate the normal force F n and tangential force F s according to the particle contact:

[0086] F n (t+Δt) = k n g s (1)

[0087] F s (t+Δt) = F n (t) - k s Δδ s (2)

[0088] In the formula, F n (t+Δt) is the normal force of the particle at time t+Δt, k n is the normal contact stiffness of the particle, g s is the normal overlap of the particle, is the tangential force of the particle at time t, k s is the tangential contact stiffness of the particle, and Δδ s is the incremental tangential displacement of the particle calculation.

[0089] Step (1.2.2), after completing the calculation of each iteration step, calculate the normalized unbalanced force F unbal,norm of the numerical model specimen to evaluate the current balance degree of the system, and its calculation is obtained by the following formula:

[0090]

[0091] In the formula, F unbal,norm is the normalized unbalanced force, F unbal,i is the sum of all external forces and contact forces received by particle i, and F contact,i is all contact forces received by particle i.

[0092] Step (1.2.3), when the normalized unbalanced force F unbal,norm of the numerical model specimen is less than the set threshold (such as 10 -5 ), it is considered that the system reaches static equilibrium. If not satisfied, return to step (1.2.1) to continue iterative calculation.

[0093] Step (1.2.4), when the system reaches static equilibrium, stop iterative calculation and output the discrete element numerical model specimen at this time. The specimen obtained at this time is recorded as the initial numerical model specimen.

[0094] Among them, the initially constructed numerical model specimens include uniaxial compression specimens, biaxial compression specimens, and Brazilian splitting specimens.

[0095] Step (2), apply confining pressure to the initial numerical model and perform servo calculations to establish the actual in-situ stress conditions. The calculation process includes steps (2.1) to (2.2):

[0096] Step (2.1), for the obtained initial numerical model specimens, apply the target stress σ target that matches the confining pressure conditions of the indoor rock specimens on the upper, lower, left, and right rigid walls to simulate the magnitude of the confining pressure exerted on the indoor rock specimens in the natural state. To achieve gradual loading of the confining pressure and stabilization at the target stress, select the servo gain coefficient k servo and dynamically calculate the moving speed v of the rigid wall according to the following formula b for servo control:

[0097] v b = k servo (σ target - σ current ) (4)

[0098] In the formula: v b is the current moving speed of the rigid wall, k servo is the servo gain coefficient, σ target is the target stress exerted on the wall, and σ current is the current stress exerted on the wall.

[0099] Step (2.2), when the stresses on the upper, lower, left, and right rigid walls all reach the target stress and tend to be stable, it is regarded as the completion of servo loading. The internal stress state of the initial numerical model specimens is regarded as the same as the in-situ stress conditions of the indoor specimens (or actual rock mass), laying a foundation for the addition of the contact model and the introduction of microcracks in the next step.

[0100] Step (3), add the Flat-joint contact model and perform static equilibrium calculations again. The calculation process includes steps (3.1) to (3.2):

[0101] Step (3.1), in the numerical model specimens after completing servo loading, modify all particle-particle contacts to the Flat-joint contact model and specify their initial mesoscopic parameters.

[0102] Step (3.2), restart the iterative calculation, and also calculate the contact force between particles using Equations (1) and (2), and measure the normalized unbalanced force using Equation (3). To ensure that the stress state of the specimen at this time matches the confining pressure environment set in Step (2), servo control can be combined during the iteration process using Equation (4): If the actual boundary stress deviates from the target value, the wall velocity v b is automatically adjusted to finely tune the boundary stress. Repeat until there is no obvious deviation between the boundary stress and the target stress, and the system unbalanced force is lower than the threshold value, which indicates that the specimen reaches a new static equilibrium under the Flat-joint contact model. At this time, the specimen after "completion of servo" is considered to have a true stress field and a more reasonable mesoscopic contact relationship.

[0103] Among them, the mesoscopic parameters specified in the Flat-joint model are: Young's modulus of the parallel bond Ratio of the normal stiffness to the shear stiffness of the parallel bond Particle friction coefficient μ, tensile strength of the flat joint Cohesion of the flat joint

[0104] Step (4), introduce the inherent microcracks of the rock. The calculation process includes Steps (4.1) to (4.4):

[0105] Step (4.1), under the static equilibrium state, retrieve the contacts between particles in the specimen and determine whether they are bonded contacts, and count the number of bonded contacts Num all . Randomly generate a value between 0 and 1 as the non-bonded contact ratio φ of the Flat-joint contact model, and calculate and convert it to the total number of non-bonded contacts Num unbond :

[0106] Num unbond = φ.Num all (5)

[0107] In the formula: φ is the non-bonded contact ratio, Num all is the total number of bonded contacts between particles in the numerical specimen model, Num unbond is the number of contacts between particles in the numerical specimen model set as non-bonded contacts.

[0108] Step (4.2), select Num unbond contacts in the specimen, change them from bonded contacts to non-bonded contacts to simulate the distribution of natural microcracks inside the rock material, and reset their normal bonding strength and tangential bonding strength to zero.

[0109] Step (4.3), to further reflect the pore characteristics inside the rock, an initial contact gap g0 is set in the non-bonded contact, and its size is given by the following formula:

[0110]

[0111] where: g0 is the initial contact gap, L is the specimen length, n is the average number of particles along the loading axis of the specimen, γ0 is the void ratio, and φ is the non-bonded contact ratio.

[0112] Step (4.4), after completing the non-bonded conversion and introducing the contact gap, start the iteration again. Calculate the inter-particle contact force using Eqs. (1) and (2), and use Eq. (3) to judge whether the unbalanced force is less than the threshold. Since introducing microcracks will affect the local force chain distribution, Eq. (4) can be continuously used to make a small dynamic adjustment to the boundary stress to avoid the confining pressure deviation caused by the contact change. When the numerical model is balanced again, the stable state after introducing microcracks is completed; at this time, the specimen is close to the initial mechanical environment of the real rock before loading.

[0113] Step (5), control the movement of the upper and lower rigid walls at a constant speed v load to simulate uniaxial compression tests, biaxial compression, and Brazilian splitting tests, and monitor and calculate the axial strain, axial stress, and lateral strain in real time:

[0114] Step (5.1), by recording the position changes of the upper and lower rigid walls, the axial strain ε of the numerical specimen during loading can be obtained using the following formula y :

[0115]

[0116] where ε y is the axial strain of the numerical specimen, h0 is the initial height of the numerical specimen, and h is the height of the numerical specimen at the current moment.

[0117] Step (5.2), by recording the sum of the forces on the upper and lower rigid walls, the axial stress σ of the numerical specimen can be obtained using the following formula y :

[0118]

[0119] where σ y is the axial stress of the numerical specimen, F is half of the sum of the forces on the upper and lower loading plates, and A is the area of the loading plate.

[0120] Step (5.3), by setting a measurement circle in the specimen to track the changes of particles in the lateral direction. By detecting the evolution of the diameter of the measurement circle over time, the lateral strain ε of the numerical specimen can be obtained using the following formula x :

[0121]

[0122] In the formula, ε x is the lateral strain of the numerical specimen, D0 is the initial diameter of the measurement circle, and D is the diameter of the measurement circle at the current moment.

[0123] Among them, the longitudinal loading simulation of the triaxial compression test indoors is carried out in an environment where the confining pressure is retained. By removing the left and right rigid walls, the uniaxial compression test indoors is simulated with the specimen being laterally free. During the entire process of the loading test, the axial strain, axial stress, and lateral strain are monitored and recorded in real time at each calculation time step.

[0124] Step (6), after the loading is completed, output the numerical test results such as the recorded stress-strain curve, crack distribution, failure mode, etc.; perform data processing on it, and extract and sort out the following mechanical parameters: elastic modulus, Poisson's ratio, uniaxial compressive strength, tensile strength, cohesion, internal friction angle, porosity, and the ratio of the initial elastic modulus to the final elastic modulus. Compare the mechanical parameters obtained from the numerical simulation with the results of the indoor physical test. If the following error range is satisfied:

[0125]

[0126] then the numerical simulation is accurate, where: X num is the mechanical parameter of the specimen of the numerical model, X exp is the mechanical parameter measured in the indoor test, and δ is the preset error threshold (usually set according to the test accuracy and material heterogeneity, such as 5% - 10%).

[0127] Among them, the extraction of each mechanical parameter is obtained through the following steps:

[0128] 1) The macroscopic elastic modulus of the specimen can be calculated by selecting the stress values σ2 and σ1, and the corresponding strain values ε2 and ε1 in the linear stage of the stress-strain curve of the uniaxial compressive test, and using the following formula:

[0129]

[0130] 2) The Poisson's ratio of the specimen can be recorded for the axial strain ε l and the lateral strain ε a at the same time step (or the same stress level) in the elastic stage of the uniaxial compressive test, and calculated using the following formula:

[0131]

[0132] 3) The compressive strength of the specimen is the peak value of the stress-strain curve of the uniaxial compression test.

[0133] 4) The tensile strength of the specimen can be obtained from the ultimate load in the Brazilian splitting test.

[0134] 5) The cohesion and internal friction angle of the specimen can be extracted from the results of the triaxial test by plotting the Mohr circle based on the failure point and using the slope and intercept of the envelope line.

[0135] 6) The void ratio of the specimen and the ratio of the initial elastic modulus to the final elastic modulus can be obtained from the intersection of the linear elastic extension line and the strain axis in the stress-strain curve of the uniaxial compression test.

[0136] Step (7): When the numerical simulation results match the laboratory test results, output the numerical test results and terminate the simulation process. If there is a deviation between the two, return to step (3) to modify the mesoscopic parameters of the Flat-joint contact model and re-conduct the loading test until the mechanical parameters obtained from the numerical test match the mechanical parameters obtained from the laboratory test.

[0137] Among them, when the mechanical parameters of the laboratory test do not match the mechanical parameters of the numerical simulation, the order and basic process of modifying the mesoscopic parameters are as follows:

[0138] a. Compare the mechanical parameter E0 / E obtained from the numerical model specimen results with the E0 / E obtained from the laboratory test. If the E0 / E obtained from the numerical model specimen is larger than the E0 / E obtained from the laboratory test, decrease the bonding ratio φ at an interval of 0.02 f , otherwise increase the bonding ratio φ at an interval of 0.02 f . Substitute the newly determined bonding ratio into Equation (10) to calculate the initial contact gap g0, and conduct the uniaxial compression numerical test again with other mesoscopic parameters unchanged, so as to obtain new numerical test results and new mechanical parameter E0 / E and compare it with the E0 / E obtained from the laboratory test again. Repeat the above process until the E0 / E obtained from the numerical test is equal to the E0 / E measured in the laboratory test. At this time, the bonding ratio φ f and the initial contact gap g0 are the required values.

[0139] b. Compare the mechanical parameter elastic modulus E obtained from the numerical model specimen results with the elastic modulus E obtained from the laboratory test. If the elastic modulus E obtained from the numerical model specimen is larger than the elastic modulus E obtained from the laboratory test, decrease the mesoscopic elastic modulus otherwise increase the mesoscopic elastic modulus Conduct the uniaxial compression numerical test again with other mesoscopic parameters unchanged, so as to obtain new numerical test results and new mechanical parameter E and compare it with the E obtained from the laboratory test again. Repeat the above process until the E obtained from the numerical test is equal to the E measured in the laboratory test. At this time, the mesoscopic elastic modulus is the required value.

[0140] c. Compare the Poisson's ratio υ obtained from the numerical model specimen results with the Poisson's ratio υ obtained from the laboratory test. If the Poisson's ratio υ obtained from the numerical model specimen is greater than the Poisson's ratio υ obtained from the laboratory test, then reduce the stiffness ratio. Conversely, increase the stiffness ratio. Conduct a uniaxial compression numerical test again with other mesoscopic parameters unchanged, so as to obtain new numerical test results and new mechanical parameter υ and compare it with the υ obtained from the laboratory test again. Repeat the above process until the υ obtained from the numerical test is equal to the υ measured from the laboratory test. Then the stiffness ratio at this time is the required value.

[0141] d. Compare the ultimate tensile strength UTS obtained from the numerical model specimen results with the ultimate tensile strength UTS obtained from the laboratory test. If the ultimate tensile strength UTS obtained from the numerical model specimen is greater than the ultimate tensile strength UTS obtained from the laboratory test, then reduce the tensile strength ratio. Conversely, increase the tensile strength ratio. Conduct a Brazilian splitting numerical test again with other mesoscopic parameters unchanged, so as to obtain new numerical test results and new mechanical parameter UTS and compare it with the UTS obtained from the laboratory test again. Repeat the above process until the UTS obtained from the numerical test is equal to the UTS measured from the laboratory test. Then the tensile strength ratio at this time is the required value.

[0142] e. Compare the mechanical parameter internal friction angle obtained from the numerical model specimen results with the internal friction angle obtained from the laboratory test. If the internal friction angle obtained from the numerical model specimen is greater than the internal friction angle obtained from the laboratory test, then reduce the friction coefficient μ. Conversely, increase the friction coefficient μ. Conduct a biaxial compression numerical test again with other mesoscopic parameters unchanged, so as to obtain new numerical test results and new mechanical parameters and compare it with the obtained from the laboratory test again. Repeat the above process until the obtained from the numerical test is equal to the measured from the laboratory test. Then the internal friction angle at this time is the required value.

[0143] f. Compare the uniaxial compressive strength UCS obtained from the numerical model specimen results with the uniaxial compressive strength UCS obtained from the laboratory test. If the uniaxial compressive strength UCS obtained from the numerical model specimen is greater than the uniaxial compressive strength UCS obtained from the laboratory test, then reduce the tangential bond strength. Conversely, the tangential bond strength increases. Under the condition that other microscopic parameters remain unchanged, a uniaxial compression numerical test is carried out again to obtain new numerical test results, obtain new mechanical parameter UCS, and compare it with the UCS obtained from the laboratory test again. Repeat the above process until the UCS obtained from the numerical test is equal to the UCS measured in the laboratory test. At this time, the compressive strength UCS is the required value.

[0144] Embodiment

[0145] Taking the Baishan marble of Jinping II Hydropower Station as an example, the simulation method based on the discrete element method for the inherent microcracks in rocks of the present invention will be described in detail.

[0146] 1. Construct a discrete element numerical model specimen according to the dimensions of the laboratory test specimen (φ50mm×100mm). The numerical specimen includes a uniaxial compression specimen, a biaxial specimen, and a Brazilian splitting specimen. When constructing the numerical specimen, the minimum particle diameter is selected as 0.2mm, the maximum particle diameter is 0.332mm, and the average particle diameter is 0.26mm. Install a linear model for the particles and perform a static equilibrium calculation to obtain the initial specimen of the numerical model. The constructed numerical model specimen is as Figure 2 shown.

[0147] 2. Apply a fixed confining pressure of 1MPa to the obtained initial specimen and perform a servo calculation to simulate the 1MPa in-situ stress suffered by the specimen in the natural state.

[0148] 3. Install a contact model for the servoed numerical model specimen, and initially set the microscopic parameters of contact bonding. The set initial microscopic parameters of contact are shown in Table 1, and perform a static equilibrium calculation.

[0149] Table 1

[0150]

[0151] 4. Modify the bonded contacts in the specimen after the static equilibrium calculation to non-bonded contacts. The proportion of the non-bonded contact model is 0.4 to simulate the inherent microcracks in the rock, and calculate the initial contact gap through Equation (2). Perform a static equilibrium calculation on the numerical model.

[0152]

[0153] g0 is the initial contact gap, L is the specimen length, n is the average number of particles along the loading axis of the specimen, γ0 is the porosity ratio, and φ is the non-bonded contact ratio.

[0154] 5. Conduct a loading test on the static equilibrium model specimen, monitor and record the test results during the loading process. The monitored results include the axial strain, axial stress, and lateral strain of the specimen.

[0155] 6. Process the recorded numerical test results to obtain the mechanical parameters of the numerical model specimen. The obtained mechanical parameters include: elastic modulus E, Poisson's ratio ν, uniaxial compressive strength UCS, tensile strength BTS, cohesion c, and internal friction angle void ratio γ0, ratio of initial elastic modulus to elastic modulus E0 / E, as shown in Table 3. Compare the mechanical parameters of the numerical model test with the mechanical parameters of the rock specimen obtained from the laboratory test to determine whether they match.

[0156] Table 2

[0157]

[0158] 7. If the mechanical parameters of the numerical simulation results match those of the laboratory test results, output the numerical simulation test results and terminate the simulation process. If they do not match, return to step (3) to modify the mesoscopic parameters of contact bonding and conduct the loading test again. The order of modifying the mesoscopic parameters is as follows:

[0159] 1) Compare the mechanical parameter E0 / E obtained from the numerical model specimen results with the E0 / E obtained from the laboratory test. If the E0 / E obtained from the numerical model specimen is larger than the E0 / E obtained from the laboratory test, reduce the bonding ratio φ at an interval of 0.02 f , otherwise increase the bonding ratio φ at an interval of 0.02 f . Substitute the newly determined bonding ratio into Equation (10) to calculate the initial contact gap g0. Conduct the uniaxial compression numerical test again with other mesoscopic parameters unchanged, so as to obtain new numerical test results and new mechanical parameter E0 / E and compare it with the E0 / E obtained from the laboratory test again. Repeat the above process until the E0 / E obtained from the numerical test is equal to the E0 / E measured in the laboratory test. At this time, the bonding ratio φ f and the initial contact gap g0 are the required values.

[0160] 2) Compare the mechanical parameter elastic modulus E obtained from the numerical model specimen results with the elastic modulus E obtained from the laboratory test. If the elastic modulus E obtained from the numerical model specimen is larger than the elastic modulus E obtained from the laboratory test, reduce the mesoscopic elastic modulus , otherwise increase the mesoscopic elastic modulus Under the condition that other mesoscopic parameters remain unchanged, the uniaxial compression numerical test is carried out again to obtain new numerical test results and new mechanical parameter E, and then E is compared with the E obtained from the laboratory test again. Repeat the above process until the E obtained from the numerical test is equal to the E measured in the laboratory test. At this time, the mesoscopic elastic modulus is the required value.

[0161] 3) Compare the Poisson's ratio υ obtained from the mechanical parameters of the numerical model specimen with the Poisson's ratio υ obtained from the laboratory test. If the Poisson's ratio υ obtained from the numerical model specimen is larger than the Poisson's ratio υ obtained from the laboratory test, then reduce the stiffness ratio Conversely, increase the stiffness ratio Under the condition that other mesoscopic parameters remain unchanged, the uniaxial compression numerical test is carried out again to obtain new numerical test results and new mechanical parameter υ, and then υ is compared with the υ obtained from the laboratory test again. Repeat the above process until the υ obtained from the numerical test is equal to the υ measured in the laboratory test. At this time, the stiffness ratio is the required value.

[0162] 4) Compare the ultimate tensile strength UTS obtained from the mechanical parameters of the numerical model specimen with the ultimate tensile strength UTS obtained from the laboratory test. If the ultimate tensile strength UTS obtained from the numerical model specimen is larger than the ultimate tensile strength UTS obtained from the laboratory test, then reduce the ultimate tensile strength ratio Conversely, increase the ultimate tensile strength ratio Under the condition that other mesoscopic parameters remain unchanged, the Brazilian splitting numerical test is carried out again to obtain new numerical test results and new mechanical parameter UTS, and then UTS is compared with the UTS obtained from the laboratory test again. Repeat the above process until the UTS obtained from the numerical test is equal to the UTS measured in the laboratory test. At this time, the ultimate tensile strength ratio is the required value.

[0163] 5) Compare the mechanical parameter internal friction angle obtained from the numerical model specimen with the internal friction angle obtained from the laboratory test. If the internal friction angle obtained from the numerical model specimen is larger than the internal friction angle obtained from the laboratory test, then reduce the friction coefficient μ. Conversely, increase the friction coefficient μ. Under the condition that other mesoscopic parameters remain unchanged, the uniaxial compression numerical test is carried out again to obtain new numerical test results and new mechanical parameter and then compare it with the obtained from the laboratory test again. Repeat the above process until the obtained from the numerical test is equal to the measured in the laboratory test. At this time, the internal friction angle It is the required value.

[0164] 6) Compare the uniaxial compressive strength UCS, a mechanical parameter obtained from the numerical model specimen results, with the uniaxial compressive strength UCS obtained from the laboratory test. If the uniaxial compressive strength UCS obtained from the numerical model specimen is greater than the uniaxial compressive strength UCS obtained from the laboratory test, then reduce the tangential bond strength. Conversely, increase the tangential bond strength. Under the condition that other mesoscopic parameters remain unchanged, conduct the uniaxial compression numerical test again to obtain new numerical test results and new mechanical parameter UCS, and then compare it with the UCS obtained from the laboratory test again. Repeat the above process until the UCS obtained from the numerical test is equal to the UCS measured in the laboratory test. At this time, the uniaxial compressive strength UCS is the required value.

[0165] Until the mechanical parameters obtained from the numerical test match those obtained from the laboratory test, finally determine the mesoscopic parameters as shown in Table 3, the obtained mechanical parameters as shown in Table 2, and the uniaxial compressive curve as Figure 3 shown.

[0166] Table 3

[0167]

Claims

1. A simulation method based on the discrete element method for the presence of inherent microcracks in rocks, characterized in that: It includes the following steps: Step (1), set up rigid walls perpendicular to each other according to the specimen size to form a boundary constraint region, and construct an initial numerical model specimen under static equilibrium: Step (1.1), adding a linear contact model to the contact discrete particles generated within the boundary constraint region, and setting the normal contact stiffness k n , the tangential contact stiffness k s and the contact damping parameter; Step (1.2), perform iterative calculations to obtain the initial numerical model specimen under static equilibrium: Step (1.2.1), after fixing the position of the rigid wall, calculate the normal force F iteratively n and the tangential force F s : F n (t+Δt) = k n g s (1) F s (t+Δt) = F n (t) -k s Δδ s (2) Where F n (t+Δt) is the normal force of the particle at time t+Δt, k n is the normal contact stiffness of the particle, g s is the normal overlap of the particle, F s t is the tangential force of the particle at time t, k s the tangential contact stiffness of the particle, Δδ s is the calculated incremental tangential displacement of the particle; Step (1.2.2), calculate the normalized unbalanced force F of the numerical model specimen unbal,norm : where F unbal,norm is the normalized unbalanced force, F unbal,i is the sum of the external force and the contact force acting on particle i, and F contact,i is all the contact forces acting on particle i; Step (1.2.3), when the normalized unbalanced force F unbal,norm is less than the set threshold, static equilibrium is reached, the calculation is stopped, and the obtained specimen is recorded as the initial numerical model specimen; otherwise, return to Step (1.2.1) for iterative calculation; Step (2): Apply a target stress to the rigid wall of the obtained initial numerical model specimen, and adjust the moving speed v of the rigid wall according to the servo gain coefficient b to make the wall stress reach the target stress; v b = k servo (σ target - σ current )(4) v b is the moving speed of the rigid wall, k servo is the servo gain coefficient, σ target is the target stress on the wall, σ current is the stress on the wall; Step (3), modify the linear contact model in the initial numerical model specimen after servo loading to a Flat-joint contact model, set mesoscopic parameters and perform static equilibrium calculations: Step (3.1), perform iterative calculations. Calculate the contact force between particles through Equations (1) and (2), and obtain the normalized unbalanced force through Equation (3); during the iteration process, adjust the moving speed v of the rigid wall through Equation (4) b until the boundary stress is equal to the target stress; when the unbalanced force of the system is lower than the threshold, the specimen reaches static equilibrium; Step (4), perform non-bond conversion and introduce inherent microcracks in the rock. The process is as follows: Step (4.1), under static equilibrium, traverse the particle contacts in the specimen and count the number of bonded contacts Num all ; generate a value as the non-bonded contact ratio φ in the Flat-joint contact model and calculate the total number Num converted to non-bonded contacts unbond : Num unbond = φ.Num all (5) where: φ is the non-bonded contact ratio, Num all is the total number of bonded contacts between particles in the numerical specimen model, Num unbond is the number of contacts set as non-bonded contacts between particles in the numerical specimen model; Step (4.2), change the bonded contact to non-bonded contact at Num unbond in the specimen, and reset the normal bonding strength and tangential bonding strength to zero; Step (4.3), set the initial contact gap g0 in non-bonded contact: In the formula: g0 is the initial contact gap, L is the specimen length, n is the average number of particles along the specimen loading axis, γ0 is the void ratio, and φ is the non-bonded contact ratio; Step (4.4), start iteration, calculate the inter-particle contact force using Equation (1) and Equation (2), and use Equation (3) to judge whether the unbalanced force is less than the threshold; adjust the boundary stress to make the numerical model balanced again using Equation (4); Step (5), control the movement of the upper and lower rigid walls at a constant speed to simulate uniaxial compression tests, biaxial compression, and Brazilian splitting tests, and monitor and calculate the axial strain, axial stress, and lateral strain in real time; Step (6), after the loading is completed, output the recorded stress-strain curve, crack distribution, and failure mode numerical test results; compare the obtained mechanical parameters with the mechanical parameters measured in the laboratory test. If Equation (10) is satisfied, output the results and terminate the simulation: Where: X num is the mechanical parameter of the numerical model specimen, X exp is the mechanical parameter measured by the laboratory test, and δ is the preset error threshold; Step (7), if Equation (10) is not satisfied, return to Step (3) to modify the mesoscopic parameters of the Flat-joint contact model and perform the loading test.

2. The simulation method based on the discrete element method for the simulation of inherent microcracks in rocks according to claim 1, characterized in that: The process of Step (2) is as follows: Step (2.1): Apply the target stress σ matching the confining pressure of the indoor rock sample to the rigid wall of the obtained initial numerical model specimen target ; Select the servo gain coefficient k servo , and calculate the moving speed v of the rigid wall b : v b = k servo (σ target - σ current )(4) Where: v b is the moving speed of the rigid wall, k servo is the servo gain coefficient, σ target is the target stress on the wall, σ current is the current stress on the wall; Step (2.2), when the stress on the rigid wall reaches the target stress, the servo loading is completed; the internal stress of the initial numerical model specimen is the same as the in-situ stress on the laboratory specimen.

3. The simulation method based on the discrete element method for the simulation of inherent microcracks in rock according to claim 1, wherein: In step (3), the mesoscopic parameters include the Young's modulus of the parallel bonds The ratio of the normal stiffness to the shear stiffness of the parallel bond The particle friction coefficient μ, the tensile strength of the smooth joint and the cohesion of the smooth joint 4. The simulation method based on the discrete element method for the simulation of inherent microcracks in rocks according to claim 1, characterized in that: In Step (4.1), under static equilibrium, generate a value between 0 and 1 as the non-bonded contact ratio φ in the Flat-joint contact model.

5. The simulation method based on the discrete element method for the presence of inherent microcracks in rocks according to claim 1, characterized in that: The process of Step (5) is as follows: Step (5.1), by recording the position changes of the upper and lower rigid walls, the axial strain ε of the numerical specimen during loading is obtained y : where ε y is the axial strain of the numerical specimen, h0 is the initial height of the numerical specimen, and h is the height of the numerical specimen at the current moment; Step (5.2), by recording the sum of the forces on the upper and lower rigid walls, the axial stress σ of the numerical specimen is obtained y : where, σ y is the axial stress of the numerical specimen, F is half of the sum of the forces applied to the upper and lower loading plates, and A is the area of the loading plate; Step (5.3), the change of particles in the lateral direction is traced by setting a measurement circle in the specimen; by detecting the evolution of the measurement circle diameter over time, the lateral strain ε of the numerical specimen is obtained. x : where ε x is the lateral strain of the numerical specimen, D0 is the initial diameter of the measuring circle, and D is the diameter of the measuring circle at the current moment.

6. The simulation method based on the discrete element method for the simulation of inherent microcracks in rocks according to claim 1, characterized in that: In Step (6), the extracted mechanical parameters include elastic modulus, Poisson's ratio, uniaxial compressive strength, tensile strength, cohesion, internal friction angle, void ratio, and the ratio of the initial elastic modulus to the final elastic modulus.

7. The simulation method based on the discrete element method for the simulation of inherent microcracks in rocks according to claim 1, characterized in that: In step (7), the process of returning to step (3) to modify the mesoscopic parameters of the Flat-joint contact model is as follows: when the E0 / E of the numerical model specimen is greater than that of the laboratory test, the bonding ratio φ is decreased at intervals; otherwise, φ is increased at intervals. f Substitute the newly determined bonding ratio into Equation (10) to calculate the initial contact gap g0. Conduct a uniaxial compression numerical test until the E0 / E of the numerical test is equal to that of the laboratory test, and obtain the modified bonding ratio φ f and the initial contact gap g0. f ​ 8. The simulation method based on the discrete element method for the presence of inherent micro-cracks in rocks according to claim 1, characterized in that: In step (7), the process of returning to step (3) to modify the mesoscopic parameters of the Flat-joint contact model is as follows: If the elastic modulus E obtained from the numerical model specimen is greater than the elastic modulus E of the laboratory test, then decrease the mesoscopic elastic modulus Conversely, increase it Conduct the uniaxial compression numerical test again until the E of the numerical model specimen is equal to the E of the laboratory test, and obtain the modified mesoscopic elastic modulus 9. The simulation method based on the discrete element method for the simulation of inherent microcracks in rocks according to claim 1, characterized in that: In step (7), the process of returning to step (3) to modify the mesoscopic parameters of the Flat-joint contact model is as follows: If the Poisson's ratio υ of the numerical model specimen is greater than the Poisson's ratio υ of the laboratory test, then reduce the stiffness ratio Conversely, increase it Conduct a uniaxial compression numerical test again until the υ of the numerical model specimen is equal to the υ measured in the laboratory test, and obtain the modified stiffness ratio 10. The simulation method based on the discrete element method for the presence of inherent microcracks in rocks according to claim 1, characterized in that: In step (7), the process of returning to step (3) to modify the mesoscopic parameters of the Flat-joint contact model is as follows: If the ultimate tensile strength (UTS) of the numerical model specimen is greater than the UTS of the laboratory test, then reduce the tensile strength ratio. Conversely, increase it. Conduct the Brazilian splitting numerical test again until the UTS of the numerical model specimen is equal to the UTS measured in the laboratory test, and obtain the modified tensile strength ratio.

Images