Sph-dem coupled system and method for preventing blowout of a screw-type excavator

Abstract

This invention discloses a blowout prevention system and method for auger excavators based on SPH-DEM coupling, belonging to the field of tunnel engineering construction technology. The system integrates an SPH fluid simulation module, a DEM particle simulation module, a bidirectional coupling interface, a parametric model of the auger excavator, and a blowout early warning module. It uses the Lagrange particle method to discretize the fluid, Hertz-Mindlin theory to solve for particle contact forces, and GPU parallel computing to achieve bidirectional exchange of forces between the fluid and particles. It also integrates multiple key physical quantities to construct a blowout risk criterion. The method includes multi-scale model construction, coupling calculation parameter setting, bidirectional coupling calculation, blowout feature extraction, risk early warning, and parameter optimization steps. This invention simulates the particle migration and blowout mechanism caused by water pressure disturbance in water-rich strata, improving the accuracy of blowout risk identification and early warning capabilities, and providing a scientific basis for the structural design and construction parameter control of auger excavators.

Description

Technical Field

[0001] This invention relates to the field of tunnel engineering construction technology, and in particular to a blowout prevention system and method for a spiral excavator based on SPH-DEM coupling. Background Technology

[0002] In complex shield tunneling projects such as subway construction and cross-river / sea tunneling, the operation of augers in water-rich sandy and gravelly strata is easily affected by pore water pressure disturbances. If local water and soil pressure imbalances occur, it often induces particle migration and detachment, leading to excavation face instability, surface subsidence, and even severe blowout accidents. Current technologies still have the following shortcomings in blowout prediction and control: Traditional empirical formulas are difficult to accurately describe the unstable migration behavior of particle groups under water pressure, resulting in large prediction biases; single CFD or discrete element method (DEM) methods cannot simultaneously simulate large fluid deformation and discrete particle motion, making it difficult to capture key causes of gushing; furthermore, physical model tests are costly and difficult to reproduce complex water and soil loading conditions; and field gushing monitoring technology is lagging behind, making it difficult to identify particle detachment causes and water pressure disturbance responses in a timely manner. Summary of the Invention

[0003] The purpose of this invention is to provide a blowout prevention system and method for a spiral excavator based on the coupling of smooth particle hydrodynamics (SPH) and discrete element method (DEM). This involves constructing a multiphysics coupled model of the spiral excavator under deep burial conditions to simulate the disturbance, migration, and detachment processes of DEM particles under water pressure loading; analyzing the influence of different construction parameters (such as excavation angle and pressure plate loading rate) on particle response patterns in the spiral structure to identify potential triggering factors for blowouts; and proposing a blowout risk criterion based on numerical simulation results to provide a scientific basis for shield machine structural design optimization and water pressure control.

[0004] To achieve the above objectives, the present invention provides a blowout prevention system for a spiral soil excavator based on SPH-DEM coupling, comprising: The SPH fluid simulation module employs the Lagrange particle method to discretize and model the fluid, treating pore water as a collection of discrete particles moving over time. Within the meshless SPH numerical framework, governing equations describing the conservation of mass and momentum in the fluid are introduced to characterize the flow and disturbance behavior of pore water under high water pressure conditions. The Wendland kernel function is used to perform weighted interpolation of physical quantities between adjacent particles, enabling the numerical solution of the governing equations within the particle system. DEM Particle Simulation Module: The discrete element method is used to model soil particles. Soil particles are represented as single rigid spheres and clumps composed of multiple spheres as needed to reflect the morphological characteristics of different particle sizes and irregular particles. Bidirectional coupling interface: used to realize bidirectional exchange of forces between fluid and particles, accelerated by GPU parallel computing, including fluid-particle momentum transfer algorithm, and real-time calculation of drag force, buoyancy and particle disturbance feedback to fluid; Parametric model of the screw conveyor: Construct an adjustable screw conveyor structural model, including variable pitch screw blades, a cylinder with pressure compensation orifice, and an adjustable gate mechanism for the soil outlet; the screw body is used to simulate the actual soil outlet path; the adjustable screw conveyor is equipped with a loading plate assembly at the bottom to realize water-soil coupling loading from the bottom up, simulating the particle group response behavior under high water pressure; Gust warning module: Integrates key physical quantities, including particle detachment rate, local pore pressure growth rate, and relative velocity disturbance amplitude, to construct a gust risk criterion.

[0005] Preferably, in the Lagrange particle method, the interaction between fluid particles is described by a meshless SPH discrete model, where the kernel function is used to define the influence range and weight distribution between adjacent particles, and its mathematical form is as follows: ; in, This is the Wendland kernel function, used to characterize the weights of physical quantity interpolation between adjacent SPH fluid particles; r The distance between the two particles. h To influence the radius; The normalization constant is d For spatial dimensions, The value of varies with the spatial dimension.

[0006] Preferably, the contact force calculation formula in the DEM particle simulation module is as follows: Normal force: ; Tangential force: ; in, F n and F t These are the normal contact force and the tangential contact force, respectively. , These are the overlapping displacements in the normal and tangential directions, respectively. n and t These are the contact normal unit vector and the tangential unit vector, respectively; , These are the components of the relative velocity at the contact point in the normal and tangential directions, respectively. , These are the normal contact stiffness and the tangential contact stiffness, respectively. , These are the normal damping coefficient and the tangential damping coefficient, respectively; Particle size distribution is expressed based on the Rosin-Rammler function, used to generate a mixed particle population. The specific expression is as follows: ;in, Y For particle size smaller d The particle mass fraction; d Particle size; d 0 represents the characteristic particle size; n It is a distribution index used to characterize the width of the particle size distribution.

[0007] Preferably, in the bidirectional coupling interface, the drag force is calculated according to the Schiller-Naumann modified model: ; in, C D The drag coefficient, ρ f For fluid density, v rel The relative velocity between the fluid and the particles. A p This is the reference area for the particles; After particles are subjected to the reaction force of fluid, the resulting reaction force is introduced into the SPH governing equation describing the conservation of fluid momentum in the form of a source term. This achieves the conservation of momentum and bidirectional coupling between the fluid phase and the particle phase without changing the kernel function interpolation form. The surge risk criterion in the surge warning module is a risk index function constructed based on weighted coefficients: ; in, R As an indicator of gushing risk; For each risk factor, specifically, As a normalized index of particle detachment rate; This is a normalized index for the local pore pressure growth rate. This is a normalization index for the amplitude of particle-fluid relative velocity disturbance. ω i For the corresponding weight coefficients, and ω 1+ ω 2+ ω 3 = 1.

[0008] Based on the above system, the present invention also provides a method for preventing blowouts in a spiral soil removal machine based on SPH-DEM coupling, comprising the following steps: S1. Multi-scale model construction: Based on geological exploration data, SPH fluid particles and DEM particle aggregates are generated to construct an adjustable spiral machine three-dimensional model. S2. Coupled Calculation Parameter Settings: Define fluid-particle physical parameters, and set the spiral conveyor rotation boundary and the outlet flow-pressure boundary conditions; S3. Two-way coupled calculation: Start parallel calculation of SPH-DEM to record fluid pressure field, particle velocity distribution and soil outlet flow data in real time; S4. Jet Feature Extraction: Identify the critical state of jetting based on fluid volume fraction and particle-fluid velocity difference, and extract key parameters; S5. Risk Warning and Parameter Optimization: Generates a parameter table for anti-gushing construction through machine learning algorithms, and outputs a control plan when the warning threshold is triggered.

