A whole-process simulation method and system for tunnel seepage and water and mud inrush

Abstract

The application provides a whole process simulation method and system for tunnel seepage and water and mud inrush, and the method comprises the following steps: constructing a numerical model containing a tunnel and a fault, building a solid particle and a fluid particle system by using a meshless particle; constructing a random fracture network and mapping to the solid particle, and assigning attributes according to the fracture filling condition; calculating the interaction between particles, introducing a seepage-stress coupling model in the fracture area, and dynamically updating the permeability characteristics to simulate the implicit seepage; judging the damage and rupture of the solid particle by using a preset criterion, and converting into a discrete element particle to represent the cracked block when the condition is met; calculating the contact force between the discrete element particles and the undamaged particles, and simulating the movement of the broken rock mass; building a two-way fluid-solid coupling, and realizing the whole process simulation of seepage incubation, rock mass rupture and water and mud inrush through time integration. The application seamlessly simulates the strong coupling process of rock mass rupture and large deformation fluid migration by using the DEM-SPH coupling.

Description

Technical Field

[0001] This invention relates to the field of tunnel engineering technology, and in particular to a method and system for simulating the entire process of tunnel seepage, water inrush, and mud inrush. Background Technology

[0002] In the construction of deep and complex geological tunnels, sudden water and mud inrush disasters are the primary safety threat due to their suddenness and high destructiveness. The essence of these disasters lies in the implicit seepage within a random network of fractures ahead of the tunnel face, which, under specific conditions, triggers explicit dynamic failure of the rock mass, resulting in transient dynamic disasters. Current technologies face three interconnected fundamental challenges when simulating this entire process:

[0003] First, the underlying mechanisms of the disaster are unclear: heterogeneous seepage caused by random fractures is difficult to characterize precisely. The occurrence, geometry, and connectivity of rock fractures are highly random, making seepage paths concealed and unpredictable. Traditional continuous media methods cannot effectively capture the subtle influence of this randomness on the seepage field, making it difficult to reveal the key physical mechanisms in the hidden seepage process. This results in an inability to accurately identify potential sources of water inrush risk and insufficient early warning accuracy.

[0004] Secondly, the catastrophic process is difficult to reproduce: the transient dynamic process of large deformation and strong coupling is difficult to simulate effectively. Once a catastrophic event occurs, it evolves into a strong fluid-structure interaction transient problem involving rock mass fracture, large displacement, and high-speed fluid and sediment transport. Traditional mesh-based methods (such as FEM) cannot handle large deformation collapse due to mesh distortion; while meshless fluid methods (such as finite element method FEM and FDM) cannot natively simulate the fracture and solid contact behavior of rock mass. Both single methods have inherent limitations.

[0005] Finally, there is a problem with methodological coordination: cross-scale multiphysics coupling covering the entire process presents technical barriers. Ideal simulations require the integration of multiple models, such as fracture seepage, rock mass fracturing, and fluid transport. However, how to efficiently and stably couple the Discrete Element Method (DEM), which excels at simulating rock mass fracturing, with the Smoothed Particle Hydrodynamics (SPH), which excels at simulating large deformation fluid flow, to achieve seamless transfer and bidirectional coupling of physical fields, remains a long-standing technical challenge.

[0006] Therefore, there is an urgent need for a full-process simulation method and system for tunnel seepage and water and mud inrush, which can achieve seamless simulation of the strongly coupled process of rock mass fracture and large deformation fluid transport through DEM-SPH coupling. Summary of the Invention

[0007] The purpose of this invention is to provide a method and system for simulating the entire process of tunnel seepage and sudden water and mud inrush, aiming to solve the technical problem that existing technologies are unable to accurately simulate the entire process of tunnel seepage from random fissures to sudden water and mud inrush.

[0008] To achieve the above objectives, in a first aspect, the present invention provides a method for simulating the entire process of tunnel seepage and sudden water and mud inrush, the steps of which include:

[0009] S1. Construct a numerical model including the tunnel structure and fault geological conditions, and establish a solid particle system to characterize the rock mass and a fluid particle system to characterize the groundwater using a meshless particle method; preferably, the solid particles and fluid particles are initialized separately.

[0010] S2. Construct a random fracture network in the numerical model, map the fractures to the solid particle system, and assign corresponding mechanical or seepage properties to the particles in the fracture region according to whether the fractures are filled with a weak medium.

[0011] S3. Calculate the interactions between fluid particles and solid particles, and introduce a seepage-stress coupling model in the fracture region. Dynamically update the fracture permeability characteristics according to the stress state of the rock mass to simulate the implicit seepage process in the fault rock mass. Preferably, the SPH-based method is used to calculate the interactions between fluid particles and solid particles.

[0012] S4. The solid particles are damaged and fractured by using a preset yield criterion. When the solid particles meet the failure conditions, the corresponding particles are transformed from continuous medium particles into discrete element particles (DEM particles) to characterize the cracking and block formation of the rock mass.

[0013] S5. Calculate the contact forces between the discrete element particles and between the discrete element particles and the undamaged solid particles to simulate the sliding, rolling and transport behavior of the rock blocks after fracturing.

[0014] S6. Establish a two-way fluid-structure interaction relationship between fluid particles, solid particles, and discrete element particles. Through time integral iterative calculation, realize the simulation of the entire process of seepage incubation, rock mass fracturing, and water and mud inrush in fault tunnels.

[0015] As a further improvement to the above scheme, in step S2, the step of establishing the random fracture network includes:

[0016] S11. Obtain the fracture characteristics of the rock mass through image recognition or probabilistic statistical methods to generate a Discrete Fracture Network (DFN) that conforms to the actual geological characteristics.

[0017] S12. Map the crack features generated in step S11 to the numerical model, mark the complete SPH particles covering the cracks, and convert them into DEM particles or remove them directly according to the crack properties.

[0018] As a further improvement to the above scheme, in step S3, the fracture permeability characteristics are dynamically updated using the cubic law, and their value is adjusted according to the change of normal stress on the fracture surface, specifically calculated by the following formula:

[0019] ;

[0020] in, Where is the permeability coefficient; g is the acceleration due to gravity. The dynamic viscosity of the fluid. and These are the maximum and minimum principal stresses, respectively. Poisson's ratio, The normal deformation coefficient is... The stiffness coefficient of the crack; For water head.

[0021] As a further improvement to the above scheme, in step S1, when initializing solid particles and fluid particles, the solid particles are given initial position, velocity, stress, material parameters (such as elastic modulus, internal friction angle, cohesion, tensile strength) and boundary conditions.

[0022] Assign initial position, velocity, and water parameters to fluid particles;

[0023] The virtual particle boundary method is used to apply a repulsive force to prevent particles from penetrating the boundary. Type I virtual particles are preferred to apply a boundary repulsive force to the internal particles to prevent them from penetrating the boundary.

[0024] As a further improvement to the above scheme, the repulsive force exerted by the boundary on the internal particles As shown below:

[0025] ;

