Tunnel blasting shock disturbance load precision simulation method and device

Abstract

The application discloses a tunnel blasting impact disturbance load precise simulation method and device, and relates to the technical field of blasting load simulation. The method comprises the following steps: obtaining rock mass physical data and joint fissure geometric characteristics of a tunnel region to be simulated, establishing a three-dimensional finite element model, and performing local grid encryption on a blast hole influence region; applying blast hole wall pressure time history data, simulating equivalent plastic strain and plastic strain rate of a grid, and defining an elastoplastic constitutive model of damage evolution based on the same; generating a blasting response data set by using a smoothed particle hydrodynamics method, and applying the same to the model grid; dynamically updating rock mass permeability and solving a gas seepage control equation, realizing coupling of blasting impact and gas driving, and completing load simulation. The method can more accurately simulate a gas seepage process, and then realize coupling of gas driving and blasting impact, effectively improves the scientificity and accuracy of tunnel blasting impact disturbance load simulation, and effectively improves the safety and economy of a tunnel project.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of blasting load simulation, in particular to a tunnel blasting impact disturbance load precise simulation method and device. BACKGROUND

[0002] Tunnel blasting technology has important application value in engineering construction, especially in underground engineering and mining, blasting operation is a common operation method. However, in the process of traditional blasting design and implementation, there are often some technical problems. For example, how to effectively predict the influence of blasting operation on the surrounding rock mass, how to reasonably evaluate and control the surface deformation and damage degree caused by blasting, how to optimize the blasting parameters and design to achieve the best economic benefit, etc. The existence of these problems often leads to uncertainty in blasting operation, which affects the safety and economy of the project.

[0003] In the field of tunnel blasting simulation, the existing technology mainly relies on theoretical models and empirical formulas. Although these methods can provide preliminary assessment of blasting effects to some extent, they often cannot truly reflect the blasting response under complex rock mass conditions. Especially in the presence of irregular structures such as joints and fractures, traditional numerical simulation methods cannot fully consider the influence of these factors on the blasting process and effect. At the same time, the influence of the gas pressure generated by blasting on the tunnel is ignored, resulting in a large error in the simulation effect. Therefore, there is a lack of a method that can accurately simulate the tunnel blasting impact disturbance load, which limits the analysis and decision-making of engineers on the blasting process.

[0004] In addition, the traditional finite element model usually has difficulty in realizing grid densification in local areas during blasting simulation, resulting in insufficient accuracy of stress distribution and dynamic response in key areas. At the same time, existing damage evolution models are mostly based on simple elastic or linear plastic models, lacking in-depth analysis of the evolution of material microstructure under complex stress states, such as the action of blasting gas, and cannot effectively capture the true damage evolution characteristics of materials during blasting. Therefore, an innovative technical solution that comprehensively considers the structural characteristics of rock mass and blasting parameters and can accurately simulate the coupling effect of impact load and gas seepage is urgently needed to improve the scientificity and rationality of tunnel blasting design.

[0005] In the prior art, a blasting disturbance impact load precision simulation device is disclosed in CN116384052A, which comprises a blasting data acquisition module, a blasting monitoring module, a main control module, a blasting impact calculation module, a blasting simulation module, a vibration prediction module, an evaluation module and a display module. The vibration prediction module uses Matlab program numerical simulation, and the evaluation module introduces the charge distance into the evaluation of the influence of blasting vibration on buildings, so as to calculate whether the peak vibration acceleration of building particles meets the blasting vibration safety control standard, and to determine the influence degree of buildings around the target blasting area in different charge distance ranges by using statistics. However, this scheme cannot directly reflect the damage of the blasting area, cannot reflect the accurate blasting simulation in the case of irregular structures such as joints and cracks, and does not consider the influence of the gas pressure generated by blasting on the tunnel, so that the accuracy and effectiveness of the blasting impact disturbance load simulation are reduced.

[0006] The above information disclosed in the background section is only intended to enhance the understanding of the background of the present disclosure, and therefore it can include information that does not constitute the prior art known to those of ordinary skill in the art. SUMMARY

[0007] The present application aims to provide a tunnel blasting impact disturbance load precision simulation method and device to solve the problems in the background art.

[0008] To achieve the above-mentioned purpose, the present application provides the following technical solutions:

[0009] A tunnel blasting impact disturbance load precision simulation method, comprising the following specific steps:

[0010] Obtain the tunnel rock mass physical data and joint fracture geometric characteristics of the tunnel area to be simulated, establish a three-dimensional finite element model of the tunnel area to be simulated based on the tunnel rock mass physical data, and implement local grid encryption in the blast hole influence area based on the joint fracture geometric characteristics;

[0011] Apply the blast hole wall pressure time history data to the established three-dimensional finite element model, simulate the equivalent plastic strain and plastic strain rate time evolution data of each grid, and define an elastoplastic constitutive model of damage evolution based on the equivalent plastic strain and plastic strain rate time evolution data;

[0012] Simulate the damage of the grid rock mass based on the elastoplastic constitutive model, characterize the blasting impact by combining the smoothed particle hydrodynamics method, generate a blasting response data set, and apply the data in the blasting response data set to the corresponding grid of the three-dimensional finite element model;

[0013] Based on the damage variable evolution process output by the elastic-plastic constitutive model, the rock mass permeability is dynamically updated and the gas seepage control equation is solved, and the obtained gas pressure is applied as an additional load on the surface of the damage element to realize the coupling effect of blasting impact and gas driving, and complete the simulation of the blasting impact disturbance load of the tunnel.

[0014] Further, the tunnel rock mass physical data specifically includes the density, elastic modulus and compressive strength of the rock mass, and the joint fracture geometric features include the length, width, depth and distribution position of the fractures on the rock mass.

[0015] Based on the tunnel rock mass physical data, a three-dimensional finite element model of the tunnel region to be simulated is established, wherein a three-dimensional modeling software is used to construct a geometric model of the tunnel to determine the shape, size and depth of the tunnel, the tunnel rock mass physical data is embedded into the geometric model, a three-dimensional finite element model is generated by the modeling software, the entire model is meshed, the grid geometry is selected, and according to the joint fracture geometric feature data, specifically the fracture length, width, depth and distribution position, local grid densification is implemented in the blast hole influence area.