[0009] Preferably, the specific content of S1 is as follows: S11. Based on geological survey data, a water phase particle system is established using the meshless SPH method described by Lagrange. The underground pore water is represented as a discrete fluid particle set moving with time, which is used to simulate the infiltration process of pore water under complex strata conditions and its disturbed propagation behavior. The water phase particles are SPH fluid particles, and their interactions are described by the SPH discrete model. In this model, the Wendland kernel function is introduced to perform weighted interpolation of physical quantities between adjacent fluid particles to meet the smoothness requirements and support domain continuity conditions in numerical computation. The mathematical form of the kernel function is as follows: ; in, This is the Wendland kernel function, used to characterize the weights of physical quantity interpolation between adjacent SPH fluid particles; r The distance between the two particles. h To influence the radius; The normalization constant is d For spatial dimensions, The value of varies with the spatial dimension; When constructing the SPH fluid particle model, the density of saline formation water is corrected using a formula: ; in The density of pure water, The density increase caused by salinity is calculated based on salinity test data; considering the effect of temperature on hydrodynamic viscosity, the Andrade formula is used: ; In the formula For temperature viscosity at that time Reference temperature viscosity at that time It is the activation energy for viscous flow. It is the gas constant; S12. Construct a DEM particle population, using clump aggregate modeling to reflect the discontinuity and contact evolution characteristics of sand and gravel particles, and support the simulation of particle migration, rolling, and peeling. The interparticle contact force is described using the Hertz–Mindlin contact model, and is realized in numerical calculations in the form of equivalent stiffness and damping. The formula for calculating the contact force is as follows: Normal force: ; Tangential force: ,and ; in, F n and F t These are the normal contact force and the tangential contact force, respectively. , These are the overlapping displacements in the normal and tangential directions, respectively. n and t These are the contact normal unit vector and the tangential unit vector, respectively; , These are the components of the relative velocity at the contact point in the normal and tangential directions, respectively. , These are the normal contact stiffness and the tangential contact stiffness, respectively. , These are the normal damping coefficient and the tangential damping coefficient, respectively; The coefficient of friction between particles is given. For sand and gravel mixed formations, the proportion of particles of different sizes is obtained through sieve analysis, and the particle size distribution is determined according to the Rosin–Rammler function. ; in, Y For particle size smaller d The particle mass fraction; d Particle size; d 0 represents the characteristic particle size; n The distribution index is used to characterize the width of the particle size distribution; S13. Import the 3D parametric structural model of the auger soil remover into the solution environment, including the variable pitch auger blades, the cylinder with pressure compensation holes, and the bottom loading plate assembly, which serve as solid boundaries and coupling interfaces. When constructing the 3D model of the auger soil remover, consider the impact of auger blade wear on soil removal efficiency and anti-gushing performance, and adopt a wear model: ; In the formula For wear volume, The wear coefficient is... For normal force, The sliding distance, The value represents the material hardness.

[0010] Preferably, S2 specifically includes the following steps: S21. Set the interaction parameters between particles and fluid, including drag force and buoyancy. The drag force calculation uses the Schiller-Naumann model: ; in, The drag coefficient, ρ f For fluid density, v rel The relative velocity between the fluid and the particles. A p This is the reference area for the particles; drag coefficient The values ​​under different flow regimes are refined based on the particle Reynolds number Re, which is used to characterize the flow state of the relative motion between the fluid and the particles; When Re ≤ 1, the flow is in the low Reynolds number region dominated by viscosity, and the Stokes drag model is adopted: ; when The transition zone adopts ; Turbulent region , ; For non-spherical particles, a shape factor is introduced. Make corrections: ; S22. Define boundary conditions, including the rotation drive of the helical blades and the discharge outlet flow control mechanism, to simulate different construction loads and soil discharge requirements. When setting the rotation boundary of the screw conveyor, an S-curve acceleration / deceleration control is used: ; in, For the rotation boundary of the screw machine at time... Angular acceleration function; The preset maximum angular acceleration is used to limit the peak acceleration during the start-up and braking processes of the screw conveyor. For calculating time; The acceleration phase is the period during which the screw conveyor speed smoothly increases from a standstill. The start time of the constant speed operation phase, in The screw conveyor maintains a constant speed within the designated area; The deceleration phase is the duration during which the screw conveyor speed smoothly decreases from the operating state to the target value. The above acceleration and deceleration control function adopts a cosine S-curve form to ensure that the angular velocity and angular acceleration change continuously during the screw conveyor rotation process, avoiding numerical instability and non-physical impact response caused by sudden loading. The outlet flow-pressure feedback control uses a PID control algorithm: ; in, The control output of the soil outlet flow-pressure feedback control is used to adjust the opening degree of the soil outlet and the equivalent control parameters to achieve dynamic control of the soil outlet process. , , These are the proportional coefficient, integral coefficient, and derivative coefficient, respectively, used to adjust the proportional effect, cumulative effect, and trend suppression effect of the control system on the error response; τ is the integral time variable. System control error, used to characterize the deviation between the target setpoint and the real-time monitored value, is defined as follows, depending on the control objective: and ; in, The target excavation volume is set. For real-time monitoring of the discharge outlet flow rate; The pre-set target excavation pressure, For real-time monitoring of the soil outlet pressure; To control the time variable; S23. Apply a pore water pressure gradient on one side of the soil chamber to reproduce the soil and water loading process under actual water-rich strata conditions; the pore water pressure gradient setting follows Darcy's law: in For seepage velocity, For penetration rate, For viscosity, This represents the pressure gradient.