[0026] in, The cutoff radius is determined by the initial spacing between the particles. This represents the position difference between two paired particles. It is the distance between the inner particle and the boundary virtual particle. ≥ At that time, the repulsive force no longer applies; and These are all empirical index parameters, typically set to 12 and 4 respectively, used to adjust the nonlinear relationship between repulsive force and distance; This represents the position vector between the internal particles and the boundary virtual particles. It is the unit direction vector of that position vector; It is the intensity coefficient of the boundary repulsion force, an empirical parameter determined by the specific problem (such as the type of fluid being simulated and the particle parameters), and its magnitude is usually equivalent to the square of the maximum particle velocity.

[0027] As a further improvement to the above scheme, in step S3, when calculating the interaction force between particles:

[0028] S31. Perform intelligent pairing of neighboring particles: Based on the linked list search algorithm and distance threshold optimization, first filter neighboring particles within the influence domain of the kernel function, and at the same time add particle type label to avoid invalid pairing;

[0029] S32. Perform force calculations by type:

[0030] Fluid-fluid interaction calculation: Based on the discretized Navier-Stokes equations, considering viscous fluids with shear stress, calculate the interaction forces between fluid particles;

[0031] Solid-solid interaction calculations: When using the SPH momentum equation for solid particles, an elastoplastic constitutive correction is added (a yield criterion is added to the stress tensor) to replace the simple viscous term, which better reflects the deformation characteristics of solids.

[0032] S33. Solve for the pressure term using the EoS equation of state:

[0033] Equations of state In this process, the reference sound velocity is dynamically adjusted with the particle density ρ (to avoid pressure distortion in low-density regions). Different adiabatic indices γ are used for fluid / solid particles; γ is a dimensionless parameter, preferably 7. ρ0 is the reference density, preferably 1000 kg / m³. 3 c0 is the calculated numerical speed of sound, which is a constant, preferably 10;

[0034] S34. Post-processing verification of force:

[0035] After the calculation is completed, the momentum conservation of the particles is checked (the forces between adjacent particles are equal in magnitude and opposite in direction). If the condition is not met, the discrete scheme of the partial derivatives of the kernel function is finely adjusted.

[0036] As a further improvement to the above scheme, in step S31, the momentum equation on which the fluid-fluid interaction calculation is based is as follows:

[0037] ;

[0038] Where a and b represent the fundamental fluid particle and its neighboring fluid particles within its influence domain, respectively. The total number of neighboring solid particles within the influence domain of the kernel function of the basic particle a; Let be the rate of change of the density of particle a with time, and α and β be the coordinate components. Artificial viscosity, Let be the relative velocity component of particle a in the α direction. Let be the partial derivative of the kernel function with respect to the α-coordinate direction of particle a. Let be the relative velocity components of particles a and b in the α direction. Let be the acceleration component of particle a in the α direction. Let be the component of the external force acting on particle a in the α direction. For kernel function, Let b be the mass of the neighboring particle. The mass of the basic fluid particle a; Let be the stress tensor of particle a. P represents particle pressure. The Kronecker delta function (1 when α=β, 0 otherwise) The dynamic viscosity of the fluid. For the strain rate tensor components.

[0039] As a further improvement to the above solution, the artificial viscosity is as follows:

[0040] ;

[0041] in, and Here, c represents the control parameter for artificial viscosity, and c represents the numerical velocity of sound. The position difference between two paired particles. The velocity difference between the two paired particles This represents the average density of the particles.

[0042] As a further improvement to the above scheme, in step S4, the Drucker-Prager damage criterion with tensile truncation is used to determine the damage and fracture of solid particles, in order to distinguish between tensile failure and shear failure of rock mass, and to determine the crack initiation mode and propagation direction.

[0043] As a further improvement to the above scheme, the preset yield criterion is as follows:

[0044] ;

[0045] In the formula For tensile strength, For the maximum principal stress, and These are the first and second invariants of the stress tensor, respectively. and All are Drucker-Prager constants;

[0046] If any one of the conditions is met, the particle is determined to have undergone tensile or shear failure.

[0047] As a further improvement to the above scheme, in step S5, the contact forces between the discrete element particles and between the discrete element particles and the undamaged solid particles are calculated based on the penalty function contact model. The penalty function contact model considers both normal contact force and shear contact force, and constrains the slippage behavior between particles according to Coulomb's friction law.

[0048] As a further improvement to the above scheme, the construction steps of the penalty function contact model are as follows: first, update the position and influence domain of each DEM particle using a linked list search algorithm, then perform particle pairing, and finally, set a contact threshold. From the formula It is confirmed that, among them, , Let d be the radius of the two DEM particles, and let d be half the initial particle spacing ΔP / 2. id The distance between two DEM particles;

[0049] When U n When the value is greater than 0, it is determined that DEM particles have come into contact, and the contact force is... Including normal contact force and shear contact force and satisfy ;

[0050] Normal contact force From the formula Calculated, where Let the normal stiffness be the influence domain of the mass. This is the unit outward normal vector of the particle at the current moment;

[0051] Shear contact force From the formula Calculated, where Where is the shear stiffness in the formula. Let ΔU be the tangential contact force acting on particle i at the current time step. s Let be the tangential component of the relative displacement increment between particles I and j.

[0052] As a further improvement to the above scheme, in step S3, the introduction of a seepage-stress coupling model to dynamically update the fracture permeability characteristics specifically includes:

[0053] Seepage velocity determined using the SPH discretization method With seepage head The evolution equation is used to update the seepage field:

[0054] ;

[0055] ;

[0056] in, Where is the water yield of the rock, and N is the total number of neighboring particles involved in the calculation. Let be the mass of the j-th neighboring particle. Let be the density of the j-th neighboring particle. Let be the seepage velocity at the i-th particle. Let be the seepage velocity at the j-th particle; and These are the seepage heads of the i-th and j-th particles, respectively; This refers to the smoothing kernel function in the SPH method. For smooth kernel function For the coordinates of the i-th particle The partial derivatives are used to calculate the interaction strength between particles;

[0057] For situations involving water level boundaries or rainfall infiltration, a specific inflow boundary condition term is superimposed on the seepage velocity calculation. The details are as follows:

[0058] ;

[0059] Flow-stress coupling is achieved by superimposing the flow force onto the Cauchy stress of the solid particles, as shown below:

[0060] ;

[0061] , and These are the corresponding deviatoric stress components; For hydrostatic pressure, This is the normal deformation coefficient.

[0062] As a further improvement to the above scheme, in step S6, a two-way fluid-structure interaction relationship is established between fluid particles, solid particles, and discrete element particles. The force exerted by the fluid on the rock mass and the reaction force exerted by the rock mass on the fluid are calculated according to Newton's third law, so as to realize explicit fluid-structure interaction simulation in the process of water and mud inrush.

[0063] In a second aspect, the present invention also provides a tunnel water inrush and mud inrush full-process simulation system for implementing the method provided in the first aspect, comprising:

[0064] The model building module is used to construct numerical models including tunnel and fault structures, and to initialize solid and fluid particles.

[0065] A fracture network generation module is used to generate and map a random fracture network to the numerical model;

[0066] The seepage-stress coupling calculation module is used to calculate the coupled evolution of rock mass stress and seepage field in fractured regions;

