Blasting process simulation method and system based on bidirectional fluid-solid coupling

Abstract

The invention discloses a blasting process simulation method and system based on bidirectional fluid-solid coupling, and relates to the technical field of rock blasting. The method comprises the following steps: firstly, modeling a solid computational domain by adopting a discrete element method, and modeling a fluid computational domain by adopting a finite volume method; and then applying a dynamic load to simulate an initial crack in a shock wave stage, mapping crack information to a fluid computational domain grid, dynamically updating local permeability, solving a compressible gas diffusion equation, applying a pressure field to the surface of a solid computational domain particle through an interpolation method, and driving the crack to continuously expand. A fluid-solid two-way feedback mechanism covering detonation gas and rock is established, through closed-loop coupling of crack-permeability-gas pressure-new crack, the driving effect of gas at the tip of the crack is accurately captured, the experimental crack form can be accurately reproduced, and the experiment precision is improved. And the influence of the initial air pressure on the crushing range, the crack number and the damage mode can be quantitatively analyzed, and the method is suitable for blasting optimization design under the condition of high-gas-yield explosives.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of rock blasting mechanics and numerical simulation technology, and particularly relates to a blasting process simulation method and system based on two-way fluid-structure coupling, which is suitable for simulating the wedging, loading and feedback behavior of high-pressure blast gas in the dynamic evolution of the fracture network in the blasting process. BACKGROUND

[0002] In the process of rock blasting, although the stress wave at the initial stage of blasting can induce dense micro-cracks, the blast gas in the subsequent millisecond quasi-static stage is also the key factor for the crack penetration, continuous energy injection and final fragmentation of the rock mass. However, the blast gas has complex characteristics such as instantaneous high pressure, strong nonlinear permeation and dynamic interaction with the fracture network, and its real-time diffusion path, local pressure distribution and driving mechanism for the crack tip in the rock mass are difficult to be directly observed through experiments, and therefore high-fidelity numerical means are urgently needed to reveal them.

[0003] At present, the blasting numerical simulation mainly adopts the finite element method (FEM), the extended finite element method (XFEM) or the discrete element method (DEM). Among them, FEM / XFEM can better describe the deformation of continuous medium, but it is difficult to spontaneously simulate the initiation, branching and penetration of complex cracks; DEM can naturally reproduce the whole process of rock failure, but the existing loading methods mostly use the method of applying decaying pressure on the blast hole wall (equivalent pressure method) or the method of expanding the center particles, ignoring the real migration process of the gas along the fracture network, resulting in no pressure loading in the fracture, and failing to reflect the gas accumulation effect at the crack tip and the dynamic correlation between gas diffusion and permeability evolution.

[0004] Although the existing methods have promoted the refinement of blasting simulation to some extent, they still have significant limitations. For example, Chinese patent CN114218831B proposes to simulate crack propagation by pre-setting joint / interface elements between solid elements, which can avoid mesh redrawing and handle complex fracture patterns, but the fracture path still depends on prior settings, and no gas diffusion or flow physical model is introduced, so the dynamic driving effect of blast gas on crack propagation cannot be reflected. For another example, Chinese patent CN115859714A constructs an FEM-DEM joint simulation framework, and applies the blast pressure time history calculated by FEM as boundary load to the DEM model, which can take into account the effects of shock wave and quasi-static pressure, but the gas action is strictly limited to the initial blast hole wall or the crushing zone, and the real-time migration process of the gas to the peripheral new fractures is not tracked, and the local pressure accumulation effect of the gas at the crack tip is not considered. For another example, Chinese patent CN109631701B focuses on the macro load superposition and vibration prediction of multi-blast hole blasting, and improves the simulation accuracy of vibration velocity by applying the blasting stress to the tunnel contour surface, but it does not involve the diffusion behavior of blast gas in the fracture network of rock mass, and lacks physical description of the millisecond quasi-static fragmentation evolution mechanism.

[0005] It can be seen that the existing method generally does not introduce the permeation flow mechanism of the detonation gas in the fracture network, and usually simplifies the gas action as a preset pressure time history on the blast hole wall, completely ignoring the physical nature of the discrete fracture as the main channel for gas diffusion. In fact, the expansion of the crack will open a new high-conductivity channel, and the gas pressure will drive the crack to further expand in the opposite direction, forming a strong nonlinear feedback loop. Ignoring this two-way coupling mechanism will lead to serious distortion of the blasting and crushing evolution process.

