A blasting process simulation method and system based on two-way fluid-structure coupling

Abstract

The application 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 firstly adopts a discrete element method to model a solid calculation domain and adopts a finite volume method to model a fluid calculation domain; subsequently, a dynamic load is applied to simulate initial cracks in a shock wave stage, crack information is mapped to a fluid calculation domain grid, local permeability is dynamically updated, a compressible gas diffusion equation is solved, and a pressure field is applied to particle surfaces of the solid calculation domain through an interpolation method to drive crack continuous expansion. The application establishes a fluid-solid bidirectional feedback mechanism covering blast gas and rock, accurately captures a driving effect of the gas at a crack tip through a closed loop coupling of cracks-permeability-gas pressure-new cracks, can not only accurately reproduce experimental crack shapes, but also quantitatively analyze influences of initial gas pressure on a broken range, crack quantity and a damage mode, and is suitable for blasting optimization design under a high gas yield explosive condition.

Description

Technical Field

[0001] This invention relates to the field of rock blasting mechanics and numerical simulation technology, specifically to a blasting process simulation method and system based on two-way fluid-structure interaction, which is suitable for simulating the wedging, loading and feedback behavior of high-pressure explosive gas in a dynamically evolving fracture network during blasting. Background Technology

[0002] In rock blasting, while the initial stress wave can induce dense microcracks, the explosive gas in the subsequent millisecond-level quasi-static stage is also a key factor in driving crack penetration, continuous energy injection, and eventual rock fragmentation. However, the explosive gas possesses complex characteristics such as instantaneous high pressure, strong nonlinear permeability, and dynamic feedback with the fracture network. Its real-time diffusion path, local pressure distribution, and driving mechanism on crack tips within the rock mass are difficult to observe directly through experiments, necessitating high-fidelity numerical methods to reveal these characteristics.

[0003] Currently, numerical simulation of blasting mainly employs the finite element method (FEM), extended finite element method (XFEM), or discrete element method (DEM). While FEM / XFEM can describe the deformation of continuous media well, they struggle to spontaneously simulate the initiation, branching, and penetration of complex cracks. DEM can naturally reproduce the entire rock fracture process, but existing loading methods often employ attenuated pressure (equivalent pressure method) applied to the borehole wall or the central particle expansion method, neglecting the actual migration process of gas along the fracture network. This results in no pressure loading inside the fracture, failing to reflect the gas accumulation effect at the crack tip and the dynamic correlation between gas diffusion and permeability evolution.

[0004] While existing methods have improved the precision of blast simulation to some extent, they still have significant limitations. For example, Chinese patent CN114218831B proposes simulating crack propagation by pre-setting joint / interface elements between solid elements. Although this avoids mesh re-division and can handle complex fracture morphologies, its fracture path still relies on prior settings and does not introduce any physical models of gas diffusion or flow, thus failing to reflect the dynamic driving effect of explosive gases on crack propagation. Another example is Chinese patent CN115859714A, which constructs an FEM-DEM co-simulation framework, applying the explosion pressure time history calculated by FEM as a boundary load to the DEM model. While this can take into account both shock wave and quasi-static pressure effects, the gas effect is strictly limited to the initial borehole wall or pulverization zone, failing to track the real-time migration process of gas to the newly formed cracks on the periphery, and also failing to consider the local pressure accumulation effect of gas at the crack tip. For example, Chinese patent CN109631701B focuses on the macroscopic load superposition and vibration prediction of multi-hole blasting. It improves the accuracy of vibration velocity simulation by applying blasting stress to the tunnel contour surface, but it does not involve the diffusion behavior of blasting gas in the fracture network inside the rock mass and lacks a physical description of the millisecond-level quasi-static fracture evolution mechanism.