[0067] The damage and fracture determination module is used to determine the fracture of rock mass particles based on a preset yield criterion and to realize the transformation of continuous particles into discrete element particles.

[0068] The contact and block motion simulation module is used to calculate the contact forces between particles after fracturing and to simulate the movement of rock blocks;

[0069] The fluid-structure interaction module is used to realize bidirectional interaction between fluid particles, solid particles, and discrete element particles;

[0070] The time integration and result output module is used to iteratively calculate and output the evolution results of the entire process of water and mud inrush.

[0071] As a further improvement to the above technical solution, the seepage-stress coupling calculation module is configured to update the fracture permeability parameters in real time when the rock mass stress changes.

[0072] As a further improvement to the above technical solution, the damage and fracture determination module and the contact and block motion simulation module achieve automatic connection between the continuous domain and the discontinuous domain through a particle type conversion mechanism.

[0073] Because the present invention adopts the above technical solutions, the beneficial effects of this application are as follows:

[0074] This invention provides a full-process simulation method for tunnel seepage and water / mud inrush. First, a random fracture network is constructed in the numerical model and the fractures are mapped to a solid particle system. At the same time, particles in the fracture region are assigned corresponding mechanical or seepage properties according to whether they are filled with weak media. Combined with the introduction of a seepage-stress coupling model in the fracture region and the dynamic updating of fracture permeability characteristics based on the stress state of the rock mass, this method can accurately reflect the influence of the randomness of fracture orientation, geometry and connectivity on the seepage field under a meshless particle framework. This allows for a more accurate revelation of the key physical mechanisms of hidden seepage and helps improve the accuracy of identifying and warning of water inrush risk sources.

[0075] Secondly, the SPH method is used to calculate the interaction between fluid particles and solid particles, which can avoid the distortion problem of traditional mesh methods under large deformations without a mesh framework. Furthermore, when solid particles meet the preset yield criterion, they are transformed into discrete element particle (DEM) particles, and the contact force between blocks after fracture is calculated from the DEM. This can natively characterize rock mass fracture, block formation and its sliding, rolling and transport behavior, thereby realizing continuous and stable simulation of the strongly coupled transient process of rock mass fracture and high-speed fluid sediment transport.

[0076] Furthermore, by establishing a two-way fluid-structure interaction relationship between fluid particles, solid particles, and discrete element particles, and performing time-integral iterative calculations, the collaborative operation and seamless transfer of physical fields between the SPH fluid module and the DEM solid fracture module were achieved. This coupling method overcomes the inherent limitations of single methods in simulating different physical processes, enabling multiple physical fields such as fracture seepage, rock mass fracturing, and fluid migration to evolve synchronously within the same numerical framework, thereby covering the entire process of water and mud inrush disaster incubation, occurrence, and evolution.

[0077] The combination of various technical features enables the present invention to complete the continuous simulation of random fissure seepage, rock mass rupture, and water and mud inrush in the same system. It has both mechanism revelation and process reproduction functions, and can provide a quantitative analysis tool for risk assessment and safety control of deep and complex geological tunnel construction, thereby improving the pertinence and reliability of disaster prediction and prevention. Attached Figure Description

[0078] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the structures shown in these drawings without creative effort.

[0079] Figure 1 This is a schematic diagram of the random discrete fracture network generation process in the simulation method disclosed in this invention;

[0080] Figure 2 This is a comparative diagram showing the traditional and improved particle neighborhood search methods disclosed in this invention in addressing the problem of "avoiding erroneous nuclear interactions between particles from different broken rock blocks";

[0081] Figure 3 This is a schematic diagram of the boundary configuration and basic particle distribution of the Zhongjiashan Tunnel model disclosed in this invention.

[0082] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

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

[0084] It should be noted that all directional indicators (such as up, down, etc.) in the embodiments of the present invention are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicator will also change accordingly.

[0085] Furthermore, in this invention, descriptions involving "first," "second," etc., are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined with "first" or "second" may explicitly or implicitly include at least one of that feature.

[0086] Furthermore, the technical solutions of the various embodiments of the present invention can be combined with each other, but only if they are based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0087] Example 1

[0088] See Figures 1-3 This invention provides a method for simulating the entire process of tunnel seepage and sudden water and mud inrush, the implementation process of which includes the following steps:

[0089] S1. Numerical Model Construction and Particle Initialization:

[0090] Based on geological survey data from actual tunnel engineering projects, a numerical model including tunnel structure and fault geological conditions was established. A meshless particle discretization method was used to generate solid particle systems representing rock mass and fluid particle systems representing groundwater. During the initialization of solid and fluid particles, initial positions, velocities, stress states, and material parameters were assigned to solid particles, including elastic modulus, internal friction angle, cohesion, and tensile strength. Initial positions, velocities, and water density and viscosity parameters were assigned to fluid particles. A virtual particle boundary method was employed to prevent particles from penetrating the boundary; Type I virtual particles used a repulsive force formula to constrain the boundaries of internal particles. Particle initialization provides physically consistent initial conditions for subsequent seepage-stress coupling and dynamic failure simulations.

[0091] Specifically, the repulsive force formula As shown below:

[0092] ;

[0093] in, The cutoff radius is determined by the initial spacing between the particles. This represents the position difference between two paired particles. It is the distance between the inner particle and the boundary virtual particle. ≥ At that time, the repulsive force no longer applies; and These are all empirical index parameters, typically set to 12 and 4 respectively, used to adjust the nonlinear relationship between repulsive force and distance; This represents the position vector between the internal particles and the boundary virtual particles. It is the unit direction vector of that position vector; It is the intensity coefficient of the boundary repulsion force, an empirical parameter determined by the specific problem (such as the type of fluid being simulated and the particle parameters), and its magnitude is usually equivalent to the square of the maximum particle velocity.

[0094] S2. Construction and attribute mapping of random fracture networks:

[0095] In numerical models, stochastic discrete fracture networks (DFNs) are generated based on field statistics or image recognition techniques. See also Figure 1 The crack search algorithm maps the crack to a solid particle system, specifically including:

[0096] S11. Obtain the fracture characteristics of the rock mass through image recognition or probabilistic statistical methods to generate a Discrete Fracture Network (DFN) that conforms to the actual geological characteristics.

[0097] S12. Map the crack features generated in step S11 to the numerical model, mark the complete SPH particles covering the cracks, and convert them into DEM particles or remove them directly according to the crack properties.

[0098] S13. Depending on whether the crack is filled with weak material, assign differentiated mechanical or seepage properties to the particles in the crack area. For example, assign lower strength when filled with weak material, and directly convert to DEM particles or remove them when not filled.

[0099] This step improves the accuracy of implicit seepage field simulation by refining the characterization of the heterogeneous distribution of fractures, providing a foundation for the evolution of seepage paths during the disaster incubation stage.

[0100] S3. Interaction Calculation and Seepage-Stress Coupling:

[0101] The SPH method is used to calculate the interaction forces between fluid particles and solid particles. The fluid momentum equation is solved based on the discretized Navier-Stokes equations, and solid deformation is simulated using the SPH form of the momentum equation. A seepage-stress coupling model is introduced in the fractured region. The fracture permeability coefficient is dynamically updated based on the rock mass stress state, and the seepage field is corrected in real time using the SPH discretized seepage velocity and head update equations. This coupling mechanism can reflect the influence of stress changes on fracture conductivity, achieving a high-precision characterization of water pressure distribution and flow evolution during implicit seepage.