[0016] Further, the logic for implementing local grid densification in the blast hole influence area based on the joint fracture geometric features is as follows: the blast hole influence radius is determined, and the blast hole influence area is determined based on the radius, wherein the blast hole influence area is a circular area divided by taking the blast hole center as the center and the blast hole influence radius as the radius.

[0017] The grid in the blast hole influence area is densified, specifically by reducing the area size of the grid in the blast hole influence area, wherein the formula for densifying the grid is as follows:

[0018]

[0019] wherein, is the area size of the grid in the blast hole influence area, is the initial grid area size, , and are the average length, average width and average depth of the fractures in the blast hole influence area, is the fracture volume reference value, is the number of fractures in the blast hole influence area, is the fracture number reference value, and are the weight coefficients of the fracture volume and the number of fractures, respectively, wherein and and are both greater than 0.

[0020] Further, an elastic-plastic constitutive model for defining the damage evolution is defined, wherein the elastic-plastic constitutive model is specifically based on the formula:

[0021]

[0022] wherein, D a (t) is a damage variable of the a-th grid at the time t, D a (0) is a damage state of the a-th grid before the blasting, e p a (t) is an equivalent plastic strain of the a-th grid at the time t, p a (t) is a plastic strain rate of the a-th grid at the time t, and D a (t) is a damage variable of the a-th grid at the time t, The compressive strength of the tunnel rock mass and the fracture characteristics of the tunnel rock mass are characterized, and the specific formula is as follows:

[0023]

[0024] wherein, D a (0) is a damage state of the a-th grid before the blasting, D (0) is the overall compressive strength of the tunnel rock mass material before the blasting, D f is a tunnel fracture influence coefficient, N a is the total number of fractures in the a-th grid;

[0025] The established three-dimensional finite element model is subjected to the blast hole wall pressure time history data, and the specific formula for calculating the blast hole wall pressure time history is as follows:

[0026]

[0027] wherein, P (t) is the blast hole wall pressure at the time t, which is the pressure received by the blast hole wall at the time t due to the blasting, P p is a blasting peak pressure, D is a decay coefficient, T is a blasting pulse width, wherein the blasting peak pressure The specific formula for calculation is as follows:

[0028]

[0029] wherein, D is a density of the explosive, V is a detonation velocity of the explosive;

[0030] The blasting pulse width The formula used for the calculation is:

[0031]

[0032] In the formula, is the charge diameter.

[0033] Further, the smooth particle hydrodynamics method is used to simulate the blasting shock to generate a blasting response dataset including the stress tensor and velocity vector of the rock element in each grid in the finite element model, wherein the formula used for the calculation of the velocity vector of the rock element in the grid is:

[0034]

[0035] In the formula is the velocity vector of the cth rock element in the grid at time t, is the mass of the bth rock element in the neighborhood of the cth rock element, is the density of the cth rock element, is the density of the bth rock element in the neighborhood of the cth rock element, is the artificial viscosity term of the cth rock element and the bth rock element in the neighborhood thereof, is the position gradient of the cth rock element, is the kernel function, and are the Cauchy stresses borne by the cth rock element and the bth rock element in the neighborhood thereof at time t, wherein b is the index of the rock element in the neighborhood of the cth rock element, wherein is the total number of rock elements in the neighborhood of the cth rock element, and c is the index of the rock element in the grid, and the neighborhood of the cth rock element is the set of rock elements adjacent to the cth rock element in the same grid;

[0036] wherein the Cauchy stress borne by the cth rock element at time t is The damage variable is specifically calculated, and the formula used for the calculation is:

[0037]

[0038] In the formula, is the damage variable of the cth rock element at time t, is the elastic stiffness tensor of the cth rock element, is the elastic strain tensor of the cth rock element, which is specifically characterized by finite element analysis, wherein the vector sum of the Cauchy stresses borne by all rock elements in the grid is taken as the stress tensor of the corresponding grid.

[0039] Further, the formula for dynamically updating the rock mass permeability is as follows:

[0040]

[0041] wherein, is the rock mass permeability of the a-th grid at time t, is the initial rock mass permeability of the a-th grid, is the damage-permeability correlation index;

[0042] The formula for solving the gas seepage control equation is as follows:

[0043]

[0044] wherein, is the gas seepage velocity vector of the a-th grid at time t, is the dynamic viscosity of the gas, is the gas pressure suffered by the a-th grid at time t, is the gas pressure gradient of the a-th grid, is the porosity of the a-th grid, is the gas density, is the molar mass of the gas, is the universal gas constant, is the gas temperature, is the divergence operator;

[0045] By solving the gas seepage control equation, the gas pressure suffered by the a-th grid at time t is obtained, and the obtained gas pressure is applied as an additional load on the surface of the damage unit, thereby realizing the coupling effect of the blasting impact and the gas driving, wherein the damage unit refers to the surface grid unit of the tunnel region to be simulated.

[0046] The application further provides a tunnel blasting impact disturbance load precise simulation device, which is used for executing the tunnel blasting impact disturbance load precise simulation method.

[0047] The tunnel analysis model establishing module is used for acquiring the tunnel rock mass physical data and the joint fracture geometric characteristics of the tunnel region to be simulated, establishing a three-dimensional finite element model of the tunnel region to be simulated based on the tunnel rock mass physical data, and implementing local grid encryption based on the joint fracture geometric characteristics in the blast hole influence region.

[0048] The blasting response data analysis module is used for applying the blast hole wall pressure time history data to the established three-dimensional finite element model, simulating the equivalent plastic strain and plastic strain rate time evolution data of each grid, and defining the elastoplastic constitutive model of damage evolution based on the equivalent plastic strain and plastic strain rate time evolution data;

[0049] Based on the elastoplastic constitutive model, the damage of the grid rock mass is simulated, and the blasting impact is characterized by combining the smoothed particle hydrodynamics method, a blasting response data set is generated, and the data in the blasting response data set is applied to the corresponding grid of the three-dimensional finite element model;