[0011] Preferably, in S3, the specific content is as follows: S31. Start coupled computing, support GPU parallel solving, and enable adaptive time step algorithm to improve computing stability and efficiency; The formula for updating the position and velocity of SPH particles is: ; ; in, and They represent the firsti The SPH fluid particles in the first n The time step and the first n +1 time step position vector; and These represent the velocity vectors of the particle at the corresponding time steps; For the first i The SPH fluid particles in the first n The acceleration at each time step is calculated by the fluid control equations in the SPH discrete framework, including pressure gradient terms, viscosity terms, and source terms resulting from particle-fluid interaction forces. This is the time step for numerical calculations, used to control the time progression of the SPH particle position and velocity. During the coupled calculation process, the interparticle contact force is updated according to the Hertz–Mindlin contact model in S12; The adaptive time step control strategy is as follows: ;in, This is the updated computation time step; The time step used in the previous calculation step; S32. Synchronously record key data, including fluid pressure field, particle velocity, and soil discharge flow rate, for subsequent extraction and analysis of gushing behavior; The fluid pressure field recorded values ​​are: ; in, Indicates the first n At each time step, located at the spatial index ( i , j , k The fluid pressure value at point () is used to characterize the spatiotemporal distribution characteristics of the SPH fluid pressure field; i , j , k It serves as a spatial index, used to identify different spatial locations within a pressure field. n Number the time steps; The particle velocity vector is represented as: ; in, For the first n The time step m The velocity vector of each particle , , The particles are respectively in x , y , z Velocity component in the direction; The discharge flow rate at the first n The calculation expression for each time step is: ; in, For the first n The outflow rate at each time step is statistically obtained, and the outlet represents the statistical area corresponding to the outflow. For the first n Within a certain time step, the first [time step] passed through the excavation port i Volumetric flux contribution of each particle For time step; S33. When the local fluid volume fraction is detected to exceed the set threshold, the system switches to high-frequency data sampling mode to improve the response accuracy to the precursor of the surge.

[0012] Preferably, the specific content of S4 is as follows: S41. Use three-dimensional data cloud map to identify the initial area of ​​the gushing, and extract the critical parameters related to the occurrence of the gushing, including fluid pressure gradient, soil velocity, and particle disturbance degree. The critical pressure-speed relationship is as follows: ; in, It is a critical parameter for determining the occurrence of gushing, used to characterize the critical characteristics of the system transitioning from a stable soil discharge state to a gushing state under given operating conditions. a , b , c These are the fitting coefficients. The critical relationship is obtained by regression fitting of multiple sets of numerical simulation results to characterize the dimensionless features of the operating conditions and control parameter levels. S42. Extract key physical quantities during a typical eruption evolution process: Local fluid volume fraction peak , is used to characterize the maximum value of fluid enrichment in a local area before and after a gushing event; The local fluid volume fraction is a function of time. For volume fraction to reach peak The corresponding time; extreme values ​​of particle-fluid velocity difference: ; Abrupt change rate of particle flow rate at the soil outlet: ; in, For the first m Each particle velocity vector; For the first i A velocity vector of SPH fluid particles; The relative velocity difference between the particle and the neighboring SPH fluid particles; This represents the maximum particle-fluid velocity difference within the statistical region. and These represent the maximum and minimum particle flow rates at the soil outlet within the statistical time window. The average particle flow rate at the soil outlet within the same time window is used to quantitatively characterize the flow rate fluctuation and abrupt changes during the soil discharge process. S43. Construct a gushing response database, summarize the simulation results under formation conditions and parameter combinations, and provide data support for gushing mechanism analysis and prediction modeling; the database includes the condition number, formation parameters, auger parameters, simulation time, peak volume fraction, extreme velocity difference, flow rate change amplitude, and whether gushing has occurred.

[0013] Preferably, the specific content of S5 is as follows: S51. Establish a surge risk prediction model based on the surge response database, specifically including: Extract time-series data samples under multiple operating conditions from the surge response database, and construct an input sequence to characterize the evolution of the system's operating state in chronological order; The input sequence includes at least formation pore pressure, spiral rotation speed and soil outlet opening, and introduces gushing characteristic parameters based on the gushing mechanism analysis, including volume fraction, particle-fluid velocity difference, and soil outlet flow rate change rate. The input time series data is preprocessed and normalized, and the continuous data is segmented according to the preset time window length to form a sample sequence for predicting the risk of gushing. Based on the above sample sequences, a Long Short-Term Memory (LSTM) neural network is used to learn the temporal characteristics of the gushing evolution process. By training and optimizing the network parameters, a gushing risk prediction model is established, and the gushing risk index and gushing risk level corresponding to the current working conditions are output. S52. Based on the gushing risk prediction model, output the risk index and use construction parameters, including spiral speed and soil discharge opening, as decision variables to reduce the probability of gushing while ensuring soil discharge efficiency. A multi-objective optimization method is used to search for combinations of construction parameters. The multi-objective optimization method is the non-dominated sorting genetic algorithm NSGA-II, which constructs a Pareto optimal parameter solution set that takes into account both safety and soil removal efficiency. S53. When real-time operating conditions are input into the surge risk prediction model, if the output surge risk index exceeds the preset threshold, the system will automatically generate a surge warning report. The early warning report includes at least the risk level, key influencing parameters, and suggested adjustment directions, and feeds back the corresponding dynamic adjustment strategy for construction parameters to the construction control system to guide the real-time control of the auger speed and the opening of the soil discharge.

[0014] Therefore, this invention employs the aforementioned anti-blowout system and method for shield tunneling augers based on SPH-DEM coupling. It uses the meshless SPH method to solve for pore water flow and the DEM method to simulate the dynamic behavior of soil particles, constructing a complete fluid-particle multiphysics coupled computational framework. This method and system are applicable to the prediction and control of blowout risks for shield tunneling augers in deeply buried, water-rich sandy and gravelly strata. It can accurately capture particle migration and blowout triggering mechanisms caused by water pressure disturbances, improving the physical realism and predictive reliability of numerical simulations. Specifically, it includes: (1) Based on the Rocky DEM platform, combined with the SPH fluid dynamics and DEM particle mechanics modules, efficient coupled simulation is achieved through GPU parallel computing. By utilizing the independence and coupling of the physical processes of each module, the computational efficiency and stability of complex water and soil coupled systems are improved, meeting the needs of large-scale engineering simulation. (2) A unified SPH-DEM coupled calculation framework was constructed to reflect the nonlinear influence of underground pore water hydrodynamics on soil particle groups, accurately simulate the detachment, migration and gushing evolution of particle groups, and improve the identification accuracy and early warning capability of gushing risk. (3) It integrates machine learning algorithms to assist in the analysis of gushing risks and the optimization of construction parameters. Combined with an LSTM neural network trained with a large amount of simulation data, it realizes dynamic prediction and real-time early warning of gushing risks, thereby improving the level of construction safety management. (4) It provides a basis for the structural design of the spiral soil removal machine and the control of shield tunneling construction parameters, which helps to optimize soil removal efficiency and reduce the incidence of gushing accidents, and enhance the construction safety and economic benefits of underground engineering.

[0015] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0016] Figure 1 This is a flowchart of an embodiment of the present invention; Figure 2 This is a schematic diagram of the spiral soil excavator tube body according to an embodiment of the present invention; Figure 3 This is a schematic diagram simulating the gushing of the spiral soil excavator according to an embodiment of the present invention; Figure 4 This is a graph showing the change of the mass flow rate at the soil outlet over time according to an embodiment of the present invention. Figure 5 This is a graph showing the cumulative discharge mass over time according to an embodiment of the present invention. Detailed Implementation

[0017] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0018] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0019] Example 1 This invention provides a blowout prevention system for a spiral soil excavator based on SPH-DEM coupling, comprising: The SPH fluid simulation module employs the Lagrange particle method to discretize and model the fluid, treating pore water as a collection of discrete particles moving over time. Within the meshless SPH numerical framework, governing equations describing the conservation of mass and momentum in the fluid are introduced to characterize the flow and disturbance behavior of pore water under high water pressure conditions. The Wendland kernel function is used to perform weighted interpolation of physical quantities between adjacent particles, enabling the numerical solution of the governing equations within the particle system. In the Lagrange particle method, the interaction between fluid particles is described by a meshless SPH discrete model, where the kernel function is used to define the influence range and weight distribution between adjacent particles, and its mathematical form is as follows: ; in, This is the Wendland kernel function, used to characterize the weights of physical quantity interpolation between adjacent SPH fluid particles; r The distance between the two particles. h To influence the radius; The normalization constant is d For spatial dimensions, The value of varies with the spatial dimension.

[0020] The DEM particle simulation module uses the discrete element method to model soil particles. Soil particles are represented as single rigid spheres or clumps composed of multiple spheres, reflecting different particle sizes and irregular particle morphologies. The contact force calculation formula in the DEM particle simulation module is as follows: Normal force: ; Tangential force: ; in, F n and F t These are the normal contact force and the tangential contact force, respectively. , These are the overlapping displacements in the normal and tangential directions, respectively. n and t These are the contact normal unit vector and the tangential unit vector, respectively; , These are the components of the relative velocity at the contact point in the normal and tangential directions, respectively. , These are the normal contact stiffness and the tangential contact stiffness, respectively. , These are the normal damping coefficient and the tangential damping coefficient, respectively; Particle size distribution is expressed based on the Rosin-Rammler function, used to generate a mixed particle population. The specific expression is as follows: ;in, Y For particle size smaller d The particle mass fraction; d Particle size; d 0 represents the characteristic particle size; n It is a distribution index used to characterize the width of the particle size distribution.