[0005] It is evident that existing methods generally fail to incorporate the infiltration and flow mechanism of explosive gases within the fracture network. They typically simplify gas action to a pre-defined pressure-time history on the borehole wall, completely ignoring the physical nature of discrete fractures as the dominant channels for gas diffusion. In reality, crack propagation opens new high-conductivity channels, while gas pressure, in turn, drives further crack propagation, forming a strongly nonlinear feedback loop. Ignoring this two-way coupling mechanism leads to severe distortion of the explosive fragmentation evolution process.

[0006] Therefore, there is an urgent need for a numerical simulation device and method that can synchronously, bidirectionally, and dynamically couple gas diffusion and rock fracture to realistically reproduce the entire process of migration, loading, and feedback of explosive gases in evolving fractures, providing scientific support for understanding blasting mechanisms and optimizing engineering parameters. Summary of the Invention

[0007] The purpose of this invention is to provide a method for simulating the blasting process based on two-way fluid-structure interaction. This invention integrates the finite volume method (FVM) and the discrete element method (DEM), utilizes a custom interface to achieve two-way data interaction, and constructs a closed-loop feedback mechanism of "crack identification—permeability update—gas solution—pressure mapping—particle loading—new crack generation," realistically simulating the diffusion and driving behavior of blasting gases in a dynamic fracture network.

[0008] This invention employs a phased loading strategy: first, a dynamic load based on Duvall pressure time history is applied to simulate the initial crack during the shock wave stage; then, a fluid-structure interaction module is activated to map crack information onto a fluid mesh, dynamically update local permeability, solve the compressible gas diffusion equation, and apply a pressure field to the particle surface via interpolation, driving continuous crack propagation. This method not only accurately reproduces experimental crack morphology but also quantitatively analyzes the impact of initial gas pressure on the fragmentation range, crack number, and failure mode, making it particularly suitable for blasting optimization design under high gas yield explosive conditions.

[0009] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0010] A method for simulating blasting processes based on two-way fluid-structure interaction includes:

[0011] S1: The discrete element method (DEM) is used to model the solid computational domain, and the rock medium is represented by densely packed spherical particles; the finite volume method (FVM) is used to model the fluid computational domain, and the diffusion behavior of gas in it is described using a grid region with the same geometry as the solid computational domain.

[0012] 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.

[0013] 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.

[0014] 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.

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

[0016] Furthermore, in step S1, the process of modeling the solid computational domain and the fluid computational 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 a parallel bond model (PBM), and key mechanical parameters such as elastic modulus, tensile strength, cohesion, and friction angle are defined independently.

[0018] S1.2: In the Finite Volume Method (FVM) solver, a mesh with the exact same geometry as the solid computational domain is created as the fluid computational domain;

[0019] 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;

[0020] S1.4: Initialize the porosity and permeability of the solid computational domain as the baseline diffusion parameters in the unruptured state.

[0021] Furthermore, in step S2, the specific process for extracting the initial fracture information of the rock during the shock wave stage is as follows:

[0022] 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:

[0023]

[0024] 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;

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

[0026] 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.

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

[0028]

[0029] In the formula, The rise time of the impact pressure is expressed by the following formula:

[0030]

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

[0032]

[0033] In the formula, These are characteristic parameters of pressure attenuation. The longitudinal wave velocity of the rock is expressed in m / s. The borehole diameter is in meters (m).

[0034] Furthermore, in step S3, solving the unsteady seepage process of the explosive gas specifically includes:

[0035] 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;

[0036] 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:

[0037]

[0038] In the formula, Initial penetration rate, Where is the initial porosity, and n is an empirical exponent; Porosity is defined as:

[0039]

[0040] 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;

[0041] 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;

[0042] 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:

[0043]

[0044] Where ρ 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:

[0045]