[0050] The blasting gas additional load simulation module is used for dynamically updating the rock mass permeability based on the damage variable evolution process output by the elastoplastic constitutive model and solving the gas seepage control equation, and the obtained gas pressure is applied as an additional load to the surface of the damage unit, so as to realize the coupling effect of blasting impact and gas driving, and complete the simulation of the tunnel blasting impact disturbance load.

[0051] Compared with the prior art, the blasting method has the beneficial effects that:

[0052] By acquiring the rock mass structure characteristic data of the region to be simulated, a three-dimensional finite element model is constructed, the effective description of the complex rock mass condition is realized, the simulation accuracy of the blasting response is improved by combining the elastoplastic constitutive model and the smoothed particle hydrodynamics method, and the dynamic behavior in the blasting process is understood. The implementation of the method can solve the limitations existing in the traditional simulation to a great extent, and provide more reliable blasting design basis.

[0053] Meanwhile, based on the elastoplastic constitutive model of damage evolution, the rock mass permeability is dynamically updated, the gas seepage process can be more accurately simulated, and the coupling effect of gas driving and blasting impact is realized, thereby effectively improving the safety and economy of the tunnel engineering. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 It is a whole method flowchart of the application;

[0055] Figure 2 It is a damage variable change point line graph before and after blasting;

[0056] Figure 3 It is a fitting curve of rock mass compressive strength and total number of fractures;

[0057] Figure 4 It is a fitting curve of rock mass compressive strength and damage state before blasting;

[0058] Figure 5 It is a fitting curve of total number of fractures and damage state before blasting;

[0059] Figure 6 The figure is a schematic diagram of the device structure of the present application. DETAILED DESCRIPTION

[0060] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with specific examples.

[0061] It should be noted that, unless otherwise defined, the technical terms or scientific terms used in the present application should be understood as the general meaning understood by those skilled in the art to which the present application belongs. The terms "first", "second" and similar terms used in the present application do not represent any order, quantity or importance, but are only used to distinguish different components. The terms "include" or "contain" and similar terms mean that the elements or objects before the terms cover the elements or objects listed after the terms and their equivalents, without excluding other elements or objects. The terms "connect" or "connected" and similar terms are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. The terms "up", "down", "left", "right" and the like are only used to represent relative positional relationships, which may change accordingly when the absolute position of the described object changes.

[0062] Embodiment:

[0063] Please refer to Figures 1-5 The present application provides a technical solution:

[0064] A precise simulation method for tunnel blasting impact disturbance load, the specific steps comprising:

[0065] Step 1: Obtain the tunnel rock mass physical data and joint fracture geometric characteristics of the tunnel area to be simulated, establish a three-dimensional finite element model of the tunnel area to be simulated based on the tunnel rock mass physical data, and implement local grid encryption based on the joint fracture geometric characteristics in the blast hole influence area.

[0066] The tunnel rock mass physical data specifically includes the density, elastic modulus and compressive strength of the rock mass, and the joint fracture geometric characteristics include the length, width, depth and distribution position of the fractures on the rock mass.

[0067] The process of obtaining the tunnel rock mass physical data and joint fracture geometric characteristics of the tunnel area to be simulated generally includes the following steps: field investigation: first, conduct detailed field investigation, obtain the geological background and structural characteristics of the region through geological mapping, geological profile analysis and other means.

[0068] Drilling and sampling: collect core samples at different depths in the tunnel area through drilling technology. The commonly used drilling methods include core drilling or rotary drilling.

[0069] Physical property testing:

[0070] Density: The density of the rock can be measured by water immersion or gas displacement methods. Elastic Modulus: The elastic modulus is determined using static and dynamic testing methods, such as uniaxial compression testing and ultrasonic velocity testing. Compressive Strength: The compressive strength of the rock is determined through uniaxial compression experiments.

[0071] The specific ways to obtain the geometric characteristics of joint fissures are: Fissure length, width, and depth: The geometric characteristics of the fissures can be obtained through visual observation and measurement of the core, or using high-resolution three-dimensional laser scanning technology, CT scanning, and other means. Distribution location: The spatial distribution of the fissures can be obtained through geological profile maps, measurement records, and three-dimensional modeling.

[0072] The measured physical data and geometric characteristics are organized to establish a database. Statistical analysis is performed on the rock physical properties and joint fissure characteristics to identify their distribution patterns for subsequent model establishment.

[0073] A three-dimensional finite element model of the tunnel area to be simulated is established based on the tunnel rock physical data. The three-dimensional modeling software is used to construct the geometric model of the tunnel, determine the shape, size, and depth of the tunnel, embed the tunnel rock physical data into the geometric model, generate a three-dimensional finite element model through the modeling software, divide the entire model into grids, select the grid geometry, and according to the geometric characteristics of the joint fissures, specifically the length, width, depth, and distribution location of the fissures, implement local grid densification in the blast hole influence area.

[0074] Commonly used three-dimensional modeling software includes ABAQUS, ANSYS, SolidWorks, COMSOL Multiphysics, etc. The three-dimensional modeling software is used to construct the geometric model of the tunnel, determine whether the tunnel is circular, rectangular, or other shapes, and input the diameter, width, height, and depth of the tunnel according to design requirements. The geometric modeling function of the modeling software is used to draw the shape of the tunnel. The sketching tool can be used, then stretched or rotated to generate a three-dimensional shape. In the modeling software, create or import a material library. Set the material properties for the rock of the tunnel, input the density, elastic modulus, and compressive strength of the rock mass, etc. According to the different regions of the rock mass, embed the corresponding physical data into the geometric model. Ensure that each material region reflects the actual physical properties.

[0075] According to the collected data, create the geometric characteristics of the fissures.

[0076] The length, width, and depth of the fissures can be represented by drawing line segments, rectangles, or using other tools. Combine the fissure geometric model with the rock mass model to ensure that the distribution location of the fissures matches the actual situation. According to the complexity of the model and the calculation requirements, select appropriate grid types such as tetrahedron, hexahedron, or mixed grid.

[0077] The entire model is preliminarily meshed using the meshing tool provided by the software.