[0006] Therefore, there is an urgent need for a numerical simulation device and method capable of synchronously, bidirectionally and dynamically coupling gas diffusion and rock mass rupture, which can truly reproduce the whole process of migration, loading and feedback of detonation gas in evolving fractures, and provide scientific support for blasting mechanism cognition and engineering parameter optimization. SUMMARY

[0007] The purpose of the present application is to provide a blasting process simulation method based on bidirectional fluid-structure coupling. The present application integrates the finite volume method (FVM) and the discrete element method (DEM), realizes bidirectional data interaction by using a self-defined interface, constructs a closed-loop feedback mechanism of "crack identification-permeability update-gas solving-pressure mapping-particle loading-new crack generation", and truly simulates the diffusion and driving behavior of detonation gas in a dynamic fracture network.

[0008] The present application adopts a staged loading strategy: first, the initial cracks in the shock wave stage are simulated by applying a dynamic load based on the Duvall pressure time history; then, the fluid-structure coupling module is started, the crack information is mapped to the fluid grid, the local permeability is dynamically updated, the compressible gas diffusion equation is solved, and the pressure field is applied to the particle surface through the interpolation method to drive the continuous expansion of the fracture. This method not only can accurately reproduce the crack morphology of the experiment, but also can quantitatively analyze the influence of the initial gas pressure on the crushing range, the number of cracks and the damage mode, and is especially suitable for blasting optimization design under the condition of high gas yield explosives.

[0009] To achieve the above purpose, the technical scheme adopted by the present application is:

[0010] A blasting process simulation method based on bidirectional fluid-structure coupling, comprising:

[0011] S1: using the discrete element method (DEM) to model the solid calculation domain, and using densely arranged spherical particles to represent the rock medium; using the finite volume method (FVM) to model the fluid calculation domain, and using a grid area consistent in size with the solid calculation domain to describe the diffusion behavior of the gas therein;

[0012] S2: Based on the solid computing domain established in S1, the transient dynamic response and dynamic fracture evolution process under the action of blasting shock wave are solved, and the initial crack information of rock in the shock wave stage is extracted synchronously, including the initiation, propagation path and spatial distribution of cracks;

[0013] S3: Based on the fluid computing domain constructed in S1 and the crack information generated by the solid computing domain in S2, the non-steady seepage process of blast gas in the crack is solved, and the spatial distribution of gas pressure field in the fluid computing domain at the current time step is obtained;

[0014] S4: Map the gas pressure field in the fluid computing domain at the current time step to the solid computing domain, and apply it as an additional force on the particles in the solid computing domain to drive the initiation and propagation of new cracks;

[0015] S5: Extract the updated crack information and repeat steps S3-S4 to realize the bidirectional coupling between the fluid computing domain and the solid computing domain through closed-loop iteration, and drive the system to evolve to the next time step until the gas pressure decays to insufficient to drive new cracks.

[0016] Further, in step S1, the process of modeling the solid computing domain and modeling the fluid computing domain includes:

[0017] S1.1: In the discrete element method (DEM) platform, spherical particles are used to model the rock sample, the spherical particles are connected by parallel bond model (PBM), and key mechanical parameters such as elastic modulus, tensile strength, cohesion, friction angle, etc. are defined independently;

[0018] S1.2: In the finite volume method (FVM) solver, a grid identical in size to the solid computing domain is established as the fluid computing domain;

[0019] S1.3: In the solid computing domain, the blast hole area is modeled as a cavity formed by removing the center particles, which is used to apply the transient stress load in the blasting shock wave stage; In the fluid computing domain, the blast hole area is set as the inlet boundary condition of blast gas to simulate the injection of blast gas; The outer boundary condition of the solid computing domain is set as free boundary or fixed constraint according to the actual working condition; The outer boundary condition of the fluid computing domain is set as free pressure outlet boundary to allow gas to flow out freely, to simulate the gas diffusion behavior in an open environment;

[0020] S1.4: Initialize the porosity and permeability of the solid computing domain as the reference diffusion parameters under the unbroken state.

[0021] Further, in step S2, the specific process of extracting the initial crack information of rock in the shock wave stage is as follows:

[0022] S2.1: In the shock wave stage, a dynamic pressure boundary condition varying with time is applied to the borehole wall surface of the solid calculation domain, and the pressure time history is described by the Duvall formula, and the specific expression is:

[0023]

[0024] wherein, is the pressure at time is the dynamic stress applied to the borehole wall surface of the solid calculation domain; and α and β are frequency-dependent attenuation constants; is a normalization factor for ensuring that the pressure curve is reasonably normalized in the rising segment; denotes the peak pressure;