[0102] S4. Rock Mass Damage and Fracture Determination and Particle Transformation

[0103] Damage assessment of solid particles is performed using the Drucker-Prager yield criterion with tensile truncation. When the particle stress state meets the failure condition, it is transformed from a continuous medium SPH particle into a discrete element particle (DEM particle) to characterize the cracking and block formation of the rock mass. This transformation mechanism achieves a seamless transition from continuous deformation to discontinuous fracture in the rock mass, overcoming the inherent limitation of the traditional SPH method in simulating fracture.

[0104] S5. Contact Force Calculation and Block Motion Simulation:

[0105] Contact forces are calculated using the penalty function method for both the transformed DEM particles and the undamaged solid particles. Collision and friction behavior between blocks are described using normal and tangential stiffness models, and Coulomb's friction law is introduced to determine the sliding state. The motion of the DEM particles is solved using Newton's second law, thereby simulating the sliding, rolling, and migration behavior of the fractured rock blocks, achieving a stable reproduction of the explicit dynamic failure process of the rock mass.

[0106] S6. Two-way fluid-structure interaction and full-process iterative simulation:

[0107] Establish a two-way coupling relationship between fluid particles and solid particles:

[0108] The force exerted by the fluid on the solid is calculated using the pressure contribution term, and the reaction force exerted by the solid on the fluid is calculated using the formula. accomplish;

[0109] To avoid confusion in density calculations, the particle density contributions in the coupling region are differentiated.

[0110] For the interaction between DEM particles and fluid, a virtual particle repulsion force model is used to ensure numerical stability.

[0111] By iteratively updating particle position, velocity, and stress state using a time integration algorithm, a dynamic simulation of the entire process from fissure seepage incubation and rock mass fracturing to water and mud inrush disasters is ultimately achieved.

[0112] In a preferred embodiment, in step S3, the fracture permeability characteristics are dynamically updated using the cubic law, and their value is adjusted according to the change in normal stress on the fracture surface. Specifically, the fracture permeability coefficient is calculated using the following formula:

[0113] ;

[0114] in, Permeability coefficient; The dynamic viscosity of the fluid. and These are the maximum and minimum principal stresses, respectively. Poisson's ratio, The normal deformation coefficient is... The stiffness coefficient of the crack; For water head; specifically, maximum principal stress. As shown in the following formula:

[0115] ;

[0116] Minimum principal stress As shown in the following formula:

[0117] ;

[0118] The fracture permeability coefficient is dynamically adjusted with normal stress, which can accurately depict the fracture opening and closing effect caused by changes in rock mass stress state, and significantly improve the simulation accuracy of implicit seepage field.

[0119] The stress-flow coupling is implemented as follows:

[0120] When calculating the rock mass stress field at each time step, the principal stress values ​​of particles in each fracture region are obtained simultaneously.

[0121] Based on the cubic law formula, the fracture permeability coefficient is updated in real time according to the current stress state.

[0122] Substitute the updated permeability coefficient into the SPH discrete form of the seepage velocity equation and the head update equation;

[0123] When considering the effect of gravity, the total head difference is corrected to Δh. ij =h i -h j -y i -y j ;

[0124] By using the seepage-stress coupling mechanism, key parameters such as water pressure distribution and flow evolution in fracture networks can be accurately simulated, providing a physical basis for the early identification of water inrush risk sources. The seepage equation based on the SPH discrete format is updated synchronously with the stress field, avoiding the problem of discontinuous physical field transmission in traditional methods and ensuring the numerical stability of the entire simulation process.

[0125] As a preferred embodiment, when calculating the interparticle interaction force in step S3, the calculation accuracy and numerical stability are improved through the following refined implementation method, the specific steps of which include:

[0126] S31, Intelligent Pairing of Neighboring Particles:

[0127] Based on a linked list search algorithm, each fundamental particle is paired with neighboring particles within its kernel function's influence domain. The search range is optimized by setting a distance threshold, and particle type labels, such as fluid particles, solid particles, and DEM particles, are added to avoid invalid pairings between different particle types. This intelligent pairing mechanism reduces computational redundancy and improves the efficiency of particle interaction force calculation.

[0128] S32. Calculation of different types of forces:

[0129] Fluid-fluid interaction: Based on the discretized Navier-Stokes equations, considering a viscous fluid with shear stress, the interaction forces between fluid particles are calculated. The momentum equation is shown below:

[0130] ;

[0131] Where a and b represent the fundamental fluid particle and its neighboring fluid particles within its influence domain, respectively. The total number of neighboring solid particles within the influence domain of the kernel function of the basic particle a; Let be the rate of change of the density of particle a with time, and α and β be the coordinate components. Artificial viscosity, Let be the relative velocity component of particle a in the α direction. Let be the partial derivative of the kernel function with respect to the α-coordinate direction of particle a. Let be the relative velocity components of particles a and b in the α direction. Let be the acceleration component of particle a in the α direction. Let be the component of the external force acting on particle a in the α direction. For kernel function, Let b be the mass of the neighboring particle. The mass of the basic fluid particle a; Let be the stress tensor of particle a; P is the particle pressure. The Kronecker delta function (1 when α=β, 0 otherwise) The dynamic viscosity of the fluid. For the strain rate tensor components.

[0132] An artificial viscosity formula is introduced to maintain numerical stability. Specifically, the artificial viscosity formula is as follows:

[0133] ;

[0134] in, and Here, c represents the control parameter for artificial viscosity, and c represents the numerical velocity of sound. The position difference between two paired particles. The velocity difference between the two paired particles This represents the average density of the particles.

[0135] Solid-solid interactions: The internal forces and deformations of solid particles are simulated using momentum equations in the form of SPH, and elastoplastic constitutive corrections are added to the stress tensor. By introducing the Drucker-Prager yield criterion with tensile truncation to replace the simple viscous term, the constitutive relation is made to better fit the elastoplastic deformation characteristics of rock masses, thus improving the physical realism of solid deformation simulation.

[0136] Specifically, the momentum equation of the SPH form is shown in the following equation:

[0137] ;

[0138] The calculation of force by type takes into account both the viscosity of fluids and the elastic-plastic properties of solids, thus enhancing the realism of physical processes.

[0139] S33. State Equation Optimization and Parameter Setting:

[0140] The pressure term is solved using the equation of state (EoS), and the formula is as follows:

[0141] ;

[0142] Where ρ0 is the reference density, which is equal to 1000 kg / m³. 3 c0 is the calculated numerical sound velocity, a constant set to 10, and γ is a dimensionless parameter set to 7. To further optimize the calculation, the reference sound velocity is dynamically adjusted with the particle density ρ to avoid pressure distortion in low-density regions. Simultaneously, a different adiabatic index γ is used for fluid and solid particles to more accurately reflect the compressibility characteristics of different media. Parameter optimization of the equation of state avoids numerical distortion in low-density regions and improves the stability of the full-domain simulation.