[0078] The mesh density is ensured to be moderate so as to retain sufficient calculation accuracy in the required area. Local mesh encryption is implemented in the blast hole influence area. By adjusting the meshing parameters, the mesh in these areas is made finer, thereby improving the calculation accuracy.

[0079] The logic for implementing local mesh encryption in the blast hole influence area based on the geometric characteristics of the joint fissure is as follows: the blast hole influence radius is determined, and the blast hole influence area is determined based on the radius, wherein the blast hole influence area is a circular area divided by taking the center of the blast hole as the center and the blast hole influence radius as the radius.

[0080] The mesh in the blast hole influence area is encrypted, specifically by reducing the area size of the mesh in the blast hole influence area, wherein the formula for encrypting the mesh is as follows:

[0081]

[0082] In the formula, A is the area size of the mesh in the blast hole influence area, A0 is the initial mesh area size, L is the average length of the fissure in the blast hole influence area, W is the average width of the fissure in the blast hole influence area, D is the average depth of the fissure in the blast hole influence area, V is the fissure volume reference value, N is the number of fissures in the blast hole influence area, N0 is the fissure number reference value, and a and b are the weight coefficients of the fissure volume and the fissure number, respectively. A is the area size of the mesh in the blast hole influence area, A0 is the initial mesh area size, , L is the average length of the fissure in the blast hole influence area, W is the average width of the fissure in the blast hole influence area, D is the average depth of the fissure in the blast hole influence area, V is the fissure volume reference value, N is the number of fissures in the blast hole influence area, N0 is the fissure number reference value, a is the weight coefficient of the fissure volume, and b is the weight coefficient of the fissure number. Both a and b are greater than 0. It should be noted that the area size of the mesh in the blast hole influence area A is used to represent the set mesh area in the blast hole influence area, and different mesh areas represent the complexity of the stress analysis of the area.

[0083] The smaller the value, the more accurate the analysis of the area. In the formula, A is the area size of the mesh in the blast hole influence area, A0 is the initial mesh area size,

[0084] L is the average length of the fissure in the blast hole influence area, W is the average width of the fissure in the blast hole influence area, D is the average depth of the fissure in the blast hole influence area, respectively, are the average length, average width and average depth of the fractures in the region, which directly reflect the geometric characteristics of the fractures, affecting the strength and stability of the rock mass. The existence and characteristics of the fractures determine the failure mode and mechanical behavior of the rock mass, and the larger and more fractures in the region represent more complex stress distribution, thus requiring finer meshes to capture these complexities.

[0085] If the total volume of the fractures is greater than the reference value, the mesh area will be reduced , making the mesh finer, and vice versa. Through this mechanism, mesh refinement corresponds to the increase in fracture volume, which can better reflect the influence of fractures on tunnel stability.

[0086] The mesh density is dynamically adjusted by the change in the number of fractures . If the number of fractures in a certain region exceeds the reference value, it will further affect the mesh area, indicating that the potential for damage in that region has increased, and finer meshes are needed to improve the accuracy of the simulation.

[0087] The volume of the fracture directly reflects the width and depth of the fracture in the rock mass, and has a significant impact on the overall strength, stiffness and deformation characteristics of the rock. Larger fracture volume may lead to a decrease in the compressive strength of the rock mass, making it more prone to failure when subjected to load. When the fracture volume increases, the strength and stability of the rock mass decrease significantly. Therefore, when meshing, the change in fracture volume should be given priority to in order to set finer meshes in that area to capture more complex mechanical behavior. In the field of tunnels and underground engineering, safety is the primary consideration. Prioritizing the influence of fracture volume on meshing can better reflect potential risk areas in the model, thereby improving the accuracy of the analysis. If the weight of the number of fractures is too high, it may lead to excessively fine meshing, increasing the calculation time and complexity without corresponding improvement in accuracy, so the weight of the number of fractures and are both greater than 0.

[0088] wherein the fracture volume reference value and the number of fractures reference value are both set by expert experience.

[0089] Step 2: Apply the blast hole wall pressure time history data to the established three-dimensional finite element model to simulate the equivalent plastic strain and plastic strain rate time evolution data of each mesh, and define the elastoplastic constitutive model of damage evolution based on the equivalent plastic strain and plastic strain rate time evolution data.

[0090] Define the elastoplastic constitutive model of damage evolution, wherein the formula of the elastoplastic constitutive model is:

[0091] ​

[0092] wherein, D a(t) is the damage variable of the a-th grid at time t, D a(0) is the damage state of the a-th grid before the blast, e a(t) is the equivalent plastic strain of the a-th grid at time t, r a(t) is the plastic strain rate of the a-th grid at time t, and D a(t) = D a(0) + (1 - D a(0)) * exp(-a * t) where a is the index of the grid in the three-dimensional finite element model, where t is the time variable during the blast simulation test, specifically the evolution of the blast load, wherein the duration of the blast simulation test is generally designed to be between 10 milliseconds and 100 milliseconds, and the starting time of the blast is taken as the starting point of the blast simulation test.

[0093] wherein D a(t) is the damage variable of the a-th grid at time t, This value reflects the degree of damage of the grid after experiencing the explosion load, and is generally in the range of 0 to 1, 0 indicating no damage and 1 indicating complete damage;

[0094] D a(0) is the damage state of the a-th grid before the blast, Indicates the initial condition, which can generally be regarded as the state of the model before being subjected to load;

[0095] e a(t) is the equivalent plastic strain of the a-th grid at time t, is a quantity that describes the plastic deformation of the material under loading, reflecting the deformation capacity of the material and the initial indication of damage;

[0096] r a(t) is the plastic strain rate of the a-th grid at time t, The plastic strain rate refers to the rate of plastic deformation of the material per unit time, reflecting the response of the material under dynamic loading conditions.

[0097] e a(t) is the equivalent plastic strain of the a-th grid at time t, is the main deformation indicator of the material after experiencing external load. An increase in this variable generally means that the internal molecular structure of the material has been damaged more, so it is selected to be included in the formula, reflecting the direct impact of the deformation state on damage. Plastic strain rate provides information about the response of the material under transient loading. High strain rates often lead to faster damage evolution, so the influence is measured in the form of , allowing the model to reflect the accelerated process of material damage under high strain rate conditions.