[0025] S2.2: Run the DEM solver to calculate the shock wave propagation and particle motion process;

[0026] S2.3: Record all the cohesive bonds in the solid calculation domain that are broken due to tensile or shear failure, and output the spatial position and direction of the cohesive bonds as the initial fracture network.

[0027] Further, the peak pressure is determined by the type of explosive, and the expression is:

[0028]

[0029] wherein, denotes the rise time of the shock pressure, which is obtained by the following formula:

[0030]

[0031] The attenuation constants α and β are determined by the following empirical relationship:

[0032]

[0033] wherein, is a pressure attenuation characteristic parameter; denotes the longitudinal wave velocity of the rock, with the unit of m / s; is the borehole diameter, with the unit of m.

[0034] Further, in step S3, solving the non-steady seepage process of the detonation gas specifically includes:

[0035] S3.1: Map the fracture spatial coordinates obtained in step S2 to the grid cells of the fluid calculation domain, and identify the affected grid regions in the fluid calculation domain;

[0036] S3.2: In each grid cell of the fluid calculation domain where the fracture is located, locally increase the porosity , and update the permeability according to the power law relationship:

[0037]

[0038] wherein, is the initial permeability, is the initial porosity, n is an empirical exponent; is the porosity, defined as:

[0039]

[0040] wherein, is the number of fractures in a single grid cell of the fluid computational domain, is the porosity increment of each fracture;

[0041] S3.3: Take the residual stress size at the borehole wall at the end of the shock wave action in step S2 as the initial gas pressure of the fluid calculation domain , and take it as the initial pressure condition of the inlet boundary of the fluid calculation domain;

[0042] S3.4: Solve the unsteady seepage process of the blast gas in the fracture medium, and the control equation is composed of the mass conservation equation and the compressible form of Darcy's law; the gas mass conservation equation is expressed as:

[0043]

[0044] wherein, ρ is the gas density, unit kg / m³; is the porosity; v is the Darcy velocity, volume flow rate per unit area, unit m / s; the gas flow obeys the modified Darcy's law:

[0045]

[0046] wherein, is the permeability, unit m²; μ is the gas viscosity, unit MPa·s; is the gas pressure, the initial gas pressure is , unit MPa;

[0047] S3.5: Set the time step of the fluid calculation domain as , and obtain the spatial distribution of the gas pressure field in the fluid calculation domain at the current time step.

[0048] Further, the driving of the initiation and expansion of new cracks in step S4 specifically includes:

[0049] S4.1: Smoothly map the pressure value at the current time step on the grid cell of the fluid calculation domain to the center position of each particle in the solid calculation domain by using the interpolation method, to obtain the equivalent gas pressure suffered by the particle;

[0050] S4.2: According to the time step The solid calculation domain is calculated by using the DEM solver, new cracks are generated, and spatial information thereof is recorded.

[0051] Further, in step S5, new crack information is extracted and returned to step S3 for gas flow solving at the next time step; steps S3-S4 are repeated to form a closed loop feedback of "crack-permeability-gas pressure distribution-new crack", and the calculation is ended until the gas pressure decays to be insufficient to drive new cracks.

[0052] A blasting process simulation system based on two-way fluid-solid coupling, comprising:

[0053] A parameter input module is configured to input all initial conditions and material properties required for simulation, including the elastic modulus, tensile strength, cohesion, internal friction angle of the solid calculation domain, the initial gas pressure, gas viscosity, gas density, initial porosity, initial permeability, and power law index of permeability evolution with porosity of the fluid calculation domain, and coupling control parameters such as time step, total simulation time, peak pressure of Duvall pressure time history load, normalization factor, and decay constant.

[0054] A granular mechanics solving module is configured to simulate the dynamic response and spontaneous crack evolution of rock under blasting load in the solid calculation domain, and extract crack information of rock in the shock wave stage;

[0055] A gas diffusion solving module is configured to map the crack information to the fluid calculation domain, solve the non-steady-state seepage process of the blast gas in the crack, and obtain the spatial distribution of the gas pressure field in the fluid calculation domain;

[0056] A fluid-solid coupling module is configured to realize real-time data exchange between the crack information obtained by the granular mechanics solving module and the gas pressure field obtained by the gas diffusion solving module, establish a fluid-solid coupling mechanism, and realize two-way coupling between the fluid calculation domain and the solid calculation domain.

[0057] A result output module is configured to output the crack network morphology, gas pressure evolution cloud map, crack propagation path, and key monitoring point response time history.