[0143] S34, Post-processing verification of force

[0144] After calculating the interaction forces, the momentum conservation of the particles is verified, i.e., whether the forces between adjacent particles satisfy the relationship of equal magnitude and opposite direction. If the momentum conservation is not satisfied, the discrete scheme of the partial derivatives of the kernel function is fine-tuned, and the physical rationality of the mechanical response is ensured through iterative correction. The momentum conservation verification ensures the physical rationality of the mechanical response, providing a reliable input for subsequent fluid-structure interaction.

[0145] In a preferred embodiment, in step S4, the Drucker-Prager damage criterion with tensile truncation is used to determine the damage and fracture of solid particles, in order to distinguish between tensile failure and shear failure of rock mass, and to determine the crack initiation mode and propagation direction.

[0146] Specifically, the preset yield criterion is achieved through the following formula:

[0147] ;

[0148] In the formula For tensile strength, The maximum principal stress is calculated using the stress tensor components, specifically by the following formula:

[0149] ;

[0150] and These are the first and second invariants of the stress tensor, respectively. Specifically, the formula for calculating the first invariant is as follows:

[0151] ;

[0152] in , as well as This refers to, for example The case where the superscripts α and β are both x, y, and z.

[0153] The formula for calculating the second invariant is as follows:

[0154] ;

[0155] and All are Drucker-Prager constants, obtained through calculations using rock mass material parameters.

[0156] , Where c represents cohesion. It is the internal friction angle;

[0157] Specifically, the damage assessment process is as follows:

[0158] After calculating the stress state of the solid particles at each time step, their stress state is calculated simultaneously. , and value;

[0159] Determine whether the particle satisfies any of the destruction conditions:

[0160] like It was determined to be a shear failure;

[0161] like It was determined to be a tensile failure;

[0162] If any condition is met, the particle is determined to be damaged, and it is transformed from a continuous medium SPH particle into a discrete element DEM particle.

[0163] By jointly judging shear failure conditions (Drucker-Prager criterion) and tensile failure conditions (maximum principal stress criterion), the failure mechanism of rock mass under different stress states can be accurately identified, providing a theoretical basis for the analysis of fracture initiation modes. The introduction of the tensile truncation mechanism makes the criterion more consistent with the actual failure characteristics of rock materials, avoiding unreasonable judgments that may occur under high hydrostatic pressure due to the simple Drucker-Prager criterion. Once the particle meets the failure condition, it is transformed into a DEM particle, realizing the natural transition of rock mass from continuous deformation to discontinuous fracture, laying the foundation for subsequent block movement simulation.

[0164] In a preferred embodiment, in step S5, the present invention calculates the contact forces between discrete element particles and between discrete element particles and undamaged solid particles based on a penalty function contact model. This model considers both normal contact force and shear contact force, and constrains the slippage behavior between particles using Coulomb's friction law.

[0165] The steps for constructing the penalty function contact model are as follows:

[0166] First, a linked list search algorithm is used to perform neighbor pairing for each particle. When the distance between two particles meets the contact condition, a contact calculation is triggered.

[0167] Traditional linked list methods are passively used to search for particles within the fluid domain. To mitigate the SPH nucleus interaction between opposite IPs in the DEM, this invention also provides an improved particle search method, such as... Figure 2 As shown. Among them. Figure 2 (a) shows the particle neighborhood search under the traditional kernel radius. The red dashed circle in the figure indicates the range covered by the traditional kernel radius h, centered on a target particle. It can be seen that when a larger kernel radius, such as h / ΔP>1 (ΔP being the initial particle spacing), is used, the search range crosses the boundaries of different "broken rocks," including particles belonging to other rocks in the neighborhood. Green spheres represent discrete particles (DPs), and gray spheres represent background fluid or other particles. Under the traditional method, a large kernel radius easily causes particles from different rocks to be incorrectly included in the same neighborhood, leading to unnecessary SPH kernel function interactions and affecting computational accuracy or efficiency.

[0168] Figure 2(b) shows the particle neighborhood search under the improved kernel radius. In this figure, the kernel radius is limited (h / ΔP≤1), and the new neighborhood range is represented by a red dashed circle. It can be seen intuitively that the reduced kernel radius will no longer "cross the boundary" into other broken rock blocks, but only covers the physically adjacent particles near the current rock block. The core purpose of this is to avoid erroneous SPH kernel interactions between particles in different broken rock blocks, so that when simulating "fluid-particle (or rock block-particle) interaction", it can both ensure the mechanical relationship within a single block and reduce cross-block interference, thereby improving the stability and computational efficiency of the algorithm. Figure 2 By comparing the coverage of particle neighborhoods under two kernel radius strategies, this invention clearly demonstrates how the improved approach of "reducing the kernel radius to h / ΔP≤1" solves the problem of "false interactions between particles across rock blocks." This invention employs a particle neighborhood search with an improved kernel radius, where the ratio of the smoothed kernel radius to the initial particle spacing is h / ΔP, and this ratio is limited to h / ΔP≤1. In simulating fluid interactions, we use a traditional kernel radius h / ΔP>1, which includes particles within a range greater than 2h.

[0169] Contact threshold From the formula It is confirmed that, among them, , Let d be the radius of the two DEM particles, and let d be half the initial particle spacing ΔP / 2. id The distance between two DEM particles;

[0170] When U n When the value is greater than 0, it is determined that DEM particles have come into contact, and the contact force is... Including normal contact force and shear contact force and satisfy ;

[0171] Normal contact force From the formula Calculated, where Let the normal stiffness be the influence domain of the mass. This is the unit outward normal vector of the particle at the current moment;

[0172] Shear contact force From the formula Calculated, where Where is the shear stiffness in the formula. Let ΔU be the tangential contact force acting on particle i at the current time step. s Let be the tangential component of the relative displacement increment between particles I and j.

[0173] Coulomb friction constraint:

[0174] The tangential contact force is constrained according to Coulomb's law of friction to prevent unreasonable sliding behavior, as shown in the following formula:

[0175] ;

[0176] Where u is the coefficient of friction. When the tangential force reaches the friction limit, the contact force is adjusted according to the sliding state.

[0177] Integrating the contact force into the momentum equation, the governing equation for the complete particle is revised as follows:

[0178] ;

[0179] The motion of DEM particles can be solved using Newton's second law, and its final governing equation is:

[0180] ;

[0181] Where g is gravity (9.81 m / s²), and the internal forces are... From the formula The calculation yields α and β, which represent the x, y, and z directions in the Cartesian coordinate system, respectively, i.e., the Einstein summation convention, applied to repeated indexes. Represented as the Cauchy stress tensor, Let P be the Dirac function, and P be the isotropic pressure. It is the deviatoric stress tensor.

[0182] ;