[0098] using the exponential function form and is to capture the nonlinear characteristics of material damage. The evolution of damage is often nonlinear, especially under high strain conditions, where the material may exhibit a sudden response. This form can more accurately simulate the damage distribution of the material under different loading conditions.

[0099] By adjusting the values of and , the sensitivity of the model to different plastic strains and plastic strain rates can be flexibly controlled. This enables the model to be optimized according to different material characteristics and actual engineering requirements, improving its adaptability. By consulting relevant literature and research in the field, recommended parameter values for material characteristics and damage models can be found. Many studies provide ranges of sensitivity coefficients for different materials, such as concrete, metals, and rocks, under specific conditions, generally between 0.01 and 0.3, generally between 0.01 and 0.5.

[0100] where the equivalent plastic strain and plastic strain rate are obtained from the established three-dimensional finite element model through finite element simulation. The specific method is as follows: after the solution is completed, the stress and strain distribution data of each grid element are extracted. The equivalent plastic strain can usually be obtained directly through the post-processing function provided by the software, or it can be calculated according to the stress state. The equivalent plastic strain is usually obtained through the definition of plastic strain and the evolution model of stress;

[0101] The plastic strain rate is calculated through post-processing data, which is the result of the time derivative. Time series analysis is required to obtain the plastic strain rate at different time points.

[0102] where is characterized by the compressive strength of the tunnel rock mass and the fracture characteristics of the tunnel rock mass. The specific formula is:

[0103]

[0104] where, is the compressive strength of the a-th grid before blasting, is the overall compressive strength of the tunnel rock mass material before blasting, is the tunnel fracture influence coefficient, is the total number of fractures in the a-th grid.

[0105] It should be noted that by calculating , the difference between the compressive strength of the a-th grid and the overall compressive strength can be obtained. This difference can reflect the local strength characteristics of the grid, which in turn affects its damage state. If the local compressive strength is significantly lower than the overall strength, it indicates that this area may have a higher risk of damage.

[0106] Using the form of relative difference , the damage state of different grids can be standardized and made comparable. This setting allows us to more easily assess the damage level between different grids without being affected by absolute values.

[0107] The presence of fissures usually significantly reduces the strength of rock mass, especially under high stress conditions. By introducing and the tunnel fissure influence coefficient , the formula can dynamically reflect the influence of fissure quantity on rock mass strength. In particular, the use of a logarithmic function can better capture the nonlinear impact of fissure quantity on damage: as the number of fissures increases, the improvement of damage state is gradually accelerated, where is set by expert experience.

[0108] Step 3: Based on the elastoplastic constitutive model, simulate the damage of the grid rock mass, and combine the smoothed particle hydrodynamics method to characterize the blasting impact, generate the blasting response dataset, and apply the data in the blasting response dataset to the corresponding grid of the three-dimensional finite element model.

[0109] Apply the target explosion load to the established three-dimensional finite element model. When the blasthole wall pressure time history needs to be determined, the specific formula for calculating the blasthole wall pressure time history is:

[0110]

[0111] In the formula, is the blasthole wall pressure at time t, which is the pressure that the blasthole wall receives at time t due to blasting, is the peak pressure of blasting, is the attenuation coefficient, is the blast pulse width. When the explosion occurs, the blasthole wall is subjected to an instantaneous pressure wave, and this pressure wave propagation has the characteristic of rapid change. The pressure generated by the explosion is not static, but rapidly reaches a peak and then decays with time. The peak pressure of the explosion is determined by the properties of the explosive, the charge amount, the geometry of the blasthole, and the characteristics of the surrounding medium. Experimental data show that at the moment of explosion, the blasthole wall will be subjected to a relatively high pressure, which can be obtained through theoretical calculation.

[0112] The attenuation coefficient represents the rate of pressure attenuation, which is usually related to the physical properties of the material, such as medium density, elastic modulus, etc., and the explosion conditions. By consulting relevant research literature in the field, reference values of attenuation coefficients under similar conditions can be obtained, generally in the range of 0.01 to .

[0113] where the peak pressure of blasting is calculated according to the formula:

[0114]

[0115] wherein, is the explosive density, is the explosive detonation velocity; wherein the borehole wall pressure applied by the blasting at time t The simulated blasting load is input into the three-dimensional finite element model as an applied load.

[0116] The blasting pulse width The formula on which the calculation is based is:

[0117]

[0118] wherein, is the charge diameter.

[0119] The smooth particle hydrodynamics method is used to simulate the blasting impact to generate a blasting response dataset, which includes the stress tensor and velocity vector of the rock element in each grid in the finite element model, wherein the formula on which the calculation of the velocity vector of the grid-rock element is based is:

[0120]

[0121] wherein is the velocity vector of the cth rock element in the grid at time t, is the mass of the bth rock element in the neighborhood of the cth rock element, is the density of the cth rock element, is the density of the bth rock element in the neighborhood of the cth rock element, is the artificial viscosity term of the cth rock element and the bth rock element in the neighborhood thereof, is the position gradient of the cth rock element, is the kernel function, and are the Cauchy stresses borne by the cth rock element and the bth rock element in the neighborhood thereof at time t, respectively, wherein b is the index of the rock element in the neighborhood of the cth rock element, wherein is the total number of rock elements in the neighborhood of the cth rock element, and c is the index of the rock element in the grid, and the neighborhood of the cth rock element is the set of rock elements adjacent to the cth rock element in the same grid; and the velocity vector of the entire grid is taken as the velocity vector of all the rock elements in the grid.

[0122] Where the viscosity term is a numerical stability mechanism introduced to handle high gradient, non-physical oscillations, etc. problems caused by blasting or shock waves in numerical simulation, which is calculated through relative velocity, relative position, density, sound speed and smoothing length, etc. parameters, and is used to handle high gradient problems such as shock waves.

[0123] Position gradient is the gradient of the kernel function with respect to the coordinates, which can be obtained by taking the derivative of the analytical expression of the kernel function used.

[0124] Kernel function is the basis of smoothed particle hydrodynamics, which is used to smooth the discrete data, and its specific form can be selected according to the type of problem, such as Gaussian kernel function or cubic spline kernel function.