[0058] The beneficial effects of the present application are:

[0059] (1) The present application realizes the simulation of the interaction between blast gas and rock in a dynamic crack network, breaks through the conventional method of traditional hole wall loading, and significantly improves the physical authenticity;

[0060] (2) The present application establishes a fluid-solid two-way feedback mechanism covering blast gas and rock, accurately captures the driving effect of gas at the crack tip through the closed loop coupling of crack-permeability-gas pressure-new crack.

[0061] (3) The application can separate the action of the shock wave and the gas, support mechanism research, and reveal the control law of the shock wave and the gas on the crushing mode through parameterized analysis.

[0062] In summary, the application provides a blasting gas-fracture coupling simulation method with high fidelity, scalability and clear physical mechanism, solves the core problems of gas loading distortion, coupling loss and inability to reflect the real diffusion path in the prior art, not only deepens the understanding of the blasting crushing mechanism, but also provides a new research means for the precise control and safety optimization of mine, tunnel and other engineering blasting. BRIEF DESCRIPTION OF DRAWINGS

[0063] Figure 1 The simulation method flowchart of the application.

[0064] Figure 2 Modeling of the flow and solid calculation domain.

[0065] Figure 3 Separating Hopkinson pressure bar (SHPB) dynamic mechanical parameter calibration.

[0066] Figure 4 Shock wave stage stress curve loading.

[0067] Figure 5 Comparison of experiments and simulation methods without gas action; wherein (a) is the stress distribution of the simulation blasting process of the application, and (b) is the strain distribution of the blasting test.

[0068] Figure 6 Blasting simulation whole process under the action of detonation gas; wherein (a) is the stress distribution, (b) is the permeability distribution, and (c) is the gas stress distribution.

[0069] Figure 7 Stress condition comparison in the simulation of the embodiment of the application; wherein (a) is the stress curve without gas action, and (b) is the stress curve under the action of detonation gas.

[0070] Figure 8 Damage mode comparison; wherein (a) is the simulation result, and (b) is the experimental result. DETAILED DESCRIPTION

[0071] The technical solutions in the embodiments of the application will be described clearly and completely below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, not all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the application.

[0072] AsFigure 1 As shown, the present application provides a blasting process simulation method based on two-way fluid-structure coupling, comprising the following steps:

[0073] S1: Discrete element method (DEM) is used to model the solid calculation domain; finite volume method (FVM) is used to model the fluid calculation domain; the specific process is as follows:

[0074] S1.1: In the discrete element (DEM) platform, spherical particles are used to model the rock sample, the spherical particles are connected through parallel bond model (PBM), and the key mechanical parameters such as elastic modulus, tensile strength, cohesion, friction angle, etc. are independently defined;

[0075] S1.2: In the finite volume method (FVM) solver, a grid completely consistent with the geometric size of the solid calculation domain is established as the fluid calculation domain;

[0076] S1.3: In the solid calculation domain, the blast hole area is modeled as a cavity formed by removing the center particle, which is used to apply the transient stress load in the blast shock wave stage; in the fluid calculation domain, the blast hole area is set as the inlet boundary condition of the blast gas to simulate the injection of the blast gas; the outer boundary condition of the solid calculation domain is set as free boundary or fixed constraint according to the actual working condition, and the outer boundary condition of the fluid calculation domain is set as free pressure outlet boundary to allow the gas to flow out freely, thereby simulating the gas diffusion behavior in open environment;

[0077] S1.4: Initialize the porosity and permeability of the solid calculation domain as the reference diffusion parameters in the unbroken state.

[0078] S2: Based on the solid calculation domain established in S1, the transient dynamic response and dynamic fracture evolution process under the action of blast shock wave are solved, and the initial crack information of the rock in the shock wave stage is extracted synchronously, including the initiation, expansion path and spatial distribution of the crack; the specific process is as follows:

[0079] S2.1: In the shock wave stage, a dynamic pressure boundary condition varying with time is applied to the blast hole wall surface of the solid calculation domain, and the pressure time history is described by Duvall formula, and the specific expression is:

[0080]

[0081] Wherein, is the dynamic stress applied to the blast hole wall surface of the solid calculation domain at time α and β are frequency-dependent attenuation constants; represents the peak pressure, which is determined by the type of explosive and can be represented as:

[0082]

[0083] where, p e is the density of explosive, unit g / cm³; D is the detonation velocity, unit m / s;

[0084] is a normalization factor, which is used to ensure the pressure curve is reasonably normalized in the rising section, and its expression is:

[0085]

[0086] where, represents the rising time of the impact pressure, which can be obtained by the following formula:

[0087]

[0088] The decay constants a and b are determined by the following empirical relationship:

[0089]

[0090] where, is the pressure decay characteristic parameter; represents the longitudinal wave velocity of rock, unit m / s; is the borehole diameter, unit m;

[0091] S2.2: Run the DEM solver to calculate the shock wave propagation and particle motion process;

[0092] S2.3: Record all the cohesive bonds in the solid calculation domain that are broken due to tensile or shear failure, and output the spatial position and direction of the cohesive bonds as the initial fracture network.

[0093] S3: Based on the fluid calculation domain constructed in S1, combined with the fracture information generated by the solid calculation domain in S2, solve the non-steady seepage process of the detonation gas, and obtain the spatial distribution of the gas pressure field in the fluid calculation domain at the current time step; Specifically as follows:

[0094] S3.1: Map the fracture spatial coordinates obtained in step S2 to the grid elements of the fluid calculation domain, and identify the affected grid regions in the fluid calculation domain;

[0095] S3.2: In each grid element of the fluid calculation domain where the fracture is located, locally increase the porosity , and update the permeability according to the power law relationship:

[0096]

[0097] where, is the initial permeability, φ0is the initial porosity, n is an empirical exponent; φ is the porosity, defined as:

[0098]

[0099] where, Nfis the number of fractures in a single grid cell of the fluid computational domain, φfis the porosity increment of each fracture;

[0100] S3.3: The residual stress magnitude at the borehole wall at the end of the shock wave action in step S2 is taken as the initial gas pressure of the fluid computational domain , and is taken as the initial pressure condition of the inlet boundary of the fluid computational domain;

[0101] S3.4: The unsteady seepage process of the blast gas in the fracture medium is solved, and the control equation is composed of the mass conservation equation and the compressible form of Darcy's law; the gas mass conservation equation is expressed as:

[0102]

[0103] where, ρ is the gas density, unit kg / m³; φ is the porosity; v is the Darcy velocity, volume flow rate per unit area, unit m / s; the gas flow obeys the modified Darcy's law:

[0104]

[0105] where, k is the permeability, unit m²; μ is the gas viscosity, unit MPa·s; P is the gas pressure, and the initial gas pressure is , unit MPa;

[0106] S3.5: The time step of the fluid computational domain is set as , and the spatial distribution of the gas pressure field in the fluid computational domain at the current time step is obtained.

[0107] S4: The gas pressure field at the current time step in the fluid computational domain is mapped to the solid computational domain, and is applied as an additional force on the particles in the solid computational domain, driving the initiation and expansion of new fractures; specifically:

[0108] S4.1: An interpolation method is used to smoothly map the pressure value P at the current time step on the grid cells of the fluid computational domain to the center position of each particle in the solid computational domain, to obtain the equivalent gas pressure on the particle;

[0109] S4.2: According to the time step , the solid domain is calculated using the DEM solver to generate new fractures and record their spatial information.

[0110] S5: extract new crack information and return to step S3 for gas flow solving at next time step; repeat steps S3-S4 to form a closed loop feedback of "crack-permeability-gas pressure distribution-new crack" until gas pressure decays to insufficient to drive new crack, and the calculation ends.

[0111] Embodiment

[0112] In the polycarbonate (PC) plate blasting test conducted by Yang et al. (Visualizing the blast-induced stress wave and blasting gas action effects using digital image correlation), the specimen has a geometry of 300 mm × 300 mm × 5 mm, a center hole diameter of 4 mm, and lead azide explosive is used for blasting. The dynamic elastic modulus of the PC material is 4.5 GPa, and the dynamic Poisson's ratio is 0.32.

[0113] To verify the effectiveness of the method of the present application, numerical simulation is carried out for the above test case. To systematically evaluate the effects of the shock wave and the blast gas on rock mass fracture, two groups of comparative working conditions are set: (1) shock wave effect - only the transient stress load in the shock wave stage is applied, the subsequent gas pressure is not considered, and only the solid calculation domain is activated for solving; (2) combined effect - both the shock wave stress and the permeation pressure of the blast gas in the crack are considered, and both the solid calculation domain and the fluid domain are activated for solving. By comparing the results of the two groups, the key contribution of the blast gas in the crack propagation and penetration process can be clearly revealed.