[0183] in It is the density of the initial particles. It is particle density. Young's modulus can be calculated from the elastic modulus E and Poisson's ratio. For anisotropic shear stress, assuming small displacement, the stress rate and strain rate are proportional. The proportionality coefficient is given by the shear modulus G = E / (2 × (1 + ... E is the elastic modulus. Given Poisson's ratio, the equation for the stress rate tensor is:

[0184] ;

[0185] in , The strain rate tensor is defined as:

[0186] ;

[0187] To align the material information with strain, a continuity equation is derived by introducing the Jaumann rate of change:

[0188] ;

[0189] in, and All are stress tensor components. and It is the torsion rate tensor, defined as:

[0190] .

[0191] The penalty function model can accurately describe the complex contact behaviors such as collision and sliding between rock blocks after fracturing, improving the realism of block motion simulation. The separate calculation of normal and tangential contact forces, combined with Coulomb friction constraints, fully reproduces the mechanical interaction mechanism between rock blocks. By setting the contact threshold reasonably and introducing the friction law, the penetration and oscillation phenomena in numerical calculation are effectively avoided. The efficient combination of the linked list search algorithm and contact detection ensures the feasibility of large-scale particle contact calculation.

[0192] In a preferred embodiment, in step S3, the introduction of the seepage-stress coupling model to dynamically update the fracture permeability characteristics is achieved through the following specific methods:

[0193] Discrete solution of SPH for seepage field update:

[0194] The evolution equations of seepage velocity and seepage head are solved using the SPH discretization method, enabling real-time updates of the seepage field.

[0195] Seepage velocity calculation:

[0196] ;

[0197] in, Let be the mass of the j-th neighboring particle. Let be the density of the j-th neighboring particle. The seepage velocity component at the i-th particle; The current fracture permeability coefficient, , This represents the head value of the corresponding particle.

[0198] For situations involving water level boundaries or rainfall infiltration, a specific inflow boundary condition term is superimposed on the seepage velocity calculation. The details are as follows:

[0199] .

[0200] Water head evolution equation:

[0201] ;

[0202] in, Let N be the water yield of the rock, and N be the total number of neighboring particles involved in the calculation within the kernel function's influence domain. Let be the time rate of change of the seepage head at the i-th particle position, describing the evolution of the head over time.

[0203] Boundary condition handling:

[0204] For cases with water level boundaries or rainfall infiltration, specific inflow boundary conditions are superimposed in the seepage velocity calculation. When considering gravity, the total head difference is corrected as follows:

[0205] ;

[0206] and These represent the position coordinates of the i-th and j-th particles in the direction of gravity, respectively.

[0207] The seepage velocity equation is adjusted accordingly to:

[0208] ;

[0209] Seepage-stress coupling is achieved as follows:

[0210] Flow-stress coupling is achieved by superimposing the flow force onto the Cauchy stress of the solid particles, as shown below:

[0211] ;

[0212] , and These are the corresponding deviatoric stress components; For hydrostatic pressure, The normal deformation coefficient is... For water density, This represents the current head value.

[0213] The SPH discrete scheme effectively describes the non-uniform seepage process in fractured rock masses and realistically reflects the water-rock interaction by combining a stress coupling mechanism. Through head difference correction and boundary term superposition, it can accurately simulate actual engineering conditions such as water level boundaries and rainfall infiltration. The seepage force is directly embedded in the stress tensor, achieving bidirectional real-time coupling between the seepage field and the stress field, avoiding numerical errors caused by physical field transmission in traditional methods. The particle method-based discrete scheme naturally adapts to large deformation problems, ensuring the numerical stability of seepage-stress coupling calculations during water and mud inrush processes. This invention, through a refined seepage-stress coupling mechanism, provides reliable technical support for simulating the incubation process of tunnel water and mud inrush disasters, significantly improving the accuracy of implicit seepage stage simulation.

[0214] As a preferred embodiment, a two-way fluid-structure interaction relationship is established between fluid particles, solid particles, and discrete element particles, thereby achieving explicit fluid-structure interaction simulation of the water and sludge inrush process in the following manner:

[0215] When fluid particles enter the influence domain of solid particles, the solid particles, acting as a deformable boundary, exert a force on the fluid particles. This force is mainly contributed by the fluid pressure. According to Newton's third law, the fluid particles simultaneously exert equal and opposite forces on the solid particles. To avoid confusion in density calculations, the density contributions of fluid and solid particles are distinguished in the coupling region. Ultimately, the conservation equations for the fluid and solid are integrated to incorporate their respective forces. The coupling force vector from the fluid particle acting on the i-th solid particle is shown in the following equation:

[0216] ;

[0217] in It is the pressure of the i-th solid particle. It is the pressure of the a-th fluid particle, given by the formula To obtain. Let be the mass of the i-th solid particle. Let ρ be the mass of the a-th fluid particle in the neighborhood. i Let ρ be the density of the i-th solid particle. a Let be the density of the a-th fluid particle; Let x be the gradient of the smooth kernel function W with respect to the coordinates of the i-th particle. ia =x i -x a It is the relative position vector between particles i and a;

[0218] According to Newton's third law, the equal reaction force of adjacent solid particles on fluid particles, the coupling force vector from the solid particles on the a-th fluid particle is shown in the following equation:

[0219] ;

[0220] in, Let W be the gradient of the smooth kernel function with respect to the coordinates of the a-th particle.

[0221] Furthermore, considering that the sum of the densities of fluid and rock particles are separated when fluid and solid particles are close to each other, this modification prevents the density of rock particles from accumulating in the sum of the densities of the fluid domain, and vice versa. These improvements enhance the numerical stability between the fluid and solid domains.

[0222] The momentum equation for solid particles is rewritten as follows:

[0223] ;

[0224] The momentum equation for fluid particles is rewritten as follows:

[0225] ;

[0226] The first term represents the equation for a complete particle, and the second term represents the equation for a fluid particle.

[0227] DEM rock mass particles are not controlled by the SPH kernel function equation; their trajectories are determined by Newton's second law. The interaction forces between the fluid and discontinuous rock mass particles differ from those of the SPH kernel function equation. According to Newton's third law and based on the penalty function theory, when particles approach each other, a repulsive force is generated due to distance. This repulsive force is similar to that of type I virtual particles. In this case, discontinuous rock particles interacting with fluid particles are considered virtual particles, and their interaction forces... and These are the reaction forces between the fluid particles and the DEM particles, respectively. These forces are considered as external forces and are expressed as follows:

[0228] ;

[0229] The momentum equation for fluid particles is rewritten as follows:

[0230] ;

[0231] The momentum equation for DEM rock mass particles is rewritten as:

[0232] ;

[0233] By using time integration algorithms, such as the frog-leapfrog algorithm, to iteratively update particle position, velocity, and stress state, a dynamic simulation of the entire process from fracture seepage incubation and rock mass fracturing to the development of water and mud inrush disasters is ultimately achieved. Specifically, the iterative update formula is as follows:

[0234] ;

[0235] In the formula, t and t0 represent the calculation time and the initial time, respectively; Δt is the calculation time step; ρ is the particle density; v is the particle velocity; and x is the particle's position coordinates. The change rate of water head for the particles.

[0236] To further illustrate the inventive concept of this invention, a tunnel project is used as an example. A numerical model is established based on actual geological and tunnel construction conditions, and the simulation method provided by this invention simulates the entire process from fissure seepage incubation to the development of water and mud inrush disasters.

[0237] Model building and parameter setting:

[0238] See Figure 3 Based on the actual conditions of the Zhongjiashan Tunnel project, a numerical model with a height of 35 m and a length of 35 m was established. The tunnel face height is 9 m, the thickness of the water-bearing fault is 15 m, and the dip angle is 81°. The density of water is set as ρ = 1000 kg / m³. 3 The initial spacing between rock and water particles, ΔP, is 0.125 m, and the length of the SPH smooth core, h, is 0.125 m. The SPH parameters of the tunnel model are shown in Table 1.

[0239] ;

[0240] Model boundary conditions are set as follows: Gravitational acceleration g = 9.81 m / s² is applied as the driving force for water inrush and rock movement; a vertical ground stress load equivalent to a 200 m overburden layer is applied to the upper boundary of the tunnel, the lower boundary is fixed, and the lateral boundary is fixed in the horizontal direction; for water bodies, a water column height of 150 m is considered, and the corresponding water pressure is applied.

[0241] The specific implementation steps are as follows:

[0242] 1. Model Construction and Particle Initialization:

[0243] The tunnel fault structure was generated, with the entire model consisting of 60,250 rock particles and 18,150 water particles within the fault. A random fracture network was also established at the tunnel face. Appropriate material parameters were assigned to the rock particles, and in the initial stage, infinite shear and tensile strengths of the rock mass were used for stress balance to prevent premature rock damage.

[0244] 2. Initial Equilibrium and Seepage Simulation:

[0245] When the expected in-situ stress reaches initial equilibrium, the actual tensile strength and shear strength are calculated, with the calculation time set to 0 seconds. A seepage boundary is applied at the contact point between the rock mass and the water body, and a non-permeable boundary is set around it. Seepage simulation is first performed with a time step of 1 min / step for about 2 days. The seepage situation in the fractures is observed, and the seepage velocity and fracture development are recorded.

[0246] 3. Preparation for dynamic calculations:

[0247] After the dominant fractures are formed, explicit hydrodynamic and fluid-structure interaction calculations are performed. The dynamic viscosity is set to 1×10⁻⁶. - 6 Pa·s, with a time step of 5 × 10⁻⁶. -5 s / step, totaling 120,000 steps, covering the entire process of water and mud inrush.

[0248] 4. Single-time-step iterative calculation:

[0249] Particle pairing and interaction calculations: Based on searching for neighboring particles with twice the smooth core radius, the pressure of water particles and the stress of rock particles are calculated.

[0250] Seepage-stress coupling update: Update the seepage head and stress state of fractured rock mass according to the coupled cubic law;

[0251] Damage and fracture determination: The Drucker-Prager criterion with tensile truncation is used to determine the fracture state of particles. Particles that meet the failure conditions are converted into DEM particles.

[0252] Fluid-structure interaction calculations: Calculate the interaction forces between fluids and intact particles and DEM particles;

[0253] Time integration: The key parameters are integrated over time using the frog-leap algorithm to obtain particle density, water pressure, velocity and displacement.

[0254] 5. Loop control and result output:

[0255] Determine whether the calculation step has reached the set maximum value. If not, repeat the iterative calculation in step 4. If the maximum calculation step has been reached, terminate the calculation and output the results such as particle position, velocity, stress and crack distribution to achieve a quantitative evaluation of the entire process of tunnel seepage failure and water and mud inrush.

[0256] This method can overcome the limitations of traditional numerical methods in dealing with large deformation and strong coupling problems, and provides a reliable technical means for the study of the mechanism and prevention of water and mud inrush disasters in tunnel engineering.

[0257] Example 2

[0258] The present invention also provides a simulation system for the entire process of tunnel water inrush and mud inrush for implementing the method described in Embodiment 1, comprising:

[0259] This invention provides a simulation system for the entire process of water and mud inrush in tunnels to implement the method described in Embodiment 1. It achieves a complete simulation process from model construction to result output through modular design. The specific composition of each module is shown below:

[0260] The model building module is used to construct a numerical model including tunnel and fault structures, and to initialize solid and fluid particles. In practice, a geometric model is established based on actual engineering geological conditions. A meshless particle discretization method is used to generate a solid particle system representing the rock mass and a fluid particle system representing groundwater. Initial positions, velocities, stress states, and material parameters such as elastic modulus, internal friction angle, and cohesion are assigned to the solid particles, while initial positions, velocities, and water parameters are assigned to the fluid particles. Through the initialization function of this module, initial states consistent with actual engineering conditions are provided for subsequent physical process simulations.

[0261] The fracture network generation module is used to generate random fracture networks and map them to numerical models. During implementation, a discrete fracture network (DFN) is generated based on probabilistic statistical distributions or image recognition techniques. A fracture search algorithm maps the fractures to a solid particle system, and the particles in the fracture region are assigned corresponding mechanical or seepage properties based on whether the fractures are filled with weak materials. This module lays the foundation for accurate simulation of implicit seepage fields by finely characterizing the heterogeneous distribution of fractures.

[0262] The seepage-stress coupling calculation module is used to calculate the coupled evolution of rock mass stress and seepage field in fractured regions. As an improvement, this module is configured to update fracture permeability parameters in real time as rock mass stress changes. Specifically, it dynamically adjusts the fracture permeability coefficient using the cubic law, achieving bidirectional coupling between the stress field and the seepage field. This real-time update mechanism accurately reflects the impact of stress changes on fracture conductivity, improving the physical realism of the seepage field simulation.

[0263] The damage and fracture determination module determines rock mass particle fracture based on a preset yield criterion (such as the Drucker-Prager criterion with tensile truncation) and converts continuous particles into discrete element particles. As an improvement, this module automatically connects the continuous and discontinuous domains with the contact and block motion simulation module through a particle type conversion mechanism. When the particle stress state meets the failure condition, it is automatically converted from SPH continuous particles to DEM discrete particles, ensuring a natural transition of the rock mass from continuous deformation to discontinuous fracture.

[0264] The contact and block motion simulation module is used to calculate the contact forces between particles after fracturing and to simulate the movement of rock blocks. Based on the penalty function contact model, it calculates the contact forces between DEM particles and between DEM particles and undisturbed solid particles, while considering both normal and shear contact forces. It also constrains the slippage behavior between particles using Coulomb's law of friction, accurately simulating the slippage, tumbling, and migration processes of rock blocks after fracturing.

[0265] The fluid-structure interaction module realizes bidirectional interactions between fluid particles, solid particles, and discrete element particles. Based on Newton's third law, it calculates the force exerted by the fluid on the rock mass and the reaction force exerted by the rock mass on the fluid. Density contribution differentiation ensures the numerical stability of the coupling region, enabling explicit fluid-structure interaction simulation of water and mud inrush processes.

[0266] The time integration and result output module uses the frog-leap algorithm for time integration iterative calculation, and outputs key parameters such as particle position, velocity, stress and fracture distribution in real time. It fully reproduces the entire process evolution from seepage incubation and rock mass fracture to water and mud inrush disasters, providing quantitative basis for disaster mechanism analysis and engineering prevention and control.

[0267] This system provides a complete simulation chain from model building to result output, covering the entire process of water and mud inrush disasters. It also ensures the physical realism of the simulation through mechanisms such as real-time coupling of seepage and stress and automatic conversion between continuous and discontinuous processes. Modular design and intelligent particle management enhance the feasibility of large-scale computing. Furthermore, the system provides rich and intuitive output parameters to support disaster risk assessment and prevention and control design.

[0268] The above are merely preferred embodiments of the present invention and do not limit the patent scope of the present invention. All equivalent structural transformations made using the contents of the present invention's specification and drawings under the inventive concept of the present invention, or direct / indirect applications in other related technical fields, are included within the patent protection scope of the present invention.

Claims

1. A method for simulating the entire process of tunnel seepage and sudden water / mud inrush, characterized in that, include: S1. Construct a numerical model including tunnel structure and fault geological conditions, and establish a solid particle system to characterize rock mass and a fluid particle system to characterize groundwater using a meshless particle method. S2. Construct a random fracture network in the numerical model, map the fractures to the solid particle system, and assign corresponding mechanical or seepage properties to the particles in the fracture region according to whether the fractures are filled with a weak medium. The steps for establishing the random fracture network include: S11. Obtain the fracture characteristics of the rock mass through image recognition or probabilistic statistical methods to generate a discrete fracture network that conforms to the actual geological characteristics; S12. Map the crack features generated in step S11 to the numerical model, label the complete SPH particles covering the cracks, and convert them into DEM particles or remove them directly according to the crack properties. S3. First, intelligent pairing of neighboring particles is performed through a linked list search algorithm and distance threshold optimization. Then, the interaction forces between fluid particles and solid particles are calculated according to their types. At the same time, a seepage-stress coupling model is introduced in the fracture region. Based on the cubic law, the fracture permeability characteristics are dynamically updated according to the stress state of the rock mass to simulate the implicit seepage process in the fault rock mass. S4. A preset yield criterion is used to determine the damage and fracture of solid particles. When the solid particles meet the failure conditions, the corresponding particles are transformed from continuous medium particles into discrete element particles to characterize the cracking and block formation of the rock mass. S5. Calculate the contact forces between discrete element particles and between discrete element particles and undamaged solid particles based on the penalty function contact model. The penalty function contact model considers both normal contact force and shear contact force, and constrains the slippage behavior between particles according to Coulomb's friction law to simulate the slippage, rolling and migration behavior of rock blocks after fracturing. S6. Establish a two-way fluid-structure interaction relationship between fluid particles, solid particles, and discrete element particles. Through time integral iterative calculation, realize the simulation of the entire process of seepage incubation, rock mass fracturing, and water and mud inrush in fault tunnels.

2. The method for simulating the entire process of tunnel seepage and sudden water / mud inrush as described in claim 1, characterized in that, In step S3, the fracture permeability characteristics are dynamically updated using the cubic law, and their values ​​are adjusted according to the change in normal stress on the fracture surface, specifically calculated by the following formula: ； in, Where is the permeability coefficient; g is the acceleration due to gravity. The dynamic viscosity of the fluid. and These are the maximum and minimum principal stresses, respectively. Poisson's ratio, The normal deformation coefficient is... The stiffness coefficient of the crack; For water head.

3. The method for simulating the entire process of tunnel seepage and sudden water / mud inrush as described in claim 1, characterized in that, In step S3, the method for calculating the interaction forces between particles is as follows: S31. Perform intelligent pairing of neighboring particles: Based on the linked list search algorithm and distance threshold optimization, first filter neighboring particles within the influence domain of the kernel function, and at the same time add particle type labels to avoid invalid pairings. S32. Perform force calculations by type: Fluid-fluid interaction calculation: Based on the discretized Navier-Stokes equations, considering viscous fluids with shear stress, calculate the interaction forces between fluid particles; Solid-solid interaction calculations: When applying the SPH momentum equation to solid particles, an elastoplastic constitutive correction is added to replace the simple viscous term, so as to better fit the deformation characteristics of solids; S33. Using state equations The pressure term P is solved, and the reference sound velocity is dynamically adjusted with the particle density ρ. Different adiabatic indices γ are used for fluid / solid particles. ρ0 is the reference density and c0 is the calculated numerical sound velocity. S34. After the calculation is completed, check the momentum conservation of the particle under force. If it is not satisfied, fine-tune the discrete format of the partial derivative of the kernel function.

4. The method for simulating the entire process of tunnel seepage and sudden water / mud inrush as described in claim 1, characterized in that, In step S4, the Drucker-Prager damage criterion with tensile truncation is used to determine the damage and fracture of solid particles, which is used to distinguish between tensile failure and shear failure of rock mass, so as to determine the crack initiation mode and propagation direction.

5. The method for simulating the entire process of tunnel seepage and sudden water / mud inrush as described in claim 4, characterized in that, The preset yield criterion is shown in the following formula: ； In the formula For tensile strength, For the maximum principal stress, and These are the first and second invariants of the stress tensor, respectively. and All are Drucker-Prager constants; If any one of the conditions is met, the particle is determined to have undergone tensile or shear failure.

6. The method for simulating the entire process of tunnel seepage and sudden water / mud inrush as described in claim 1, characterized in that, The steps for constructing the penalty function contact model are as follows: First, the position and influence domain of each DEM particle are updated using a linked list search algorithm, then particle pairing is performed, and a contact threshold is set. From the formula It is confirmed that, among them, , Let d be the radius of the two DEM particles, and let d be half the initial particle spacing ΔP / 2. id The distance between two DEM particles; When U n When the value is greater than 0, it is determined that DEM particles have come into contact, and the contact force is... Including normal contact force and shear contact force and satisfy ; Normal contact force From the formula Calculated, where Let the normal stiffness be the influence domain of the mass. This is the unit outward normal vector of the particle at the current moment; Shear contact force From the formula Calculated, where Where is the shear stiffness in the formula. Let ΔU be the tangential contact force acting on particle i at the current time step. s Let be the tangential component of the relative displacement increment between particles I and j.

7. The method for simulating the entire process of tunnel seepage and sudden water / mud inrush as described in claim 1, characterized in that, In step S6, a two-way fluid-structure interaction relationship is established between fluid particles, solid particles, and discrete element particles. The force exerted by the fluid on the rock mass and the reaction force exerted by the rock mass on the fluid are calculated according to Newton's third law, so as to realize explicit fluid-structure interaction simulation in the process of water and mud inrush.

8. A simulation system for the entire process of tunnel water and mud inrush as described in any one of claims 1-6, characterized in that, include: The model building module is used to construct numerical models including tunnel and fault structures, and to initialize solid and fluid particles. A fracture network generation module is used to generate and map a random fracture network to the numerical model; The seepage-stress coupling calculation module is used to calculate the coupled evolution of rock mass stress and seepage field in fractured regions; The damage and fracture determination module is used to determine the fracture of rock mass particles based on a preset yield criterion and to realize the transformation of continuous particles into discrete element particles. The contact and block motion simulation module is used to calculate the contact forces between particles after fracturing and to simulate the movement of rock blocks; The fluid-structure interaction module is used to realize bidirectional interaction between fluid particles, solid particles, and discrete element particles; The time integration and result output module is used to iteratively calculate and output the evolution results of the entire process of water and mud inrush.

Images