[0125] Where the rock element is a representative of the material points in the grid element. Each grid element can be divided into multiple rock elements, the specific steps are as follows: define the volume of each element, the volume of the element usually depends on the calculation accuracy and simulation requirements, which can be set as a small component of the grid element, and the rock elements are evenly divided by the total volume of each grid element, and the rock element containing the grid center point coordinates is taken as the target element for analysis.

[0126] Smoothed particle hydrodynamics is a meshless numerical method suitable for simulating the dynamic behavior of fluids and solids, especially in the face of complex boundaries and large deformations. This method uses particles or elements to represent materials, and calculates forces and movements through interactions. The formula reflects the local force balance relationship. The rate of change of velocity is determined by the mass, stress and density of adjacent elements, reflecting the force transmission and action mechanism.

[0127] At the same time, Cauchy stress is the main way to describe the internal force of the material, which can fully reflect the stress state of the material. By combining stress and density, the dynamic response of the material can be more accurately described.

[0128] Where the Cauchy stress on the cth rock element at time t is Specifically calculated by the damage variable, the specific formula is:

[0129]

[0130] In the formula, is the damage variable of the cth rock element at time t, is the elastic stiffness tensor of the cth rock element, is the elastic strain tensor of the cth rock element, which is characterized by finite element analysis, where the vector sum of the Cauchy stress of all rock elements in the grid is taken as the stress tensor of the corresponding grid.

[0131] It is noted that this formula is based on the theory of damage mechanics, which is used to describe and analyze the changes in macroscopic mechanical properties of materials as the microscopic structure within the material gradually deteriorates under increasing loads during loading. By introducing a damage variable, the stiffness and strength of the material can be dynamically adjusted to reflect its actual performance.

[0132] The part of the formula is a traditional elastic mechanics formula, following Hooke's law. In the undamaged state, the stress and strain of the material are linearly related. Through the modification of the damage variable, the behavior of the material after damage can be modeled.

[0133] Step 4: Based on the damage variable evolution process output by the elastoplastic constitutive model, the rock permeability is dynamically updated and the gas seepage control equation is solved. The obtained gas pressure is applied as an additional load on the surface of the damage element to realize the coupling effect of blasting impact and gas driving, and complete the simulation of tunnel blasting impact disturbance load.

[0134] The formula used to dynamically update the rock permeability is as follows:

[0135]

[0136] where, is the rock permeability of the a-th grid at time t, is the initial rock permeability of the a-th grid, is the damage-permeability correlation index;

[0137] It is noted that this formula describes the dynamic updating of rock permeability during damage evolution. Through the change of the damage variable, the permeability also changes, thereby reflecting the seepage characteristics of the rock mass under the action of stress and damage. It is the inherent property of rock mass material and is usually measured by experiment.

[0138] In rock mass material, damage can lead to the generation and propagation of microcracks. The formation of these microcracks not only affects the bearing capacity of the rock mass, but also changes its pore structure, thereby affecting the permeability. Damage increases usually leads to an increase in permeability, because microcracks provide more flow channels, so the rock permeability of the a-th grid at time t is proportional to , and The introduction of the ratio is used to describe the proportional relationship between the damaged and undamaged parts. As the damage increases, the ratio gradually increases, meaning that the effect of damage on permeability is more significant.

[0139] where the damage-permeability correlation index Specifically, through laboratory tests, permeability measurements are conducted on rock samples with different damage levels. Common permeability test methods can be used, such as conventional water permeability tests and gas permeability tests. By applying different stress levels or cyclic loading to the samples, damage is induced, and permeability is measured at each stage. Record the damage level, such as crack propagation, porosity change, etc., and correlate it with the change in permeability to obtain the damage-permeability correlation index ; or in geotechnical engineering and materials science, many documents and researches have provided empirical relationships between damage level and permeability. These empirical relationships are usually based on regression analysis of a large amount of experimental data, which can provide a reference for the selection of ; generally 0.1 to 0.2, due to the simulation analysis of the grid area of the tunnel, as the length, width and height of the tunnel increase, the number of grids will increase with the volume of the entire area, and generally the number of grids for tunnel analysis will reach millions, which cannot be fully displayed, therefore only part of the grid damage simulation quantity is shown here, which is used to reflect the test analysis basis; Table 1 shows the rock permeability calculation statistical data of part of the grid.

[0140] Table 1: Rock permeability calculation statistics table of part of the grid

[0141]

[0142] In this data analysis, the compressive strength of rock mass in different grid areas, the total number of fractures and their related parameters are systematically evaluated. The sample data table includes the compressive strength, fracture number, pre-blasting damage state, process damage variable and rock permeability of each detection area. The setting of these data aims to reflect the mechanical properties and stability of the rock mass.

[0143] By observing the compressive strength of each sample, it can be found that the compressive strength of most samples is between 20.5 MPa and 30.0 MPa, showing that the compressive capacity of the overall rock mass is diversified. For example, the compressive strength of sample 10 is 30.0 MPa, although the total number of fractures is 2, but its pre-blasting damage state is low, only 0.0089, indicating that the sample performs well under external load.

[0144] The analysis of the process damage variable shows that the process damage variable of most areas in the sample is low, especially sample 1 and sample 4, which are 0.011 and 0.016 respectively. This reflects that the corresponding area has good stability, and regular monitoring of such areas is needed to ensure their compressive capacity. The process damage variable of sample 5 is 0.107, although its compressive strength is 21.9 MPa, but the high damage state suggests that this area needs attention.

[0145] The formula used to solve the gas seepage control equation is:

[0146]

[0147] In the formula, is the gas seepage velocity vector of the a-th grid at t, is the gas dynamic viscosity, is the gas pressure suffered by the a-th grid at t, is the gas pressure gradient at the a-th grid, is the porosity of the a-th grid, is the gas density, is the gas molar mass, is the universal gas constant, is the gas temperature, is the divergence operator;

[0148] The formula comprehensively considers the basic physical mechanism of gas seepage in porous media, and provides a comprehensive framework for describing gas flow by combining Darcy's law, mass conservation and the ideal gas state equation. The introduction of dynamically updated permeability and porosity enables the model to reflect the actual behavior of the rock mass during the damage process, and is suitable for complex geological engineering and gas seepage research.