[0114] The Duvall pressure time history load is used in the shock wave stage, and the relevant parameters are as follows: peak pressure P m = 293 MPa, normalized factor ξ = 3.47, decay constant α = 0.0362 μs -1 , β = 0.0808 μs -1 . The numerical simulation time step is set to 1×10 -9 s, and the total time is 100 μs, so as to accurately capture the initiation process of the initial microcrack under the action of the shock wave. The shock wave stress time history curve is shown in Figure 4 . After entering the blast gas action stage, the initial gas pressure = 27 MPa (corresponding to the residual pressure at the end of the shock wave decay, see Figure 4 ) is taken as the inlet boundary condition of the fluid domain; the gas viscosity is taken as μ = 3.0×10 -5 ​Pa·s, initial density ρ0 = 40 kg / m³, and dynamic coupling of density and pressure is achieved using the ideal gas isothermal assumption; the initial porosity of the solid matrix is ​​set to = 0.02, initial penetration rate = 1×10 -12 m², the porosity increment caused by cracks is taken as = 0.005, penetration rate evolution follows a power law relationship. (Exponent n = 4).

[0115] First, a solid-state computational domain model is performed. (See attached image) Figure 2 As shown, a plate-shaped sample with dimensions of 300 mm × 300 mm × 5 mm was constructed in the solid computational domain using the discrete element method (DEM). It contained 34,547 spherical particles and a cylindrical borehole with a diameter of 4 mm and a height of 5 mm in the center. The mechanical interactions between the particles were analyzed using... Figure 2 The parallel bond model (PBM) shown in (c) is used for description. Key mechanical parameters of the parallel bond model (PBM) are obtained through... Figure 3 The split Hopkinson bar (SHPB) test was used for calibration, and the final material parameters were determined as follows: elastic modulus of 3.1 GPa, tensile strength of 61 MPa, cohesion of 85 MPa, and internal friction angle of 50°. To monitor the dynamic response caused by the blast, five monitoring points were set up at 25 mm intervals along the main crack propagation direction at the center of the borehole. Free boundary conditions were applied around the model to simulate an unconstrained free surface.

[0116] Then, the fluid computational domain was modeled. The fluid computational domain was established using the finite volume method (FVM), with its geometry identical to the solid sample (300 mm × 300 mm × 5 mm), and the mesh distribution as shown... Figure 2 As shown in (b), the inner wall of the borehole is defined as the inlet boundary of the explosive gas, and the initial pressure conditions determined by the shock wave stage are applied. The outer boundary of the model is set as the pressure outlet boundary, allowing the gas to flow out freely to simulate the gas diffusion behavior in an open environment. During the simulation, the explosive gas is injected from the inlet, propagates in the dynamically evolving fracture network, and forms an unsteady pressure field, thereby achieving fluid-structure interaction with the solid computational domain.

[0117] Figure 5 This illustrates the stress distribution under stress wave action alone. Figure 5 (a) is the stress cloud diagram obtained by numerical simulation in this embodiment, with compressive stress as positive and tensile stress as negative; Figure 5The stress distribution in the elastic vibration range measured by Yang et al. in the experiment is converted according to Hooke's law, and it can be known that the stress distribution is approximately in the range of 0-12 MPa. It can be seen that the experimental results and the numerical simulation of the embodiment are in good agreement in the stress distribution form and amplitude.

[0118] Figure 6 The blasting simulation results under the combined action of stress wave and explosive gas are shown. It can be seen that under the action of gas, the distribution of compressive stress is not much different from that without gas, and the main difference is in the distribution of tensile stress. In the initial stage of gas action, the maximum tensile stress in the sample reaches about 100 MPa, at this time the permeability in the crack region increases significantly, and the gas pressure also reaches 27 MPa. Under the continuous loading of explosive gas, the crack expands from the blast hole outward, the tensile stress is concentrated at the tip, and the permeability in the crack region further increases. At the same time, due to energy dissipation and crack expansion, the gas pressure gradually decreases, and when the gas pressure is not enough to continue to drive crack expansion, the crack expansion stops.

[0119] Figure 7 The stress-time response law of the sample along the axial direction under different working conditions is shown. Among them, Figure 7 In (a), only the stress evolution curve considering the action of stress wave is shown, and it can be seen that the stress peak value is 51.37 MPa, which rapidly decays within 1.5 ms, indicating that this stage is mainly controlled by the propagation of shock wave, the duration is short, and the energy release is rapid. In contrast, Figure 7 In (b), the stress response process after coupling the action of explosive gas is shown, and the axial stress peak value rises to 83.59 MPa, which is much higher than that without gas action, and the stress decay process is significantly prolonged, and there is still a high residual stress after 6 ms.