[0046] In the formula, Permeability is expressed in m²; μ is the gas viscosity, expressed in MPa·s. Let $\frac{ ... Unit: MPa;

[0047] 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.

[0048] Furthermore, step S4, which drives the initiation and propagation of new cracks, specifically includes:

[0049] 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;

[0050] 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.

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

[0052] A blasting process simulation system based on two-way fluid-structure interaction includes:

[0053] The parameter input module is used to input all the initial conditions and material properties required for the simulation, including 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 the 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.

[0054] 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.

[0055] 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.

[0056] 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.

[0057] 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.

[0058] The beneficial effects of this invention are:

[0059] (1) This invention realizes the simulation of the interaction between explosive gas and rock in dynamic fracture network, breaks through the conventional method of traditional pore wall loading, and significantly improves physical realism;

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

[0061] (3) This invention can separate the interaction between shock wave and gas, support mechanism research, and reveal the control law of shock wave and gas on the fragmentation mode through parametric analysis.

[0062] In summary, this invention provides a high-fidelity, scalable, and physically clear simulation method for explosive gas-fracture coupling, which solves the core problems in existing technologies such as gas loading distortion, missing coupling, and inability to reflect the real diffusion path. It not only deepens the understanding of the explosive breaking mechanism, but also provides a new research tool for the precise control and safety optimization of blasting in mining, tunnel and other engineering projects. Attached Figure Description

[0063] Figure 1 This is a flowchart of the simulation method of the present invention.

[0064] Figure 2 Modeling of fluid and solid computational domains.

[0065] Figure 3 Calibrate the dynamic mechanical parameters of the split Hopkinson bar (SHPB).

[0066] Figure 4 The stress curve is shown for the shock wave stage.

[0067] Figure 5 The experiment and simulation methods are compared under the condition of no gas action; (a) shows the stress distribution of the simulated blasting process of the present invention, and (b) shows the strain distribution of the blasting test.

[0068] Figure 6 The simulation shows the entire process of blasting under the action of explosive gases; where (a) is the stress distribution, (b) is the permeability distribution, and (c) is the gas stress distribution.

[0069] Figure 7 The following is a comparison of stress conditions in the simulation of the embodiments of the present invention; wherein (a) is the stress curve under no gas action, and (b) is the stress curve under the action of explosive gas.

[0070] Figure 8 The comparison shows the destruction modes; (a) is the simulation result and (b) is the experimental result. Detailed Implementation

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

[0072] like Figure 1 As shown, this invention provides a method for simulating the blasting process based on two-way fluid-structure interaction, comprising the following steps:

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

[0074] S1.1: In the Discrete Element Method (DEM) platform, spherical particles are used to model the rock sample. The spherical particles are connected by a parallel bond model (PBM), and key mechanical parameters such as elastic modulus, tensile strength, cohesion, and friction angle are defined independently.

[0075] S1.2: In the Finite Volume Method (FVM) solver, a mesh with the exact same geometry as the solid computational domain is created as the fluid computational domain;

[0076] 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, while the outer boundary condition of the fluid computational domain is set as a free pressure outlet boundary, allowing gas to flow out freely, thereby simulating gas diffusion behavior in an open environment;

[0077] S1.4: Initialize the porosity of the solid computational domain With penetration rate , which serves as the baseline diffusion parameter in the unbroken state.

[0078] S2: Based on the solid computational 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 fracture information of rock during the shock wave stage is extracted simultaneously, including fracture initiation, propagation path and spatial distribution; the specific process is as follows:

[0079] 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:

[0080]

[0081] in, For at any time Dynamic stress applied to the borehole wall in the solid computational domain; α and β are frequency-dependent attenuation constants; Peak pressure, determined by the type of explosive, can be expressed as:

[0082]

[0083] Where, ρ e ρ is the density of the explosive, in g / cm³; D is the detonation velocity, in m / s;

[0084] It is a normalization factor used to ensure that the pressure curve is reasonably normalized in the rising segment. Its expression is:

[0085]

[0086] in, The rise time of the impact pressure can be obtained using the following formula:

[0087]

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

[0089]

[0090] in, These are characteristic parameters of pressure attenuation. The longitudinal wave velocity of the rock is expressed in m / s. The borehole diameter is in meters (m).

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

[0092] 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.

[0093] S3: Based on the fluid computational domain constructed in S1, and combined with the crack information generated by the rupture of the solid computational domain in S2, the unsteady seepage process of the explosive gas is solved, and the spatial distribution of the gas pressure field in the fluid computational domain at the current time step is obtained; as follows:

[0094] 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;

[0095] 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:

[0096]

[0097] In the formula, Initial penetration rate, Where is the initial porosity, and n is an empirical exponent; Porosity is defined as:

[0098]

[0099] 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;

[0100] 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;

[0101] 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:

[0102]

[0103] 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:

[0104]

[0105] In the formula, k is the permeability, in m²; μ is the gas viscosity, in MPa·s; P is the gas pressure, initially the gas pressure is... Unit: MPa;

[0106] 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.

[0107] S4: Map the gas pressure field of the current time step in the fluid computational domain to the solid computational domain, applying it as an additional force to the particles in the solid computational domain to drive the initiation and propagation of new cracks; specifically:

[0108] S4.1: Using an interpolation method, the pressure value P 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;

[0109] S4.2: According to time step The solid domain is calculated using a DEM solver, new cracks are generated, and their spatial information is recorded.

[0110] S5: Extract new fracture information and return to step S3 to solve for gas flow in the next time step; repeat steps S3 to S4 to form a closed-loop feedback of "fracture-permeability-gas pressure distribution-new fracture" until the gas pressure decays to a level that is insufficient to drive new fractures, and the calculation ends.

[0111] Example

[0112] In the blast test of polycarbonate (PC) plates conducted by Yang et al. (Visualizing the blast-induced stress wave and blasting gasaction effects using digital image correlation), the specimen geometry was 300 mm × 300 mm × 5 mm, the central bore diameter was 4 mm, and lead azide was used for blasting. Explosives. The dynamic elastic modulus of 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 invention, corresponding numerical simulations were carried out for the above-mentioned experimental cases. To systematically evaluate the influence of shock wave and explosive gas on rock mass fracture, two sets of comparative working conditions were set up: (1) shock wave action - only the transient stress load of the shock wave stage is applied, the subsequent gas pressure is not considered, and only the solid computational domain is activated for solution; (2) combined action - both shock wave stress and the seepage pressure of explosive gas in the fracture are considered, and both the solid computational domain and the fluid domain are activated for solution. By comparing the results of the two sets, the key contribution of explosive gas in the crack propagation and penetration process can be clearly revealed.

[0114] The Duvall pressure time-history load was used during the shock wave phase, and the relevant parameters were as follows: peak pressure P m = 293MPa, normalization factor ξ = 3.47, attenuation constant α = 0.0362 μs -1 β = 0.0808 μs -1 The time step for the numerical simulation is set to 1×10. -9 The time history of the shock wave was measured in seconds, with a total duration of 100 μs, to accurately capture the initiation process of initial microcracks under the action of the shock wave. The shock wave stress-time history curve is shown below. Figure 4 As shown. After entering the explosive gas action stage, the initial gas pressure is taken. = 27 MPa (corresponding to the residual pressure at the end of shock wave decay, see...) Figure 4 () is used as the inlet boundary condition for the fluid domain; the gas viscosity is taken as μ = 3.0 × 10⁻⁶. -5Pa·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 increase caused by the crack 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 5In Figure (b), the strain distribution within the elastic vibration range measured by Yang et al. in their experiment is shown. After conversion using Hooke's Law, the stress distribution is approximately in the range of 0–12 MPa. It can be seen that the experimental results and the numerical simulation in this embodiment show good consistency in terms of stress distribution shape and amplitude.

[0118] Figure 6 The simulation results of the explosion under the combined action of stress wave and explosive gas are presented. It can be seen that the compressive stress distribution is not significantly different when gas is involved compared to when no gas is present; the main difference lies in the tensile stress distribution. In the initial stage of gas action, the maximum tensile stress inside the sample reaches approximately 100 MPa, at which point the permeability in the crack region increases significantly, and the gas pressure reaches 27 MPa. Under continuous loading of the explosive gas, the crack propagates outward from the borehole, with tensile stress concentrating at the tip, and the permeability further increases in the crack region. Simultaneously, due to energy dissipation and crack propagation, the gas pressure gradually decreases. When the gas pressure is insufficient to continue driving crack propagation, crack propagation ceases.

[0119] Figure 7 The stress-time response along the axial direction inside the specimen is shown under different working conditions. Figure 7 Figure (a) shows the stress evolution curve considering only the stress wave action. It can be seen that the maximum stress peak is 51.37 MPa, which decays rapidly within 1.5 ms, indicating that this stage is mainly controlled by shock wave propagation, has a short duration, and releases energy rapidly. In contrast, Figure 7 (b) shows the stress response process after the coupled explosive gas action. The peak axial stress rises to 83.59 MPa, which is much higher than the case without gas action. Moreover, the entire stress decay process is significantly prolonged, and there is still a high residual stress after 6 ms.

[0120] Furthermore, as the measuring point moved further away from the explosion source, the peak stress showed a significant attenuation trend, but the residual stress under the action of gas remained at a high level, reflecting that the explosive gas continuously provided loading during the quasi-static stage, effectively extending the time window for crack propagation. This delayed loading behavior facilitated the propagation, penetration, and connection of existing cracks, thereby promoting larger-scale structural failure. Therefore, the gas action not only increased the tensile stress level at the crack tip but also enhanced the driving force for crack propagation.

[0121] Furthermore, the rock mass failure modes obtained from experiments by Yang et al. and numerical simulations in this embodiment were compared, such as... Figure 8 As shown, the two exhibit high consistency in crack initiation, propagation direction, and failure mode. The results demonstrate that the numerical simulation method proposed in this invention can accurately simulate the synergistic mechanism of stress wave propagation and gas-induced fracturing, verifying its applicability in simulating the entire blasting process.

[0122] Finally, it should be noted that the above embodiments are intended to illustrate the technical solutions of the present invention and do not constitute any limitation on the present invention. Those skilled in the art should fully understand that modifications to the technical solutions described in the foregoing embodiments or equivalent substitutions for any part or all of the technical features are entirely feasible. Such modifications or substitutions, as long as they do not depart from the scope of protection defined by the claims of the present invention, should be considered reasonable extensions of the present invention.

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 computation domain and the solid computation 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. 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 under the unruptured state; 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, The initial porosity, n It is an experience index; 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 gas mass conservation equation coupled with Darcy's law in compressible form; the gas mass conservation equation is expressed as: In the formula, ρ This refers to the density of a gas, expressed in kg / m³. 3 ; Porosity; v Darcy velocity, volumetric velocity per unit area, in m / s; gas flow obeys a modified Darcy's law: In the formula, Permeability, in meters (m) 2 ; μ The viscosity is expressed in MPa·s. P Let $\frac{ ... P 0 Unit: 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.

2. The method for simulating the blasting process based on bidirectional fluid-structure interaction according to claim 1, 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 β It is a frequency-dependent attenuation constant; 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.

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

4. The blasting process simulation method based on bidirectional fluid-structure interaction according to claim 2, 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 constant α and β This is determined by the following empirical relationship: In the formula, These are characteristic parameters of pressure attenuation. The longitudinal wave velocity of the rock is expressed in m / s. The borehole diameter is in meters (m).

5. A method for simulating a blasting process based on bidirectional fluid-structure interaction according to claim 2, 3, or 4, 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.

6. A blasting process simulation system based on bidirectional fluid-structure interaction, used to implement the blasting process simulation method based on bidirectional fluid-structure interaction as described in any one of claims 1 to 5, 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.

7. The blasting process simulation system based on bidirectional fluid-structure interaction according to claim 6, 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