[0021] Bidirectional coupling interface: Used to realize bidirectional exchange of forces between fluid and particles, accelerated by GPU parallel computing, including fluid-particle momentum transfer algorithm, and real-time calculation of drag force, buoyancy, and particle disturbance feedback to the fluid; drag force is calculated according to the Schiller-Naumann modified model. ; in, C D The drag coefficient, ρ f For fluid density, v rel The relative velocity between the fluid and the particles. A p This is the reference area for the particles; After particles are subjected to the reaction force of fluid, the resulting reaction force is introduced into the SPH governing equation describing the conservation of fluid momentum in the form of a source term. Without changing the kernel function interpolation form, the momentum conservation and bidirectional coupling between the fluid phase and the particle phase are realized.

[0022] Parametric model of the screw conveyor: Construct an adjustable screw conveyor structural model, including variable pitch screw blades, a cylinder with pressure compensation orifice, and an outlet gate adjustment mechanism; the screw body is used to simulate the actual soil discharge path; the adjustable screw conveyor is equipped with a loading plate assembly at the bottom to realize water-soil coupling loading from the bottom up, simulating the particle group response behavior under high water pressure.

[0023] The surge warning module integrates key physical quantities, including particle detachment rate, local pore pressure growth rate, and relative velocity disturbance amplitude, to construct a surge risk criterion. The surge risk criterion is a risk index function constructed based on weighted coefficients. ; in, R As an indicator of gushing risk; For each risk factor, specifically, As a normalized index of particle detachment rate; This is a normalized index for the local pore pressure growth rate. This is a normalization index for the amplitude of particle-fluid relative velocity disturbance. ω i For the corresponding weight coefficients, and ω 1+ ω 2+ ω 3 = 1.

[0024] The process flow of the anti-blowout method for the auger excavator established using the above system is as follows: Figure 1 As shown, it includes the following steps: S1. Multi-scale model construction: Based on geological exploration data, SPH fluid particles and DEM particle assemblies are generated to construct an adjustable spiral mechanism 3D model; the specific content is as follows: S11. Based on geological survey data, a water phase particle system is established using the meshless SPH method described by Lagrange. The underground pore water is represented as a discrete fluid particle set moving with time, which is used to simulate the infiltration process of pore water under complex strata conditions and its disturbed propagation behavior. Water phase particles are SPH fluid particles, and their interactions are described by the SPH discrete model. In this model, the Wendland kernel function is introduced to perform weighted interpolation of physical quantities between adjacent fluid particles to meet the smoothness requirements and support domain continuity conditions in numerical computation. The mathematical form of the kernel function is: ; in, This is the Wendland kernel function, used to characterize the weights of physical quantity interpolation between adjacent SPH fluid particles; r The distance between the two particles. h To influence the radius; The normalization constant is d For spatial dimensions, The value of varies with the spatial dimension.

[0025] When constructing the SPH fluid particle model, the density of saline formation water is corrected using a formula: ; in The density of pure water, The density increase caused by salinity is calculated based on salinity test data; considering the effect of temperature on hydrodynamic viscosity, the Andrade formula is used: ; In the formula For temperature viscosity at that time Reference temperature viscosity at that time It is the activation energy for viscous flow. is the gas constant.

[0026] S12. Construct a DEM particle population, using clump aggregate modeling to reflect the discontinuity and contact evolution characteristics of sand and gravel particles, and support the simulation of particle migration, rolling, and peeling. The interparticle contact force is described using the Hertz–Mindlin contact model, and is realized in numerical calculations in the form of equivalent stiffness and damping. The formula for calculating the contact force is as follows: Normal force: ; Tangential force: ,and ; in, F n and F t These are the normal contact force and the tangential contact force, respectively. , These are the overlapping displacements in the normal and tangential directions, respectively. n and t These are the contact normal unit vector and the tangential unit vector, respectively; , These are the components of the relative velocity at the contact point in the normal and tangential directions, respectively. , These are the normal contact stiffness and the tangential contact stiffness, respectively. , These are the normal damping coefficient and the tangential damping coefficient, respectively; The coefficient of friction between particles is given. For sand and gravel mixed formations, the proportion of particles of different sizes is obtained through sieve analysis, and the particle size distribution is determined according to the Rosin–Rammler function. ; in, Y For particle size smaller d The particle mass fraction; d Particle size; d 0 represents the characteristic particle size; n The distribution index is used to characterize the width of the particle size distribution; S13. Import the 3D parametric structural model of the auger soil remover into the solution environment, including the variable pitch auger blades, the cylinder with pressure compensation holes, and the bottom loading plate assembly, which serve as solid boundaries and coupling interfaces. When constructing the 3D model of the auger soil remover, consider the impact of auger blade wear on soil removal efficiency and anti-gushing performance, and adopt a wear model: ; In the formula For wear volume, The wear coefficient is... For normal force, The sliding distance, The value represents the material hardness.

[0027] This embodiment takes the calculation of the gushing response of a subway shield tunnel passing through a water-rich sand layer as an example. According to the engineering survey report, the geological conditions of the study area are clearly defined as medium-coarse sand layer with a natural water content of 35% and a dominant particle size d. 70 The thickness is 0.028 mm, and the underground pore water pressure is approximately 0.35 MPa. Based on the above formation parameters, a water-soil two-phase coupled numerical model is constructed, wherein: The aqueous phase was modeled using the meshless SPH method, with approximately 1.2 million SPH fluid particles deployed in the computational domain to simulate the pore water seepage process and the propagation of water pressure disturbances. The soil phase was modeled using the DEM method, with the sandy medium discretized into approximately 80,000 DEM particles to characterize the contact, migration, and detachment behavior between particles.

[0028] In this embodiment, the three-dimensional structural model of the spiral soil discharger is imported into the coupled computational environment. The parameters of the spiral soil discharger are set as follows: screw pitch is 630mm, blade thickness is 30mm, and spiral diameter is 688mm. An adjustable discharge port structure is set to simulate the soil discharge behavior under different discharge opening conditions. The schematic diagram of the spiral soil discharger pipe body is shown below. Figure 2 As shown, the spiral soil excavator participates in the bidirectional coupled calculation of SPH-DEM as a rigid boundary.

[0029] S2. Coupled Calculation Parameter Settings: Define fluid-particle physical parameters, and set the screw conveyor rotation boundary and outlet flow-pressure boundary conditions; specifically including the following steps: S21. Set the interaction parameters between particles and fluid, including drag force and buoyancy. The drag force calculation uses the Schiller-Naumann model: ; in, The drag coefficient, ρ f For fluid density, v rel The relative velocity between the fluid and the particles. A p This is the reference area for the particles; drag coefficient The values ​​for different flow regimes are refined based on the particle Reynolds number Re, which is used to characterize the flow state of the relative motion between the fluid and the particles. When Re ≤ 1, the flow is in the low Reynolds number region dominated by viscosity, and the Stokes drag model is adopted: ; when The transition zone adopts ; Turbulent region , ; For non-spherical particles, a shape factor is introduced. Make corrections: ; S22. Define boundary conditions, including the rotation drive of the helical blades and the discharge outlet flow control mechanism, to simulate different construction loads and soil discharge requirements. When setting the rotation boundary of the screw conveyor, an S-curve acceleration / deceleration control is used: ; in, For the rotation boundary of the screw machine at time... Angular acceleration function; The preset maximum angular acceleration is used to limit the peak acceleration during the start-up and braking processes of the screw conveyor. For calculating time; The acceleration phase is the period during which the screw conveyor speed smoothly increases from a standstill. The start time of the constant speed operation phase, in The screw conveyor maintains a constant speed within the designated area; The deceleration phase is the duration during which the screw conveyor speed smoothly decreases from the operating state to the target value. The above acceleration and deceleration control function adopts a cosine S-curve form to ensure that the angular velocity and angular acceleration change continuously during the screw conveyor rotation process, avoiding numerical instability and non-physical impact response caused by sudden loading. The outlet flow-pressure feedback control uses a PID control algorithm: ; in, The control output of the soil outlet flow-pressure feedback control is used to adjust the opening degree of the soil outlet and the equivalent control parameters to achieve dynamic control of the soil outlet process. , , These are the proportional coefficient, integral coefficient, and derivative coefficient, respectively, used to adjust the proportional effect, cumulative effect, and trend suppression effect of the control system on the error response; τ is the integral time variable. System control error, used to characterize the deviation between the target setpoint and the real-time monitored value, is defined as follows, depending on the control objective: and ; in, The target excavation volume is set. For real-time monitoring of the discharge outlet flow rate; The pre-set target excavation pressure, For real-time monitoring of the soil outlet pressure; To control the time variable; S23. Apply a pore water pressure gradient on one side of the soil chamber to reproduce the soil and water loading process under actual water-rich strata conditions; the pore water pressure gradient setting follows Darcy's law: ; in For seepage velocity, For penetration rate, For viscosity, This represents the pressure gradient.