[0149] By solving the gas seepage control equation set, the gas pressure suffered by the a-th grid at t is obtained, and the obtained gas pressure is applied as an additional load to the surface of the damage element, thereby realizing the coupling of blasting shock and gas driving, wherein the damage element refers to the surface grid element of the tunnel region to be simulated.

[0150] The coupling of blasting shock and gas driving is realized, wherein the specific coupling mode of blasting shock and gas driving is that the stress tensor, velocity vector and gas pressure in each grid are coupled in the numerical model, specifically, the stress tensor acts as an instantaneous load on the damage element, and the gas pressure provides a background load, and the two act together on the dynamic load response of the grid element. The stress tensor velocity vector and the gradient direction of the gas pressure represent the stress direction under the dynamic load response, and the dynamic response characteristics of the material under the action of the dynamic load response are accurately predicted by using numerical calculation methods such as finite element analysis (FEM).

[0151] Please refer to Figure 6 The application also provides a tunnel blasting shock disturbance load precise simulation device, which is used for executing the tunnel blasting shock disturbance load precise simulation method.

[0152] The tunnel analysis model establishing module is configured to acquire tunnel rock mass physical data and joint fissure geometric characteristics of a tunnel region to be simulated, establish a three-dimensional finite element model of the tunnel region to be simulated based on the tunnel rock mass physical data, and implement local grid encryption based on the joint fissure geometric characteristics in a blast hole influence region;

[0153] The blasting response data analysis module is configured to apply blast hole wall pressure time history data to the established three-dimensional finite element model, simulate equivalent plastic strain and plastic strain rate time evolution data of each grid, and define an elastoplastic constitutive model of damage evolution based on the equivalent plastic strain and plastic strain rate time evolution data.

[0154] Based on the elastoplastic constitutive model, the damage of the grid rock mass is simulated, and the blasting impact is characterized by combining the smoothed particle hydrodynamics method to generate a blasting response data set, and data in the blasting response data set is applied to corresponding grids of the three-dimensional finite element model.

[0155] The blasting gas additional load simulation module is configured to dynamically update the rock mass permeability based on the damage variable evolution process output by the elastoplastic constitutive model and solve a gas seepage control equation, apply the obtained gas pressure as an additional load to the surface of the damage unit, realize the coupling effect of blasting impact and gas driving, and complete the simulation of the tunnel blasting impact disturbance load.

[0156] The above formulas are all de-dimensioned to calculate their numerical values. The formulas are obtained by collecting a large amount of data to simulate a formula closest to the actual situation. The preset parameters in the formulas are set by a person skilled in the art according to the actual situation.

[0157] The above embodiments can be realized wholly or partially by software, hardware, firmware or any other combination. When realized by software, the above embodiments can be realized wholly or partially in the form of a computer program product. Those skilled in the art can realize that the units and algorithm steps of the examples described in connection with the embodiments disclosed herein can be realized by electronic hardware or a combination of computer software and electronic hardware. Whether the functions are realized by hardware or software methods depends on the specific application and design constraints of the technical solutions.

[0158] The units described as separate components can or can not be physically separated, and the components shown as units can or can not be physical units, which can be located in one place or distributed on multiple network units. Part or all of the units can be selected according to actual needs to achieve the purpose of the embodiments.

[0159] The above merely provides the specific implementation of the present application, but the protection scope of the present application is not limited to this. Any person skilled in the art can easily think of the changes or replacements within the technical range disclosed by the present application, which should be covered in the protection scope of the present application.

Claims

1. A method for accurately simulating tunnel blasting impact disturbance loads, characterized in that, The specific steps include: Obtain the tunnel rock mass physical data and joint and fracture geometric characteristics of the tunnel area to be simulated. Based on the tunnel rock mass physical data, establish a three-dimensional finite element model of the tunnel area to be simulated. In the area affected by the blast holes, implement local mesh refinement based on the joint and fracture geometric characteristics. Time history data of borehole wall pressure were applied to the established three-dimensional finite element model to simulate the time evolution data of equivalent plastic strain and plastic strain rate of each mesh. Based on the time evolution data of equivalent plastic strain and plastic strain rate, an elastoplastic constitutive model of damage evolution was defined. Damage to meshed rock mass is simulated based on an elastoplastic constitutive model, and the blasting impact is characterized by a smooth particle hydrodynamic method. A blasting response dataset is generated, and the data in the blasting response dataset is applied to the corresponding mesh of the three-dimensional finite element model. Based on the damage variable evolution process output by the elastoplastic constitutive model, the rock permeability is dynamically updated and the gas seepage control equation is solved. The obtained gas pressure is applied as an additional load to the surface of the damaged element to realize the coupling effect of blasting impact and gas drive, and to complete the simulation of tunnel blasting impact disturbance load. In the borehole-affected area, the logic for implementing local mesh refinement based on the geometric characteristics of joints and fissures is as follows: determine the borehole-affected radius, and determine the borehole-affected area based on this radius. The borehole-affected area is a circular region divided by taking the borehole center as the center and the borehole-affected radius as the radius. The mesh within the area affected by the blast hole is refined by reducing the size of the mesh within that area. The specific formula used for mesh refinement is as follows: In the formula, The size of the grid area within the influence zone of the borehole. This is the initial grid area size. , and These represent the average length, average width, and average depth of the fracture within the borehole's influence area, respectively. This is a reference value for the fracture volume. The number of cracks within the area affected by the borehole. This is a reference value for the number of cracks. and These are the weighting coefficients for fracture volume and fracture number, respectively. and and All are greater than 0; The specific formula upon which the dynamic updating of rock mass permeability is based is: In the formula, Let be the rock permeability of the a-th grid at time t. Let be the initial rock mass permeability of the a-th grid. This is a damage-penetration correlation index. Let be the damage variable of the a-th grid at time t; The specific formula used to solve the gas seepage control equation is as follows: In the formula, Let be the gas seepage velocity vector at the a-th grid point at time t. For gas dynamic viscosity, Let be the gas pressure experienced by the a-th grid at time t. Let be the gas pressure gradient at the a-th grid. Let be the porosity of the a-th grid. For gas density, For gas molar mass, This is the universal gas constant. For gas temperature, It is a divergence operator; By solving the gas seepage control equations, the gas pressure on the a-th grid at time t is obtained, and the obtained gas pressure is applied as an additional load to the surface of the damaged element to achieve the coupling effect of blasting impact and gas drive. The damaged element specifically refers to the surface grid element of the tunnel region to be simulated.