[0120] In addition, as the measurement point position gradually moves away from the explosion source, the stress peak value shows a significant decay trend, but the residual stress under the action of gas remains at a high level, reflecting that the explosive gas continues to provide loading in the quasi-static stage, effectively prolonging the time window of crack expansion. This delayed loading behavior helps the expansion, penetration and connection of the original crack, and further promotes the destruction of the structure on a larger scale. As can be seen, the gas action not only increases the tensile stress level at the crack tip, but also enhances the driving force for crack expansion.

[0121] In addition, the rock mass damage modes obtained by the experiment of Yang et al. and the numerical simulation of the embodiment are compared, as shown in Figure 8 It can be seen that both have high consistency in crack initiation, expansion direction and damage form. The results show that the numerical simulation method proposed in the embodiment can accurately simulate the cooperative mechanism of stress wave propagation and gas fracturing, and verifies its applicability in simulating the whole blasting process.

[0122] Finally, it should be noted that the above examples are intended to illustrate the technical solutions of the present application and do not constitute any form of limitation on the present application. Those skilled in the art should fully understand that it is completely feasible to modify the technical solutions described in the foregoing examples or to make equivalent replacements on any part or all of the technical features. These modifications or replacements should be considered as reasonable extensions of the present application as long as they do not deviate from the protection scope defined by the claims of the present application.

Claims

1. A method for simulating blasting processes based on two-way fluid-structure interaction, characterized in that, include: S1: The discrete element method is used to model the solid computational domain, and the rock medium is represented by densely arranged spherical particles; The finite volume method is used to model the fluid computational domain, and a mesh region with the same geometry as the solid computational domain is used to describe the diffusion behavior of the gas in it; S2: Based on the solid computational domain established by S1, solve the transient dynamic response and dynamic fracture evolution process under the action of blasting shock wave, and simultaneously extract the initial fracture information of rock during the shock wave stage, including fracture initiation, propagation path and spatial distribution. S3: Based on the fluid computational domain constructed in S1 and the crack information generated by the rupture of the solid computational domain in S2, solve the unsteady seepage process of the explosive gas in the crack, and obtain the spatial distribution of the gas pressure field in the fluid computational domain at the current time step. S4: Map the gas pressure field of the current time step in the fluid computational domain to the solid computational domain, and apply it as an additional force to the particles in the solid computational domain to drive the initiation and propagation of new cracks. S5: Extract the updated crack information, repeat steps S3~S4, achieve bidirectional coupling between the fluid computational domain and the solid computational domain through closed-loop iteration, and evolve to the next time step until the gas pressure decays to a level insufficient to drive new cracks.

2. The method for simulating the blasting process based on bidirectional fluid-structure interaction according to claim 1, characterized in that, Step S1, the process of modeling the solid computational domain and the fluid computational domain includes: S1.1: In the discrete element platform, spherical particles are used to model the rock sample. The spherical particles are connected by a parallel bonding model, and key mechanical parameters including elastic modulus, tensile strength, cohesion, and friction angle are defined independently. S1.2: In the finite volume method solver, a mesh with the exact same geometry as the solid computational domain is established as the fluid computational domain; S1.3: In the solid computational domain, the borehole region is modeled as a cavity formed by removing the central particle, used to apply transient stress loads during the blast shock wave phase; in the fluid computational domain, the borehole region is set as the inlet boundary condition for explosive gas to simulate the injection of explosive gas; the outer boundary condition of the solid computational domain is set as a free boundary or a fixed constraint according to the actual working conditions; the outer boundary condition of the fluid computational domain is set as a free pressure outlet boundary, allowing gas to flow out freely to simulate gas diffusion behavior in an open environment; S1.4: Initialize the porosity and permeability of the solid computational domain as the baseline diffusion parameters in the unruptured state.

3. The method for simulating the blasting process based on bidirectional fluid-structure interaction according to claim 2, characterized in that, In step S2, the specific process of extracting the initial fracture information of the rock during the shock wave stage is as follows: S2.1: During the shock wave phase, a time-varying dynamic pressure boundary condition is applied to the borehole wall in the solid computational domain. The pressure time history is described using the Duvall formula, specifically: In the formula, For at any time Dynamic stress applied to the borehole wall in the solid computational domain; α and β are frequency-dependent attenuation constants; It is a normalization factor used to ensure that the pressure curve is reasonably normalized in the rising segment; Indicates peak pressure; S2.2: Run the DEM solver to calculate the shock wave propagation and particle motion process; S2.3: Record all bonded bonds that break due to tensile or shear failure in the solid computational domain, and output the spatial position and orientation of the bonded bonds as the initial fracture network.