[0030] In this embodiment, to ensure computational stability, the system coupling time step is set to Δt = 0.08 ms to meet the CFL stability condition of the SPH method and the DEM particle contact stiffness requirement. Under initial conditions, the helical blade speed is set to 13 rpm, the soil outlet opening degree is set to 60%, and the SPH-DEM bidirectional coupling calculation is started.

[0031] S3. Two-way Coupled Calculation: Initiate parallel calculation of SPH-DEM to record fluid pressure field, particle velocity distribution, and outlet flow rate data in real time; specific content includes: S31. Start coupled computing, support GPU parallel solving, and enable adaptive time step algorithm to improve computing stability and efficiency; The formula for updating the position and velocity of SPH particles is: ; ; in, and They represent the first i The SPH fluid particles in the first n The time step and the first n +1 time step position vector; and These represent the velocity vectors of the particle at the corresponding time steps; For the first i The SPH fluid particles in the first n The acceleration at each time step is calculated by the fluid control equations in the SPH discrete framework, including pressure gradient terms, viscosity terms, and source terms resulting from particle-fluid interaction forces. This is the time step for numerical calculations, used to control the time progression of the SPH particle position and velocity. During the coupled calculation process, the interparticle contact force is updated according to the Hertz–Mindlin contact model in S12; The adaptive time step control strategy is as follows: ;in, This is the updated computation time step; The time step used in the previous calculation step; S32. Synchronously record key data, including fluid pressure field, particle velocity, and soil discharge flow rate, for subsequent extraction and analysis of gushing behavior; The fluid pressure field recorded values ​​are: ; in, Indicates the first n At each time step, located at the spatial index ( i , j , k The fluid pressure value at point () is used to characterize the spatiotemporal distribution characteristics of the SPH fluid pressure field; i , j , k It serves as a spatial index, used to identify different spatial locations within a pressure field. n Number the time steps; The particle velocity vector is represented as: ; in, For the first n The time step m The velocity vector of each particle , , The particles are respectively in x , y ,z Velocity component in the direction; The discharge flow rate at the first n The calculation expression for each time step is: ; in, For the first n The outflow rate at each time step is statistically obtained, and the outlet represents the statistical area corresponding to the outflow. For the first n Within a certain time step, the first [time step] passed through the excavation port i Volumetric flux contribution of each particle For time step; S33. When the local fluid volume fraction is detected to exceed the set threshold, the system switches to high-frequency data sampling mode to improve the response accuracy to the precursor of the surge.

[0032] In the bidirectional coupled calculation process, the statistical area of ​​the soil outlet is used as the monitoring section. The particle flux passing through this area in each time step is accumulated to obtain the process data of the change of the soil outlet mass flow rate and the cumulative discharged mass over time.

[0033] Figure 4 The figures show the curves of the mass flow rate at the outlet versus time under different water pressure conditions in this embodiment. The results show that under low water pressure conditions (0.1–0.5 MPa), the mass flow rate at the outlet is mainly distributed in the range of 70–120 t / h during the steady stage, with an average value of about 90–100 t / h, and the flow rate fluctuation is relatively small. When the water pressure is increased to a medium level (0.8–1.6 MPa), the average mass flow rate at the outlet increases to 130–180 t / h, and the instantaneous peak value can reach 200–220 t / h, and the flow rate fluctuation begins to increase significantly.

[0034] Under high water pressure conditions (2.0 MPa and above), the average level of the mass flow rate at the soil outlet has limited further improvement, remaining within the range of 170–200 t / h. However, its instantaneous peak value increases significantly, exceeding 300 t / h at most. Meanwhile, the minimum flow rate can drop to 20–40 t / h in local time periods. The mass flow rate at the soil outlet exhibits obvious strong pulsation characteristics, reflecting a significant increase in the non-stationarity of the soil outlet process time series.

[0035] Based on the time history data of the mass flow rate at the discharge outlet, time integration is performed to obtain the relationship between the cumulative discharged mass and time, such as... Figure 5As shown in the figure. The results indicate that, throughout the entire calculation time range, the cumulative discharged mass under different water pressure conditions generally exhibits an approximately linear increasing trend with time. Specifically, under low water pressure conditions, the cumulative discharged mass within 20 seconds is approximately 0.42–0.45 t; under medium water pressure conditions, the cumulative discharged mass increases to 0.48–0.55 t; while under high water pressure conditions (2.0 MPa and 4.0 MPa), the cumulative discharged mass is approximately 0.58 t and 0.64 t, respectively, with the difference between the two decreasing significantly.

[0036] comprehensive Figure 4 and Figure 5 The results show that as the water pressure level increases, the overall discharge capacity of the system gradually approaches a stable range, and further increasing the water pressure has limited effect on improving the cumulative discharged mass. However, at the same time, the instantaneous fluctuation amplitude and peak response of the discharge outlet mass flow rate are significantly amplified over time, and this type of abnormal fluctuation occurs before the formation of the gushing state. The above-mentioned joint change characteristics of the discharge outlet mass flow rate and the cumulative discharged mass can serve as an important process signal characterizing the abnormal evolution of the discharge process, providing a quantitative basis for the construction of subsequent gushing early warning criteria and the triggering of high-frequency data sampling.

[0037] S4. Gushing Feature Extraction: Identifying the critical state of gushing based on fluid volume fraction and particle-fluid velocity difference, and extracting key parameters; specific content includes: S41. Use three-dimensional data cloud map to identify the initial area of ​​the gushing, and extract the critical parameters related to the occurrence of the gushing, including fluid pressure gradient, soil velocity, and particle disturbance degree. The critical pressure-speed relationship is as follows: ; in, It is a critical parameter for determining the occurrence of gushing, used to characterize the critical characteristics of the system transitioning from a stable soil discharge state to a gushing state under given operating conditions. a , b , c These are the fitting coefficients. The critical relationship is obtained by regression fitting of multiple sets of numerical simulation results to characterize the dimensionless features of the operating conditions and control parameter levels. S42. Extract key physical quantities during a typical eruption evolution process: Local fluid volume fraction peak , is used to characterize the maximum value of fluid enrichment in a local area before and after a gushing event; The local fluid volume fraction is a function of time. For volume fraction to reach peak The corresponding time; extreme values ​​of particle-fluid velocity difference: ; Abrupt change rate of particle flow rate at the soil outlet: ; in, For the first m Each particle velocity vector; For the first i A velocity vector of SPH fluid particles; The relative velocity difference between the particle and the neighboring SPH fluid particles; This represents the maximum particle-fluid velocity difference within the statistical region. and These represent the maximum and minimum particle flow rates at the soil outlet within the statistical time window. The average particle flow rate at the discharge port within the same time window is used to quantitatively characterize the flow rate fluctuations and abrupt changes during the discharge process; a schematic diagram of the auger discharger's jetting simulation is shown below. Figure 3 As shown.

[0038] S43. Construct a gushing response database, summarize the simulation results under formation conditions and parameter combinations, and provide data support for gushing mechanism analysis and prediction modeling; the database includes the condition number, formation parameters, auger parameters, simulation time, peak volume fraction, extreme velocity difference, flow rate change amplitude, and whether gushing has occurred.

[0039] In the gushing feature extraction stage, the gushing state is quantitatively determined based on process response quantities such as the mass flow rate and fluid volume fraction at the S32 outlet.

[0040] Taking the mass flow rate at the excavation outlet as an example, a statistical analysis of its time series is performed, and the average mass flow rate within the statistical window is defined. Maximum mass flow rate and minimum mass flow rate And construct a quality flow rate fluctuation index: ; Simulation results show that during the stable soil excavation stage, the mass flow rate fluctuation index... I m Typically less than 0.3; when obvious abnormal evolutionary characteristics appear during the excavation process, I m It rapidly rises to the 0.5–0.7 range and, in most operating conditions, is satisfied before the gushing state determination conditions are met.

[0041] Therefore, a mass flow rate fluctuation index exceeding 0.5 can be used as one of the important reference conditions for determining the gushing state, and it can be used in conjunction with fluid volume fraction and particle-fluid velocity difference for gushing feature identification.

[0042] S5. Risk warning and parameter optimization: Generate a surge prevention construction parameter table through machine learning algorithms, and output a control plan when the warning threshold is triggered. During the gushing warning and control phase, the soil excavation process is monitored in real time based on the S4 gushing judgment index, and an early warning response is triggered before the judgment conditions are about to be met.

[0043] In this embodiment, under high water pressure conditions, multiple simulation results show that when the mass flow rate fluctuation index... I m When the preset threshold is exceeded for the first time, there is still a time interval of about 5-12 seconds before the system meets the conditions for determining the gushing state.

[0044] Within the analyzed operating conditions, the early warning strategy constructed based on the above indicators has a success rate of over 80% in identifying the gushing state in advance, providing an effective time window for adjusting the operating parameters of the auger and implementing anti-gushing control measures.

[0045] The specific content is as follows: S51. Establish a surge risk prediction model based on the surge response database, specifically including: Extract time-series data samples under multiple operating conditions from the surge response database, and construct an input sequence to characterize the evolution of the system's operating state in chronological order; The input sequence includes at least formation pore pressure, spiral rotation speed and soil outlet opening, and can introduce gushing characteristic parameters such as volume fraction, particle-fluid velocity difference, and soil outlet flow rate change rate according to the gushing mechanism analysis; The input time series data is preprocessed and normalized, and the continuous data is segmented according to the preset time window length to form a sample sequence for predicting the risk of gushing. Based on the above sample sequences, a Long Short-Term Memory (LSTM) neural network is used to learn the temporal characteristics of the gushing evolution process. By training and optimizing the network parameters, a gushing risk prediction model is established, and the gushing risk index and gushing risk level corresponding to the current working conditions are output.

[0046] S52. Based on the risk index output by the gushing risk prediction model, construction parameters such as spiral rotation speed and soil discharge opening are used as decision variables to reduce the probability of gushing while ensuring soil discharge efficiency. A multi-objective optimization method is used to search for combinations of construction parameters. The preferred multi-objective optimization method is the non-dominated sorting genetic algorithm NSGA-II, which constructs a Pareto optimal parameter solution set that balances safety and soil removal efficiency.

[0047] S53. When real-time operating conditions are input into the surge risk prediction model, if the output surge risk index exceeds the preset threshold, the system will automatically generate a surge warning report. The early warning report should include at least the risk level, key influencing parameters and suggested adjustment directions, and feed back the corresponding dynamic adjustment strategy for construction parameters to the construction control system to guide the real-time control of the auger speed and the opening of the soil discharge.

[0048] During the simulation, the fluid volume fraction, particle-fluid relative velocity, and outlet flow rate changes within the soil chamber area were monitored in real time. Simulation results showed that when the calculation time reached approximately 180 seconds, the local fluid volume fraction within the soil chamber exceeded 30%, and the particle-water phase velocity difference reached 1.8 m / s, exceeding the preset safety threshold. This triggered the surge warning module, indicating a high risk of surge under this condition. After the surge risk was triggered, the construction parameters were adjusted, reducing the screw speed to 10 rpm and adjusting the outlet opening to 40%. Coupled simulation calculations were then performed again under the same geological conditions. The optimized calculation results showed that the fluid volume fraction within the soil chamber was stably controlled below 25%, the cooperative motion characteristics of particles and fluid were enhanced, the outlet flow rate fluctuation amplitude was controlled within ±10%, and the overall soil removal process tended to be stable.

[0049] In the actual shield tunneling process of this embodiment, on-site monitoring results showed that the incidence of blowouts decreased from approximately 15% under the original conditions to approximately 3%, and the maximum surface settlement was controlled within ±15mm, meeting the settlement control requirements of the subway project. Further evaluation of the online monitoring and early warning capabilities of this method showed that the SPH-DEM coupling system can output the soil disturbance range and pore water seepage path in real time. Combined with the blowout early warning module, spatial positioning of blowout risks can be achieved with an accuracy of approximately ±0.5m and an early warning response time of less than 2 minutes. Compared with traditional experience-based judgment methods, the accuracy of blowout risk prediction is improved by more than 35%. Simultaneously, using numerical simulation results to guide the adjustment of construction parameters can effectively avoid misjudgments and over-excavation, shortening the actual construction period by approximately 15 days and reducing overall construction costs by approximately 40%.

[0050] Therefore, this invention adopts the above-mentioned anti-blowout system and method for shield tunneling augers based on SPH-DEM coupling. It solves the pore water flow using the meshless SPH method and simulates the dynamic behavior of soil particles using the DEM method, thus constructing a complete fluid-particle multiphysics coupling calculation framework. This method and system are applicable to the prediction and control of blowout risk of shield tunneling augers in deeply buried water-rich sandy and gravelly strata. It can accurately capture the particle migration and blowout triggering mechanism caused by water pressure disturbance, thereby improving the physical realism and prediction reliability of numerical simulation.

[0051] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A blowout-preventing system for a screw-type excavator based on SPH-DEM coupling, characterized in that, Comprise: SPH fluid simulation module: the Lagrangian particle method is used to discretely model the fluid, and the pore water is regarded as a set of discrete particles moving over time; under the numerical framework of meshless SPH, control equations describing the conservation of fluid mass and momentum are introduced to depict the flow and disturbance behavior of pore water under high water pressure conditions; wherein, the physical quantities between adjacent particles are weighted and interpolated by Wendland kernel function, and the numerical solution of the control equation in the particle system is realized; DEM particle simulation module: the discrete element method is used to model the soil particles, and the soil particles are represented as single rigid spheres and clump aggregates composed of multiple spheres according to the needs, so as to reflect the characteristics of different particle sizes and irregular particle shapes; Two-way coupling interface: used to realize the two-way exchange of fluid and particle interaction forces, accelerated by GPU parallel computing, including fluid-particle momentum transfer algorithm, real-time calculation of drag force, buoyancy and disturbance feedback of particles to fluid; Spiral machine parameterized model: a variable pitch spiral blade, a cylinder with pressure compensation hole and a gate adjustment mechanism of the soil outlet are constructed to form a variable spiral machine structure model; the spiral body is used to simulate the actual soil outlet path; the variable spiral machine is provided with a loading plate assembly at the bottom to realize water-soil coupling loading from bottom to top, and to simulate the response behavior of particle groups under high water pressure environment; Gushing early warning module: key physical quantities including particle detachment rate, local pore pressure growth rate and relative velocity disturbance amplitude are fused to construct a gushing risk criterion.

2. The SPH-DEM coupled based anti-blowout system of the screw discharging machine according to claim 1, wherein, In the Lagrangian particle method, the interaction between fluid particles is described by the meshless SPH discrete model, in which the kernel function is used to define the influence range and weight distribution between adjacent particles, and its mathematical form is as follows: ； wherein, is a Wendland kernel function used to characterize the weight of the interpolation of physical quantities between neighboring SPH fluid particles; r is the distance between two particles, h is the influence radius; is a normalization constant, d is the spatial dimension, the value of varies with the spatial dimension.

3. The SPH-DEM coupled based anti-burst system of the screw discharging machine according to claim 1, wherein, The contact force calculation formula in the DEM particle simulation module is as follows: Normal force: ; Tangential force: ; wherein, F n and F t are the normal and tangential contact forces, respectively; , are the normal and tangential overlap displacements, respectively, n and t are the contact normal and tangential unit vectors, respectively; , are the normal and tangential components of the contact point relative velocity; , are the normal and tangential contact stiffnesses; , are the normal and tangential damping coefficients. The particle size distribution is expressed based on a Rosin-Rammler function for generating the mixed particle population, and the specific expression is as follows: ; wherein, Y is the mass fraction of particles with a particle size less than d ; and d is the particle size; d 0 is the characteristic particle size; n is the distribution index for characterizing the width of the particle size distribution.

4. The SPH-DEM coupled based anti-surge system of the screw discharging machine according to claim 1, wherein, In the bidirectional coupling interface, the drag force is calculated according to the Schiller-Naumann correction model: ; wherein, C D is the drag coefficient, After the particles are acted on by the fluid force, the reaction force generated by the particles is introduced into the SPH control equation describing the conservation of fluid momentum in the form of a source term, without changing the interpolation form of the kernel function, so as to realize the momentum conservation and two-way coupling linkage between the fluid phase and the particle phase; f is the fluid density, v rel is the relative velocity between the fluid and the particle, A p is the particle reference area; The gushing risk criterion in the gushing early warning module is a risk index function constructed based on the weighting coefficient: The gushing risk criterion in the gushing early warning module is a risk index function constructed based on the weighting coefficient: ； wherein, R is a gushing risk indicator; is each risk factor, in particular, is a normalized indicator of particle detachment rate; is a normalized indicator of local pore pressure growth rate; is a normalized indicator of particle-fluid relative velocity fluctuation amplitude; The gushing risk criterion in the gushing early warning module is a risk index function constructed based on the weighting coefficient: i is a corresponding weight coefficient, and Comprise the following steps: 1+ S1, multi-scale model construction: generate SPH fluid particles and DEM particle aggregates based on geological survey data, and construct a variable spiral machine three-dimensional model; 2+ S2, coupling calculation parameter setting: define the fluid-particle physical parameters, set the spiral machine rotating boundary and the soil outlet flow-pressure boundary conditions; 3 = 1.

5. The method of claim 1-4, wherein the method is a method of preventing gushing of a spiral excavator based on SPH-DEM coupling. S3, two-way coupling calculation: start SPH-DEM parallel calculation, and record the fluid pressure field, particle velocity distribution and soil outlet flow data in real time; S4, gushing feature extraction: identify the gushing critical state based on the fluid volume fraction and particle-fluid velocity difference, and extract the key parameters; S5, risk early warning and parameter optimization: generate a construction parameter table for preventing gushing through a machine learning algorithm, and output the control scheme when the early warning threshold is triggered. The specific content in S1 is: ​ ​ 6. The method of claim 5, wherein the method is based on SPH-DEM coupling. ​ S11. Based on geological survey data, a water phase particle system is established using the meshless SPH method described by Lagrange. The underground pore water is represented as a discrete fluid particle set moving with time, which is used to simulate the infiltration process of pore water under complex strata conditions and its disturbed propagation behavior. The water phase particles are SPH fluid particles, and their interactions are described by the SPH discrete model. In this model, the Wendland kernel function is introduced to perform weighted interpolation of physical quantities between adjacent fluid particles to meet the smoothness requirements and support domain continuity conditions in numerical computation. The mathematical form of the kernel function is as follows: ； wherein, is a Wendland kernel function used to characterize the weight of the interpolation of physical quantities between neighboring SPH fluid particles; r is the distance between two particles, h is the influence radius; is a normalization constant, d is the spatial dimension, the value of varies with the spatial dimension; When constructing the SPH fluid particle model, the density of saline formation water is corrected using a formula: ； wherein is the pure water density, is the density increment due to salinity, calculated from salinity test data; the effect of temperature on the water dynamic viscosity is taken into account, using the Andrade formula: ； wherein T is the temperature is the viscosity at time t, Trefis the reference temperature is the viscosity at time t, Ea is the viscous flow activation energy, R is the gas constant; S12. Construct a DEM particle population, using clump aggregate modeling to reflect the discontinuity and contact evolution characteristics of sand and gravel particles, and support the simulation of particle migration, rolling, and peeling. The interparticle contact force is described using the Hertz–Mindlin contact model, and is realized in numerical calculations in the form of equivalent stiffness and damping. The formula for calculating the contact force is as follows: Normal force: ; Tangential force: , and ; where, F n and F t are the normal and tangential contact forces, respectively; , are the normal and tangential overlap displacements, respectively, n and t are the contact normal and tangential unit vectors, respectively; , are the normal and tangential components of the relative velocity at the contact point; , are the normal and tangential contact stiffnesses; , are the normal and tangential damping coefficients; is the inter-particle friction coefficient; for sandy-pebbly mixed strata, the proportion of particles of different sizes is obtained through sieve tests, and the particle size distribution is determined according to the Rosin-Rammler function: ； wherein, Y is the mass fraction of particles having a particle size less than d ; d is the particle size of the particles; d 0 is the characteristic particle size; n is the distribution index, which is used to characterize the width of the particle size distribution; S13, the three-dimensional parameter structure model of the screw excavator is introduced into a solving environment, including a variable pitch screw blade, a cylinder with pressure compensation holes, and a bottom loading plate assembly, which are used as a solid boundary and a coupling interface; when the three-dimensional model of the screw excavator is constructed, the influence of the wear of the screw blade on the unearthing efficiency and the anti-surge performance is considered, and a wear model is adopted: ; wherein is the wear volume, is the wear coefficient, is the normal force, is the sliding distance, is the material hardness.

7. The method of claim 6, wherein, S2 specifically includes the following steps: S21, setting the interaction parameters between the particles and the fluid, including the drag force, the buoyancy, the drag force calculation using the Schiller-Naumann model: ; wherein, is the drag coefficient, ρ f is the fluid density, v rel is the relative velocity between the fluid and the particle, A p is the particle reference area; Drag coefficient The values at different flow regimes are refined according to the particle Reynolds number Re, which is used to characterize the flow state of the relative movement between the fluid and the particles. When Re < 1, the flow is in the low Reynolds number regime dominated by viscosity, and the Stokes drag model is used: ; When , the transition zone adopts ; turbulent zone , ; For non-spherical particles, a shape factor is introduced Amendments made: ； S22. Define boundary conditions, including the rotation drive of the helical blades and the discharge outlet flow control mechanism, to simulate different construction loads and soil discharge requirements. When setting the rotation boundary of the screw conveyor, an S-curve acceleration / deceleration control is used: ； wherein, is the angular acceleration function of the screw at time ; is the preset maximum angular acceleration used to limit the acceleration peak value during the screw startup and braking process; is the calculation time; is the acceleration phase duration, during which the screw speed is smoothly increased from static to target value; is the start time of the constant speed running phase, during which the screw maintains constant speed; ; and is the deceleration phase duration, during which the screw speed is smoothly decreased from running state to target value; the above acceleration and deceleration control functions adopt cosine type S-curve form to ensure the continuous change of angular velocity and angular acceleration during the screw rotation process, avoiding numerical instability and non-physical impact response caused by sudden load. The outlet flow-pressure feedback control uses a PID control algorithm: ； wherein, is the control output quantity of the discharge opening flow-pressure feedback control, used to adjust the opening degree of the discharge opening and the equivalent control parameter, so as to realize dynamic control of the unearthing process; , , are respectively proportional coefficient, integral coefficient and differential coefficient, used to adjust the proportional action, cumulative action and change trend suppression action of the control system to the error response; τ is the integral time variable; is the system control error, used to represent the deviation between the target set value and the real-time monitoring value, and its definition is set according to the control target: and ; wherein, is a preset target unearthing flow rate, is a real-time monitored unearthing flow rate; is a preset target unearthing pressure, is a real-time monitored unearthing pressure; is a control time variable; S23. Apply a pore water pressure gradient on one side of the soil chamber to reproduce the soil and water loading process under actual water-rich strata conditions; the pore water pressure gradient setting follows Darcy's law: wherein is the seepage velocity, is the permeability, is the viscosity, is the pressure gradient.

8. The method of claim 7, wherein, In S3, the specific content is as follows: S31. Start coupled computing, support GPU parallel solving, and enable adaptive time step algorithm to improve computing stability and efficiency; The formula for updating the position and velocity of SPH particles is: ； ； wherein, and denote the position vector of the i thSPH fluid particle at the n thtime step and the n +1thtime step, respectively; and denote the velocity vector of the particle at the corresponding time step, respectively; is the acceleration of the i thSPH fluid particle at the n thtime step, which is calculated from the fluid governing equations in the SPH discretization framework, including the pressure gradient term, the viscous term, and the source term due to the particle-fluid interaction force; is the numerical time step length used to control the time advancement of the SPH particle position and velocity. During the coupled calculation process, the interparticle contact force is updated according to the Hertz–Mindlin contact model in S12; The adaptive time step control strategy is: ; wherein, is the updated time step for the calculation; is the time step used for the previous calculation step; S32. Synchronously record key data, including fluid pressure field, particle velocity, and soil discharge flow rate, for subsequent extraction and analysis of gushing behavior; The fluid pressure field recorded values are: ; wherein, represents the fluid pressure value at spatial index ( i , j , k ) at the n th time step, for characterizing the space-time distribution feature of the SPH fluid pressure field; i 、 j 、 k is a spatial index for identifying different spatial positions in the pressure field, n is a time step number; The particle velocity vector is expressed as: ; wherein is the velocity vector of the i-th particle at the j-th time step, n is the velocity vector of the i-th particle at the j-th time step, m is the velocity vector of the i-th particle at the j-th time step, , , is the velocity component of the particle in the x-direction, x , y , z is the velocity component of the particle in the x-direction. The calculation expression of the outflow rate at the nth time step is: n ​ ； wherein, is the outlet flow rate at the time step, n is the outlet flow rate at the time step, is the volume flux contribution of the particle at the time step, n is the volume flux contribution of the particle at the time step, i is the volume flux contribution of the particle at the time step, is the time step. S33. When the local fluid volume fraction is detected to exceed the set threshold, the system switches to high-frequency data sampling mode to improve the response accuracy to the precursor of the surge.

9. The method of claim 8, wherein, The specific content in S4 is as follows: S41. Use three-dimensional data cloud map to identify the initial area of ​​the gushing, and extract the critical parameters related to the occurrence of the gushing, including fluid pressure gradient, soil velocity, and particle disturbance degree. The critical pressure-speed relationship is as follows: ； wherein, is a critical judgment parameter for gushing occurrence, used to represent the critical characteristics of the system transition from stable state to gushing state under given working conditions; a , b , c is a fitting coefficient, is a dimensionless characteristic quantity representing the working condition and control parameter level, and the critical relationship is obtained by regression fitting on multiple sets of numerical simulation results; S42. Extract key physical quantities during a typical eruption evolution process: local fluid volume fraction peak , the maximum value of the fluid enrichment degree in the local area before and after the gushing occurs; is the local fluid volume fraction as a function of time, is the time corresponding to the peak value of the volume fraction; is the particle-fluid velocity difference extreme value; is the particle flow rate mutation rate of the outlet; ; wherein, is the velocity vector of the particle; m is the velocity vector of the SPH fluid particle; is the relative velocity difference between the particle and the adjacent SPH fluid particle; i is the velocity vector of the SPH fluid particle; is the maximum value of the particle-fluid velocity difference within the statistical region; is the maximum value of the particle-fluid velocity difference within the statistical region; and are the maximum and minimum values of the particle flow rate at the discharge port within the statistical time window, respectively, is the average particle flow rate at the discharge port within the same time window, used to quantitatively characterize the fluctuation and mutation characteristics of the discharge process flow. S43. Construct a gushing response database, summarize the simulation results under formation conditions and parameter combinations, and provide data support for gushing mechanism analysis and prediction modeling; the database includes the condition number, formation parameters, auger parameters, simulation time, peak volume fraction, extreme velocity difference, flow rate change amplitude, and whether gushing has occurred.

10. The method of claim 8, wherein, The specific content of S5 is as follows: S51. Establish a surge risk prediction model based on the surge response database, specifically including: Extract time-series data samples under multiple operating conditions from the surge response database, and construct an input sequence to characterize the evolution of the system's operating state in chronological order; The input sequence includes at least formation pore pressure, spiral rotation speed and soil outlet opening, and introduces gushing characteristic parameters based on the gushing mechanism analysis, including volume fraction, particle-fluid velocity difference, and soil outlet flow rate change rate. The input time series data is preprocessed and normalized, and the continuous data is segmented according to the preset time window length to form a sample sequence for predicting the risk of gushing. Based on the above sample sequences, a Long Short-Term Memory (LSTM) neural network is used to learn the temporal characteristics of the gushing evolution process. By training and optimizing the network parameters, a gushing risk prediction model is established, and the gushing risk index and gushing risk level corresponding to the current working conditions are output. S52. Based on the gushing risk prediction model, output the risk index and use construction parameters, including spiral speed and soil discharge opening, as decision variables to reduce the probability of gushing while ensuring soil discharge efficiency. A multi-objective optimization method is used to search for combinations of construction parameters. The multi-objective optimization method is the non-dominated sorting genetic algorithm NSGA-II, which constructs a Pareto optimal parameter solution set that takes into account both safety and soil removal efficiency. S53. When real-time operating conditions are input into the surge risk prediction model, if the output surge risk index exceeds the preset threshold, the system will automatically generate a surge warning report. The early warning report includes at least the risk level, key influencing parameters, and suggested adjustment directions, and feeds back the corresponding dynamic adjustment strategy for construction parameters to the construction control system to guide the real-time control of the auger speed and the opening of the soil discharge.

Images