2. The method for accurately simulating tunnel blasting impact disturbance load according to claim 1, characterized in that: The physical data of the tunnel rock mass specifically include the density, elastic modulus and compressive strength of the rock mass, and the geometric characteristics of the joints and fissures include the length, width, depth and distribution location of the fissures on the rock mass; A three-dimensional finite element model of the tunnel area to be simulated is established based on the tunnel rock mass physical data. The geometric model of the tunnel is constructed using three-dimensional modeling software to determine the shape, size and depth of the tunnel. The tunnel rock mass physical data is embedded into the geometric model, and a three-dimensional finite element model is generated by the modeling software. The entire model is meshed, the mesh geometry is selected, and local mesh refinement is implemented in the area affected by the blast holes based on the geometric characteristics of joints and fissures, specifically the length, width, depth and distribution of the fissures.

3. The method for accurately simulating tunnel blasting impact disturbance load according to claim 1, characterized in that: Define an elastoplastic constitutive model for damage evolution, wherein the specific formula upon which the elastoplastic constitutive model is based is: In the formula, This represents the damage state of the a-th grid before the blast. Let be the equivalent plastic strain of the a-th grid at time t. Let be the plastic strain rate of the a-th grid at time t. and The value represents the damage sensitivity coefficient, where a is the mesh index in the three-dimensional finite element model, t is the time variable during the blasting simulation test, and the equivalent plastic strain and plastic strain rate are obtained from the established three-dimensional finite element model through finite element simulation. The tunnel rock mass is characterized by its compressive strength and fracture characteristics, specifically based on the following formula: In the formula, Let be the rock mass compressive strength of the a-th grid before blasting. The overall compressive strength of the tunnel rock mass material before blasting. The tunnel fracture influence coefficient. This represents the total number of cracks within the a-th grid. The borehole wall pressure time history data is applied to the established three-dimensional finite element model. The specific formula used to calculate the borehole wall pressure time history is as follows: In the formula, Let be the pressure on the borehole wall at time t, which is the pressure exerted on the borehole wall by the blast at time t. For the peak blast pressure, The attenuation coefficient is... The burst pulse width, where the peak burst pressure is... The specific formula used for the calculation is as follows: In the formula, For the density of the explosive, For the detonation velocity of the explosive; Explosive pulse width The specific formula used for the calculation is as follows: In the formula, This refers to the diameter of the propellant charge.

4. The method for accurately simulating tunnel blasting impact disturbance load according to claim 3, characterized in that: A smoothed particle hydrodynamics method is used to simulate blasting impact, generating a blasting response dataset. This dataset includes the stress tensor and velocity vector of each rock mass micro-element within the finite element model. The specific formula used to calculate the velocity vector of a rock mass micro-element within a given mesh is as follows: In the formula Let be the velocity vector of the c-th rock mass element within the grid at time t. Let b be the mass of the b-th rock mass element within the neighborhood of the c-th rock mass element. Let c be the density of the c-th rock mass element. Let b be the density of the b-th rock mass element within the neighborhood of the c-th rock mass element. Let c be the artificial viscosity term between the c-th rock mass element and the b-th rock mass element in its neighborhood. Let c be the position gradient of the c-th rock mass micro-element. For kernel function, and Let be the Cauchy stresses experienced by the c-th rock mass element and the b-th rock mass element in their neighborhood at time t, respectively, where b is the index of the rock mass element in the neighborhood of the c-th rock mass element. ,in is the total number of rock mass micro-elements in the neighborhood of the c-th rock mass micro-element, where c is the index of the rock mass micro-element within the grid. The neighborhood of the c-th rock mass micro-element is the set of rock mass micro-elements adjacent to the c-th rock mass micro-element within the same grid. Wherein, at time t, the Cauchy stress experienced by the c-th rock mass element is... Specifically, the calculation is performed using damage variables, and the formula used is as follows: In the formula, Let be the damage variable of the c-th rock mass micro-element at time t. Let be the elastic stiffness tensor of the c-th rock mass element. Let be the elastic strain tensor of the c-th rock mass micro-element, specifically characterized by finite element analysis, where the vector sum of the Cauchy stresses experienced by all rock mass micro-elements within the grid is used as the stress tensor of the corresponding grid.

5. A device for accurately simulating tunnel blasting impact disturbance loads, characterized in that: The tunnel blasting impact disturbance load precision simulation device is used to execute the tunnel blasting impact disturbance load precision simulation method according to any one of claims 1-4, including: The tunnel analysis model building module is used to acquire the tunnel rock mass physical data and joint and fracture geometric features of the tunnel area to be simulated. Based on the tunnel rock mass physical data, a three-dimensional finite element model of the tunnel area to be simulated is built, and local mesh refinement is implemented in the area affected by the blast holes based on the joint and fracture geometric features. The blast response data analysis module is used to apply borehole wall pressure time history data to the established three-dimensional finite element model, simulate the time evolution data of equivalent plastic strain and plastic strain rate of each grid, and define the elastoplastic constitutive model of damage evolution based on the time evolution data of equivalent plastic strain and plastic strain rate. Damage to meshed rock mass is simulated based on an elastoplastic constitutive model, and the blasting impact is characterized by a smooth particle hydrodynamic method. A blasting response dataset is generated, and the data in the blasting response dataset is applied to the corresponding mesh of the three-dimensional finite element model. The blasting gas additional load simulation module is used to dynamically update the rock permeability and solve the gas seepage control equation based on the damage variable evolution process output by the elastoplastic constitutive model. The obtained gas pressure is applied as an additional load to the surface of the damaged element to realize the coupling effect of blasting impact and gas drive, and complete the simulation of tunnel blasting impact disturbance load.

Images