4. The blasting process simulation method based on bidirectional fluid-structure interaction according to claim 3, characterized in that, The peak pressure is determined by the type of explosive, and is expressed as follows: In the formula, ρ e ρ is the density of the explosive, in g / cm³; D is the detonation velocity, in m / s.

5. The blasting process simulation method based on bidirectional fluid-structure interaction according to claim 3, characterized in that, The normalization factor The expression is: In the formula, The rise time of the impact pressure is expressed by the following formula: The attenuation constants α and β are determined by the following empirical relationship: In the formula, These are characteristic parameters of pressure decay; The longitudinal wave velocity of the rock is expressed in m / s. The borehole diameter is in meters (m).

6. A method for simulating a blasting process based on bidirectional fluid-structure interaction according to claim 3, 4, or 5, characterized in that, Step S3, solving for the unsteady seepage process of the explosive gas, specifically includes: S3.1: Map the fracture spatial coordinates obtained in step S2 to the fluid computation domain grid cells to identify the affected grid regions in the fluid computation domain; S3.2: Locally increase porosity in each grid cell of the fluid computational domain containing the fracture. And update the penetration rate according to the power law relationship: In the formula, Initial penetration rate, Where is the initial porosity, and n is an empirical exponent; Porosity is defined as: In the formula, This represents the number of cracks in a single grid cell within the fluid computational domain. The porosity increment for each fracture; S3.3: The magnitude of the residual stress at the borehole wall at the moment the shock wave action ends in step S2 is taken as the initial gas pressure of the fluid computation domain. This is used as the initial pressure condition at the inlet boundary of the fluid computational domain; S3.4: Solve the unsteady seepage process of explosive gas in fractured media. The governing equation is composed of the mass conservation equation coupled with Darcy's law in compressible form; the gas mass conservation equation is expressed as: In the formula, ρ is the gas density, in kg / m³; Porosity; v is Darcy velocity, volumetric flow rate per unit area, in m / s; gas flow follows a modified Darcy's law: In the formula, Permeability, in m²; μ, gas viscosity, in MPa·s; P, gas pressure, initially P0, in MPa; S3.5: Let the time step of the fluid computational domain be... This allows us to obtain the spatial distribution of the gas pressure field within the fluid computational domain at the current time step.

7. The blasting process simulation method based on bidirectional fluid-structure interaction according to claim 6, characterized in that, In step S4, the initiation and propagation of new cracks specifically include: S4.1: Using an interpolation method, the pressure value of the current time step on the grid cell of the fluid computation domain is smoothly mapped to the center position of each particle in the solid computation domain to obtain the equivalent gas pressure on the particle; S4.2: According to time step The solid computational domain is calculated using a DEM solver, new cracks are generated, and their spatial information is recorded.

8. A blasting process simulation system based on two-way fluid-structure interaction, used to implement the blasting process simulation method based on two-way fluid-structure interaction as described in any one of claims 1 to 7, characterized in that, include: The parameter input module is used to input all the initial conditions and material properties required for the simulation. The particle mechanics solution module is used to simulate the dynamic response and spontaneous crack evolution of rocks under blasting loads in the solid computation domain, and to extract the fracture information of rocks during the shock wave stage. The gas diffusion solution module is used to map fracture information to the fluid computational domain, solve the unsteady seepage process of explosive gas in the fracture, and obtain the spatial distribution of the gas pressure field in the fluid computational domain. The fluid-structure interaction module uses scripts to exchange real-time data of the crack information obtained by the particle mechanics solution module and the gas pressure field obtained by the gas diffusion solution module, thereby establishing a fluid-structure interaction mechanism and realizing bidirectional coupling between the fluid computation domain and the solid computation domain. The results output module is used to output the fracture network morphology, gas pressure evolution cloud map, crack propagation path, and response time history of key monitoring points.

9. A blasting process simulation system based on bidirectional fluid-structure interaction according to claim 8, characterized in that, The parameters input by the parameter input module include the elastic modulus, tensile strength, cohesion, and internal friction angle of the solid computational domain, and the initial gas pressure, gas viscosity, gas density, initial porosity, initial permeability, and power-law exponent of permeability evolution with porosity of the fluid computational domain. It also includes coupling control parameters such as time step, total simulation duration, peak pressure of Duvall pressure time-history load, normalization factor, and decay constant.

Images