Microscopic simulation method for soil erosion and pavement subsidence induced by non-pressure pipeline leakage

Abstract

The invention discloses a mesoscopic simulation method for soil erosion and pavement subsidence induced by non-pressure pipeline leakage, and relates to the technical field of municipal engineering and underground structure catastrophe simulation, and the mesoscopic simulation method comprises the following steps: S1, establishing a discrete element particle model; s2, establishing an underground water flow field model; s3, carrying out particle flow-pore seepage coupling solution; s4, constructing a non-pressure pipeline-soil body around the pipeline-pavement structure model; s5, performing model verification through a physical model test; according to the mesoscopic simulation method for soil erosion and pavement subsidence induced by non-pressure pipeline leakage, the problems that in the prior art, microscopic information is difficult to obtain in a physical model test, the efficiency is low, and software incompatibility exists in numerical simulation are solved.

Description

Technical Field

[0001] This invention relates to the field of municipal engineering and underground structure disaster simulation technology, and in particular to a microscopic simulation method for soil erosion and road subsidence induced by unpressurized pipeline leakage. Background Technology

[0002] With the acceleration of urbanization, the scale of underground unpressurized drainage pipe networks is constantly expanding. Pipe leakage caused by factors such as pipe aging, construction defects, geological changes, and external loads is becoming increasingly prominent. As soil and water are lost at the leakage points, soil erosion gradually spreads upwards, eventually forming "erosion voids" beneath the road. When erosion voids reach a certain scale, the load on the overlying road surface will not be adequately supported by the soil. Especially under traffic loads or other external loads, the road may undergo plastic deformation or even fracture, exacerbating the potential risk of road collapse.

[0003] Existing research primarily analyzes underground pipeline leakage through physical model tests and numerical simulations. Physical model tests can directly reflect the formation of erosion cavities and surface subsidence, but they struggle to obtain information on soil particle-level movement, and the testing cycles are long and the types of operating conditions are limited. In terms of numerical simulation, researchers typically use CFD-DEM fluid-structure interaction numerical simulation to study the soil erosion mechanism caused by pipeline leakage. However, this method requires the integration of multiple software programs, which may lead to compatibility issues in the interface function implementations between these programs, resulting in inconsistent data transmission or loss of accuracy, thus affecting the effectiveness of fluid-structure interaction. Summary of the Invention

[0004] The purpose of this invention is to provide a microscopic simulation method for soil erosion and road subsidence induced by leakage from unpressurized pipelines, which solves the problems of existing technologies, such as the difficulty in obtaining microscopic information and low efficiency of physical model tests, and the software incompatibility of numerical simulations.

[0005] To achieve the above objectives, this invention provides a microscopic simulation method for soil erosion and road subsidence induced by leakage in unpressurized pipelines, comprising the following steps: S1. Establish a discrete element particle model; S2. Establish a groundwater flow field model; S3. Perform particle flow-pore flow coupling solution; S4. Construct a structural model of the unpressurized pipeline, surrounding soil, and road surface. S5. Verify the model through physical model experiments.

[0006] Therefore, the present invention employs the above-mentioned microscopic simulation method for soil erosion and road subsidence induced by unpressurized pipeline leakage, which has the following beneficial effects: (1) The coupled solution of particle flow and pore flow is completed within the same program framework, which effectively avoids the data transmission and compatibility problems caused by the multiple software interfaces of traditional CFD-DEM. The calculation process is clear and the stability is good. (2) Achieve integrated microscopic simulation of unpressurized pipeline-soil-road surface within a unified framework; adopt the continuous-discontinuous element method to model the unpressurized drainage pipeline, surrounding soil and road surface structure in the same numerical platform, realize the coupled solution of seepage field, particle migration field and road surface structure, and intuitively reflect the impact of unpressurized pipeline leakage on the stability of the overlying road surface. (3) Detailed description of the entire process of seepage erosion and road subsidence; Through particle-level drag force calculation, bond element failure judgment and particle trajectory marking, detailed description of the soil erosion development process is obtained, and the morphological evolution of erosion cavities and road subsidence is captured simultaneously, providing detailed evidence for disaster mechanism research. (4) Support parameterized assessment and risk assessment of typical working conditions; examine the influence of key factors such as leakage geometry, groundwater level, soil gradation, pipeline burial depth and traffic load on soil erosion and road subsidence through a unified numerical platform system, and provide quantitative basis for risk classification of urban road diseases, selection of treatment timing and optimization of repair scheme.

[0007] 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

[0008] Figure 1 This is an overall flowchart of a microscopic simulation method for soil erosion and road subsidence induced by unpressurized pipeline leakage according to the present invention. Figure 2 This is a numerical analysis model and particle arrangement diagram of pipeline leakage according to an embodiment of the present invention; Figure 3 This is a basic solution flowchart for the particle discrete element method according to an embodiment of the present invention; Figure 4 These are schematic diagrams of the connection key model in different dimensions according to embodiments of the present invention; where (a) is a two-dimensional model and (b) is a three-dimensional model; Figure 5 This is a schematic diagram of the connection key model according to an embodiment of the present invention; Figure 6 The first embodiment of the present invention i The particle and the first j A schematic diagram of the contact forces between particles; where (a) is a diagram of the contact distribution between particles, and (b) is a diagram of the action of microscopic contact forces; Figure 7 This is a flowchart illustrating the seepage solution process according to an embodiment of the present invention. Figure 8 This is a flowchart illustrating the solution process for the coupling of particle flow and pore seepage in an embodiment of the present invention. Figure 9 This is a schematic diagram illustrating the spatial relationship between particles and seepage grid units in an embodiment of the present invention. Figure 10 This is a schematic diagram illustrating the calculation of particle drag force according to an embodiment of the present invention; Figure 11 This is a schematic diagram of a physical model device for drainage pipe leakage according to an embodiment of the present invention; Figure 12 This is a schematic diagram simulating the development process of soil seepage and erosion in an embodiment of the present invention. The unit in the diagram is 10,000 steps. Figure 13 This is a schematic diagram of the soil seepage and erosion development process at different times in the experiment of this embodiment of the invention; Figure 14 This is a statistical comparison of soil seepage in the experiments and simulations in the embodiments of the present invention; Figure 15 This is a cloud map showing the development of road surface stress in an embodiment of the present invention. Detailed Implementation

[0009] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0010] Please see Figure 1 A microscopic simulation method for soil erosion and road subsidence induced by leakage from unpressurized pipelines includes the following steps: S1. Establish a discrete element particle model; S11. The particles are arranged in a staggered and adjacent manner to generate a model field to ensure a more uniform porosity distribution among the particles. Each particle is assigned corresponding physical properties and a corresponding mechanical model is selected. Please see Figure 2 Based on the CDEM platform, a numerical analysis model of pipeline leakage with dimensions of 4.2m × 3.2m was established. The concrete pipeline diameter was 0.8m, the groundwater level depth was 0.4m, and the overburden depth was 1.8m. From top to bottom, the model consisted of an asphalt pavement and a sand layer, with the asphalt pavement and sand particles arranged in a staggered and adjacent manner. In the numerical simulation, the strain of the pipeline due to external stress variations was minimal, so the strain and wall thickness effects during leakage were ignored, and the pipeline model was represented by arc segments. If the sand particle size were generated according to the actual size, the number of particles would reach millions, making computer calculations difficult. Therefore, the model was enlarged by approximately 15 times the actual particle size, resulting in a homogeneous particle size of 1.5cm. The initial leak point was located directly above the pipeline, with a defect width of 12cm, eight times the particle size, ensuring unimpeded leakage of soil particles.

[0011] A standard axle load of 100kN was selected for the dual-wheel set, applied directly above the pipeline. The ground pressure of the vehicle tires was 0.7MPa, and the grounding length was 0.3m.

[0012] The seepage of the foundation soil can be considered as a discrete, unbonded mass. The strength of the soil particles is disregarded; only the resistance to inter-particle movement and unreasonable collisions during transport and sliding are considered. Therefore, the soil particle damping is set to 0.1. For the asphalt pavement, the tensile strength is considered to be 5 × 10⁶ Pa, and the damping is 0.8. Simultaneously, pore seepage is applied to soil units below the groundwater level. Specific structural material parameters and initial pore seepage parameters are shown in Table 1.

[0013] Table 1 Numerical simulation parameter values

[0014] S12. At the beginning of each time step, based on the current particle position and rigid surface position, perform a contact search for each particle and rigid surface to establish particle contact relationships. Please see Figure 4-5 A connection key model is used to construct the contact model between particles. The connection key is similar to the beam element in finite element method, with a certain size and shape, used to transmit the force and torque between two contacting particles. The "key" between two-dimensional contacting particles is a rectangle, one side of which is the sum of the radii of the two contacting particles, and the other side of which is the diameter of the smaller particle. The "key" between three-dimensional contacting particles is a cylinder, the height of which is the sum of the radii of the two contacting particles, and the radius of which is the radius of the smaller particle. In two-dimensional computation, the length of the connecting bond and width The expressions are as follows: ; ; In the formula, Indicates the first The particle and the first The sum of the radii of all particles; Indicates the first The radius of each particle; Indicates the first The radius of the first particle; the radius of the second particle. The particle and the first Each particle consists of two particles in contact with each other; In 3D calculations, the radius of the connecting bond and height The expressions are as follows: ; ; Cross-sectional area of ​​the connecting key The calculation expression is: ; In the formula, This indicates that a two-dimensional calculation is being performed. This indicates that a three-dimensional calculation is being performed. S13. Apply force-displacement constitutive models to each contact and update the contact forces; Please see Figure 3 The particle discrete element method employs an explicit iterative approach, advancing the timeline in small iterations. Contact search is a crucial step in the discrete element solution process, encompassing both coarse and fine detection. Coarse detection, including subspace and dynamic bounding box methods, primarily maps particles and other geometric elements (rods, rigid surfaces, etc.) onto a background lattice based on their spatial positions, providing local particle search targets for fine detection. Fine detection mainly determines whether two geometric elements are in contact, such as whether particles are in contact with each other or with rigid surfaces. Once contact relationships are established between particles, between particles and rigid surfaces, or between particles and rods, contact constitutive models can be added for mechanical calculations.

[0015] Please see Figure 6 By using the elastic modulus and shear modulus of the connecting bond, the required normal stiffness of the particle discrete element is calculated. and tangential stiffness The expressions are as follows: ; ; In the formula, This represents the average elastic modulus of two particles in contact with each other. This represents the average shear modulus of two particles in contact with each other. The connecting key model can be used to calculate the normal force, tangential force, bending moment, and torque between two contacting particles; the normal contact force of the connecting key. and tangential contact force The calculation expressions are as follows: ; ; In the formula, Indicates the current moment; Indicates the time step; This represents the normal relative displacement increment of two particles in contact. This represents the tangential relative displacement increment between two particles in contact. By correcting the normal and tangential contact forces using the brittle Mohr-Coulomb criterion and the maximum tensile stress criterion, we obtain: ; In the formula, This represents the average tensile strength of two particles in contact. , and These represent the tensile strengths of the two particles in contact with each other; This represents the average cohesive force between two particles in contact. , and These represent the cohesive forces of two particles in contact with each other; This represents the average internal friction angle between two particles in contact. , and These represent the internal friction angles of the two particles in contact with each other; Bending moment of connecting key and torque The calculation expressions are as follows: ; ; In the formula, Represents the moment of inertia; Represents the polar moment of inertia; Indicates the increment of the bending angle; Indicates the increment of the torsion angle; , , and satisfy: ; ; ; ; In the formula, and These represent the radii of the two particles in contact with each other; , and Angle vector The three components satisfy ; Represents the corner transformation matrix; Indicates the first angular velocity of each particle; Indicates the first angular velocity of each particle; When the bonding bond model fails, the bonding bonds between particles cannot transmit bending moment and torque, but they can transmit normal contact force and tangential contact force, satisfying the following expression: ; In the formula, Indicates the initial tensile strength; Indicates the angle of internal friction; Indicates initial cohesion; S14. Calculate the translation and rotation of the particles based on the net external force and net external torque acting on them. Based on the fundamental principles of rigid body dynamics, the first The expressions for calculating the translation and rotation of each particle are as follows: ; ; In the formula, Indicates the first The mass of each particle; Indicates the first The displacement of the center of mass of each particle; Indicates the first The net external force on each particle; Indicates the first The moment of inertia of a single particle; Indicates the first The rotational acceleration of each particle; Indicates the first The net external torque of each particle; satisfy: ; In the formula, Indicates the application to the first External force on each particle; Indicates the relationship with the first The number of particles in contact with each other; Indicates the first The particle and the first Contact force between particles; Indicates contact damping force; Represents the global damping force; satisfy: ; In the formula, Indicates the application to the first External torque on each particle; Indicates the first The particle and the first Contact vector between particles; Represents the global damping torque; In the display algorithm, the translational acceleration and rotational acceleration of the particle are first calculated, with the following expressions: ; ; In the formula, Indicates the first The translational acceleration of each particle; Indicates the first The net external force on each particle; Indicates the first The rotational acceleration of each particle; Indicates the first The net external torque of each particle; The velocity and displacement of the particles are calculated using the Euler forward interpolation method, and the calculation expressions are as follows: ; ; In the formula, Indicates the first The first time step Translational displacement of each particle; Indicates the first The first time step Translational displacement of each particle; Indicates the first The first time step The rotational speed of each particle; Indicates the first The first time step The rotational speed of each particle; No. Incremental displacement at each time step and incremental turning angle The calculation expressions are as follows: ; ; No. The first time step Displacement of individual particles New coordinates of the centroid and corners The calculation expressions are as follows: ; ; ; In the formula, Indicates the first The first time step The displacement of each particle; Indicates the first The first time step The centroid coordinates of each particle; Indicates the first The first time step The corner of each particle.

[0016] S2. Establish a groundwater flow field model. Please refer to [link / reference]. Figure 7 ; S21. Construct a groundwater flow field model using the finite element Blkdyn module in CDEM. Generate a triangular mesh using the command flow method and orthogonalize the mesh elements. The mesh elements are the seepage mesh elements. S22. Define the physical properties of the fluid, the initial conditions and boundary conditions for fluid flow; S23. Calculate the seepage velocity and flow rate of the seepage grid cells; Assume that the transport of water in rock and soil follows Darcy's law, expressed as: ; In the formula, Indicates the seepage grid element number 1 Fluid velocity in each direction; The permeability coefficient of a pore flow grid cell is expressed in units of m. 2 / Pa / s; Represents the relative permeability coefficient. , This represents the average saturation of the seepage grid cells. , This indicates the number of nodes that make up the seepage grid cell. Indicates the first The saturation of each node; Indicates the total pressure of the fluid at the node; Indicates the first y Spatial coordinates in one direction; Simplifying Darcy's law formula using Gauss's divergence theorem, we obtain the following expression: ; In the formula, Indicates the volume of the seepage grid cell; This indicates the total number of faces within a percolation grid cell; Represents the total pressure of the fluid at the nodes. In the The average value within each face; Indicates the first The unit outward normal of each facet Component of direction; Indicates the first The area of ​​each face; Based on the flow velocity of the seepage grid cells, calculate the flow velocity applied by the seepage grid cells to the first... The flow of each node is expressed as: ; In the formula, Indicates the first Traffic received by each node; Indicates the relationship with the first The number of element surfaces associated with each node; Represents the velocity vector of the seepage grid cell; Indicates the relationship with the first The node related to the first The unit outward normal of each face; Indicates the first The total number of nodes on each face; If multiple seepage grid cells exist, the velocity and flow rate at the common node are superimposed; let a certain node be... For a common node of a seepage grid cell, the expressions for the average velocity and total flow rate at the common node are as follows: ; ; In the formula, This represents the average flow velocity at the common nodes after the seepage grid cells are stacked. Indicates the first The flow velocity of each seepage grid cell at a common node; This represents the total flow rate at the common node after the seepage grid cells are stacked; Indicates the first The flow rate of each seepage grid cell at a common node; S24, Calculate the total pressure at the node; If a node is a pressure boundary condition applied to it, then the pore water pressure at that node is the pressure given by the external environment; otherwise, the node saturation at the current time step is calculated using the following expression: ; In the formula, Indicates saturation; Indicates the flow boundary; Indicates porosity; Represents the total volume of the nodes; Indicates the calculation time step; like If, then the pore water pressure at that node is set to 0; if The pore water pressure at the nodes can then be calculated. The expression is: ; In the formula, Indicates the bulk modulus of a fluid; Based on the average saturation of the seepage grid cells The total pressure of the fluid at any node of the seepage grid cell The expression is: ; In the formula, Indicates fluid density; , and These are the three components of the global coordinates of a node in a seepage grid cell; , and These are the three components of global gravitational acceleration; S3. Perform a coupled solution for particle flow and pore seepage. Please refer to [link / reference]. Figure 8-10 ; S31. Based on the pore flow calculation, the seepage field is obtained, the drag force at the location of the particle is calculated, and the drag force is applied as an external force to the center of mass of the particle. S311. Determine the location of the particles; When determining whether the centroid of a particle is located within a certain seepage grid cell, let the particle... i The coordinates are , Let be the total number of faces within the seepage mesh element, and let be the outward normal vector of a given face. Face center coordinates Calculate the direction vectors between the particle's center of mass and its face center. The expression is: ; The direction vector is quantized to obtain: ; In the formula, Represents the unit direction vector between the particle's centroid and face center; Calculate the outward normal direction vector and The dot product; if the dot product value of all surfaces is not less than 0, it indicates that the particle is located within the seepage grid cell; the expression for calculating the dot product is: ; In the formula, Represents the dot product of the surfaces; If the centroid of a particle is located inside a pore flow grid cell, a coupling relationship between the particle and the pore flow grid cell is established. The spatial location of the particle and the flow grid cell is quickly searched using the spatial lattice method. A background network is set with the maximum size of the pore flow grid cell as the lattice size, covering all flow grid cell regions. The combination of all flow grid cells forms a flow grid. The flow grid cells are mapped to the background network based on their centroid locations. The background grid cell number of the particle is calculated based on its centroid coordinates. Using the background grid cell containing the particle as the central grid, potential flow grid cells are searched from adjacent cells (9 grids in 2D and 27 grids in 3D, including the central grid cell itself). S312. Calculate the drag force at the location of the particle; Suppose a certain particle Located in a certain seepage grid cell Inside, calculate the seepage grid cells. In particles The interpolated seepage velocity at the centroid is expressed as: ; In the formula, This represents the seepage velocity component at the interpolation point; Indicates the direction of velocity; Indicates the seepage grid cell in the first... Seepage velocity components at each node; Indicates the first Shape functions on each node; Calculate the flow velocity applied to the particles based on the seepage velocity at the interpolation point. The drag force on the surface is expressed as: ; In the formula, This indicates the total drag force. Indicates the drag coefficient; Indicates fluid density; This represents the sum of the flow velocities at the interpolation point; Indicates the cross-sectional area of ​​the particle; and satisfy: ; ; In the formula, This represents the seepage velocity component at interpolation point 1; This represents the seepage velocity component at interpolation point 2; This represents the seepage velocity component at interpolation point 3; Indicates the radius of the connecting key; The drag force components are calculated based on the resultant drag force acting at the particle's center of mass. The expression is as follows: ; In the formula, This represents the components of the drag force of the particle at its center of mass. Adding the drag force component to the external force of the particle, we obtain the resultant force of the particle, expressed as: ; In the formula, This represents the new particle resultant force; This represents the original net force of the particles; S32. Based on the distribution of particles in the pore flow grid cells, calculate the porosity and permeability coefficient to realize the influence of particle transport on the pore flow characteristics. The specific process for calculating porosity and permeability is as follows: It has The particle is located in the pore flow grid cell. Inside the pore flow grid cell, the porosity is calculated. The expression is: ; In the formula, Indicates a certain particle Pore ​​seepage grid unit Internal volume; This represents the volume of a pore flow mesh cell; Calculate the permeability coefficient of the pore flow grid cell. The expression is: ; In the formula, Indicates the permeability coefficient; Indicates the characteristic size of the voids between particles. ; Indicates the dynamic viscosity of a fluid; The porosity values ​​used in two-dimensional and three-dimensional numerical simulation programs differ slightly. When comparing with indoor model tests, the porosity needs to be converted. By correcting for density, the three-dimensional porosity of particles of equal size is converted to two-dimensional porosity. The expression is as follows: ; ; In the formula, Indicates two-dimensional porosity; Indicates three-dimensional porosity; Indicates the relative density of sandy soil; The permeability coefficient of the flow field is determined according to Darcy's law; the Biot coefficient has the following range: ; In the formula, Represents the Biot coefficient; Indicates two-dimensional porosity; S4. Construct a structural model of the unpressurized pipeline, surrounding soil, and road surface. S41. Construct an integrated geometric model of the unpressurized pipeline, surrounding soil, and road surface structure based on the working conditions. S42. Divide the integrated geometric model into regions and clarify the structural layering of the surrounding soil and road in order to assign corresponding material properties and parameters; The method for dividing the integrated geometric model into regions is as follows: the boundary of the pipeline and the integrated geometric model is established as a continuous element, the soil around the pipe and the road surface are established as discrete particle elements, and the groundwater flow field is established as a finite flow element. The finite flow element of the groundwater flow field is consistent with the discrete particle element of the soil around the pipe in dividing the region, so that the groundwater flow field can effectively act on every particle in the soil around the pipe. S43. Based on relevant literature, set the material parameters of the pavement layer, the soil around the pipe, and the unpressurized pipeline, and adjust the material parameters so that the micro parameters can accurately reflect the macroscopic characteristics of the material, and ensure that the particle flow transport solution meets the calculation accuracy requirements, including elastic modulus, Poisson's ratio, density, compaction degree, internal friction angle of the soil, and void ratio. S44. Set the boundary conditions for the integrated geometric model; Set displacement boundary conditions and load boundary conditions according to actual working conditions, including vehicle loads on the road surface and lateral constraints on the soil. Initial conditions are set for the fluid component in the integrated geometric model, including the initial pressure and velocity of the fluid. S45. Solve and calculate the integrated geometric model of the unpressurized pipeline, surrounding soil, and pavement structure, and analyze the microscopic processes of soil erosion and pavement settlement. Please refer to [link / reference needed]. Figure 15 ; During the model calculation process, the integrity of the pipeline is first maintained, and particle stress balance calculations are performed, imposing gravitational acceleration on the particles to allow them to deposit and consolidate under their own weight. Secondly, the seepage force generated by the groundwater flow field is calculated, forming a groundwater stress flow field. Finally, the particle equilibrium stress state is coupled with the groundwater stress flow field, and the leakage defect at the top of the pipeline is opened, forming a numerical simulation model of pipeline leakage-induced soil erosion under the action of groundwater.

[0017] Under traffic loads, road subsidence occurred directly above the leak point, was vertically symmetrical about the leak point, and increased in size with increasing leakage time. During the subsidence process, the road surface fractured under traffic loads, with fractures first appearing at the bottom layer and later at both ends of the surface subsidence area. Horizontal stress development mainly occurred above the leak point, showing near-symmetrical verticality. The surface directly above the leak point exhibited compressive stress, while the area bordering the subgrade soil exhibited tensile stress; the surface on both sides of the subsidence area exhibited tensile stress, while the area bordering the subgrade soil exhibited compressive stress. Under traffic loads, the tensile stress areas of the road surface all reached their ultimate tensile strength and fractured, confirming that the overloaded tensile strength of the bottom layer was the main cause of road fracture. This pattern provides important evidence for the rapid location of pipeline leaks and early warning of road surface defects in engineering projects.

[0018] S5. Validate the model through physical model experiments. Please refer to [link / reference]. Figure 11-14 ; S51. Macro-scale tests are conducted using a physical model device for pipeline leakage and erosion. In the physical model test, various test parameters are precisely controlled to ensure the accuracy and reliability of the test results. The physical model setup consists of a soil sample box, two side water tanks, and pipes. The soil sample box is designed according to a 1:10 scale of the numerical model, with dimensions of 420mm × 300mm × 420mm (length × width × height). The two side water tanks are each 100mm long, and the pipe diameter is 80mm. The model box is constructed of transparent plexiglass, and a camera is placed at the front to record the leakage development process. A permeable filter plate separates the soil sample box from the water tanks, ensuring free water flow between them. The pipes are made of PVC, with circumferential-longitudinal defects measuring 12mm × 50mm to simulate longitudinal cracks at the top of the pipe. The specific defect size and location can be adjusted according to the experimental conditions.

[0019] The experiment employed a semi-structural model approach, with the pipe defect positioned directly against the front of the model box. In this configuration, leakage observed on the front of the model box could be identified as planar erosion. The center of the pipe in the model was 110mm from the bottom edge, and the pipe's soil cover depth was... The water level above the pipe is 150mm. With soil cover depth The parameters are equal, used to simulate the particle loss pattern in water-rich soil. To clearly observe the erosion development pattern of the soil, a thin layer of red soil sample is laid every 50mm on the soil cover over the pipeline. In the model, the soil cover depth Hs of the pipeline is equal to the water level height. The values ​​are equal, both being 150mm, used to simulate the particle loss pattern in water-rich soil.

[0020] The test soil samples were selected from sandy soil. Impurities and silt particles were removed by sieving, ensuring that the sand particle size was evenly distributed within the range of 0.6–0.9 mm, in order to explore the seepage and erosion mechanism of relatively homogeneous soil. The internal friction angle of the soil was used as a macroscopic mechanical index to calibrate the microscopic parameters of the particles, and the shear strength of the soil was determined by direct shear test.

[0021] The experiment was considered complete when a stable erosion pit formed by soil seepage occurred, and the volume of the collected sand-water mixture was measured. and quality , and according to the formula The soil seepage rate was calculated. The increasing trends of soil seepage in the experiment and simulation were similar, both showing an approximately linear increase. At the end of 64 seconds in the experiment, the total soil seepage was 1309 g, with an average seepage rate of 20.46 g / s. Statistical comparison of soil seepage rates in the experiment and simulation was conducted. By comparing and analyzing the relationship between soil erosion development characteristics and time in the experiment and simulation, the 100,000 iterative cycle calculation time step in the simulation corresponds to 4 seconds in the experiment. The consistency between the simulated and experimental soil seepage development indicates that the simulation method has good simulation performance. The high degree of agreement between the simulated and experimental soil seepage erosion development process and characteristics demonstrates that the CDEM fluid-structure interaction numerical simulation method can accurately reflect the development process of soil erosion induced by pipeline leakage.

[0022] S52. Compare the physical model test results with the simulation results to evaluate the accuracy of the integrated geometric model; By comparing the development process of soil erosion induced by pipeline leakage in simulations and experiments, it was found that the erosion development characteristics of the soil exhibit a high degree of similarity, all going through three stages: initial leakage, expansion of void width, and expansion of void depth. In the initial leakage stage, the erosion disturbance occurs within the foundation soil, without affecting the surface or pavement layers; at this stage, no erosion voids have yet appeared. As the erosion disturbance continues to expand upwards, the foundation soil and pavement layers begin to separate, and erosion voids appear beneath the pavement. At this point, soil erosion enters the void width expansion stage, with the width of the erosion voids gradually increasing until it no longer changes significantly. In the void depth expansion stage, the width of the void induced by soil erosion remains relatively stable, while the depth of the erosion voids increases as soil erosion continues. Both experiments and simulations confirm the three stages of soil erosion development, exploring the evolutionary process of soil erosion induced by pipeline leakage.

[0023] Therefore, this invention employs the aforementioned microscopic simulation method for soil erosion and pavement subsidence induced by unpressurized pipeline leakage. It completes the coupled solution of particle flow and pore seepage within the same program framework, resulting in a clear calculation process and good stability. Within this unified framework, it achieves integrated microscopic simulation of the unpressurized pipeline, soil, and pavement, intuitively reflecting the impact of unpressurized pipeline leakage on the stability of the overlying pavement. It microscopically depicts the entire process of leakage erosion and pavement subsidence, providing microscopic evidence for disaster mechanism research. Furthermore, it supports parametric assessment and risk determination for typical working conditions, providing quantitative basis for risk classification of urban road diseases, selection of remediation timing, and optimization of repair schemes.

[0024] 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 microscopic simulation method for soil erosion and road subsidence induced by leakage in unpressurized pipelines, characterized in that, Includes the following steps: S1. Establish a discrete element particle model; S2. Establish a groundwater flow field model; S3. Perform particle flow-pore flow coupling solution; S4. Construct a structural model of the unpressurized pipeline, surrounding soil, and road surface. S5. Verify the model through physical model experiments.

2. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 1, characterized in that, In S1: S11. The particles are arranged in a staggered and adjacent manner to generate a model field, and each particle is assigned physical properties and a mechanical model is selected. S12. At each time step, based on the current particle position and rigid surface position, perform a contact search for each particle and rigid surface to establish particle contact relationships. A connection key model is used to construct the contact model between particles, where the connection key is used to transmit the force and torque between two contacting particles. In two-dimensional computation, the length of the connecting bond and width The expressions are as follows: ； ； In the formula, Indicates the first The particle and the first The sum of the radii of all particles; Indicates the first The radius of each particle; Indicates the first The radius of the first particle; the radius of the second particle. The particle and the first Each particle consists of two particles in contact with each other; In 3D calculations, the radius of the connecting bond and height The expressions are as follows: ； ； Cross-sectional area of ​​the connecting key The calculation expression is: ； In the formula, This indicates that a two-dimensional calculation is being performed. This indicates that a three-dimensional calculation is being performed. S13. Apply force-displacement constitutive models to each contact and update the contact forces; S14. Calculate the translation and rotation of the particle based on the net external force and net external torque acting on it.

3. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 2, characterized in that, In S13: The required normal stiffness of the particle discrete element is calculated by using the elastic modulus and shear modulus of the connecting bond. and tangential stiffness The expressions are as follows: ； ； In the formula, This represents the average elastic modulus of two particles in contact with each other. This represents the average shear modulus of two particles in contact with each other. Normal contact force of the connecting key and tangential contact force The calculation expressions are as follows: ； ； In the formula, Indicates the current moment; Indicates the time step; This represents the normal relative displacement increment of two particles in contact. This represents the tangential relative displacement increment between two particles in contact. By correcting the normal and tangential contact forces using the brittle Mohr-Coulomb criterion and the maximum tensile stress criterion, we obtain: ； In the formula, This represents the average tensile strength of two particles in contact. , and These represent the tensile strengths of the two particles in contact with each other; This represents the average cohesive force between two particles in contact. , and These represent the cohesive forces of two particles in contact with each other; This represents the average internal friction angle between two particles in contact. , and These represent the internal friction angles of the two particles in contact with each other; Bending moment of connecting key and torque The calculation expressions are as follows: ； ； In the formula, Represents the moment of inertia; Represents the polar moment of inertia; Indicates the increment of the bending angle; Indicates the increment of the torsion angle; , , and satisfy: ； ； ； ； In the formula, and These represent the radii of the two particles in contact with each other; , and Angle vector The three components satisfy ; Represents the corner transformation matrix; Indicates the first angular velocity of each particle; Indicates the first angular velocity of each particle; When the bonding bond model fails, the bonding bonds between particles cannot transmit bending moment and torque, but they can transmit normal contact force and tangential contact force, satisfying the following expression: ； In the formula, Indicates the initial tensile strength; Indicates the angle of internal friction; This represents the initial cohesion.

4. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 3, characterized in that, In S14: Based on the fundamental principles of rigid body dynamics, the first The expressions for calculating the translation and rotation of each particle are as follows: ； ； In the formula, Indicates the first The mass of each particle; Indicates the first The displacement of the center of mass of each particle; Indicates the first The net external force on each particle; Indicates the first The moment of inertia of a single particle; Indicates the first The rotational acceleration of each particle; Indicates the first The net external torque of each particle; satisfy: ； In the formula, Indicates the application to the first External force on each particle; Indicates the relationship with the first The number of particles in contact with each other; Indicates the first The particle and the first Contact force between particles; Indicates contact damping force; Represents the global damping force; satisfy: ； In the formula, Indicates the application to the first External torque on each particle; Indicates the first The particle and the first Contact vector between particles; Represents the global damping torque; The expressions for calculating the translational and rotational accelerations of a particle are as follows: ； ； In the formula, Indicates the first The translational acceleration of each particle; Indicates the first The net external force on each particle; Indicates the first The rotational acceleration of each particle; Indicates the first The net external torque of each particle; The velocity and displacement of the particles are calculated using the Euler forward interpolation method, and the calculation expressions are as follows: ； ； In the formula, Indicates the first The first time step Translational displacement of each particle; Indicates the first The first time step Translational displacement of each particle; Indicates the first The first time step The rotational speed of each particle; Indicates the first The first time step The rotational speed of each particle; No. Incremental displacement at each time step and incremental turning angle The calculation expressions are as follows: ； ； No. The first time step Displacement of individual particles New coordinates of the centroid and corners The calculation expressions are as follows: ； ； ； In the formula, Indicates the first The first time step The displacement of each particle; Indicates the first The first time step The centroid coordinates of each particle; Indicates the first The first time step The corner of each particle.

5. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 4, characterized in that, In S2: S21. Construct a groundwater flow field model using the finite element Blkdyn module in CDEM, and generate triangular meshes and orthogonalize the mesh elements using command flow. S22. Define the physical properties of the fluid, the initial conditions and boundary conditions for fluid flow; S23. Calculate the seepage velocity and flow rate of the seepage grid cells; Assume that the transport of water in rock and soil follows Darcy's law, expressed as: ； In the formula, Indicates the seepage grid element number 1 Fluid velocity in each direction; This represents the permeability coefficient of a pore flow mesh cell; Represents the relative permeability coefficient. , This represents the average saturation of the seepage grid cells. , This indicates the number of nodes that make up the seepage grid cell. Indicates the first The saturation of each node; Indicates the total pressure of the fluid at the node; Indicates the first y Spatial coordinates in one direction; Simplifying Darcy's law formula using Gauss's divergence theorem, we obtain the following expression: ； In the formula, Indicates the volume of the seepage grid cell; This indicates the total number of faces within a seepage grid cell; Represents the total pressure of the fluid at the nodes. In the The average value within each face; Indicates the first The unit outward normal of each facet Component of direction; Indicates the first The area of ​​each face; Based on the flow velocity of the seepage grid cells, calculate the flow velocity applied by the seepage grid cells to the first... The flow of each node is expressed as: ； In the formula, Indicates the first Traffic received by each node; Indicates the relationship with the first The number of element surfaces associated with each node; Represents the velocity vector of the seepage grid cell; Indicates the relationship with the first The node related to the first The unit outward normal of each face; Indicates the first The total number of nodes on each face; If multiple seepage grid cells exist, the velocity and flow rate at the common node are superimposed; Let a certain node be For a common node of a seepage grid cell, the expressions for the average velocity and total flow rate at the common node are as follows: ； ； In the formula, This represents the average flow velocity at the common nodes after the seepage grid cells are stacked. Indicates the first The flow velocity of each seepage grid cell at a common node; This represents the total flow rate at the common node after the seepage grid cells are stacked; Indicates the first The flow rate of each seepage grid cell at a common node; S24, Calculate the total pressure at the node; If a node is a pressure boundary condition applied to it, then the pore water pressure at that node is the pressure given by the external environment; otherwise, the node saturation at the current time step is calculated using the following expression: ； In the formula, Indicates saturation; Indicates the flow boundary; Indicates porosity; Represents the total volume of the nodes; Indicates the calculation time step; like If, then the pore water pressure at the node is 0; if The pore water pressure at the nodes can then be calculated. The expression is: ； In the formula, Indicates the bulk modulus of a fluid; Based on the average saturation of the seepage grid cells The total pressure of the fluid at any node of the seepage grid cell The expression is: ； In the formula, Indicates fluid density; , and These are the three components of the global coordinates of a node in a seepage grid cell; , and These are the three components of global gravitational acceleration.

6. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 5, characterized in that, In S3: S31. Based on the pore flow calculation, the seepage field is obtained, the drag force at the location of the particle is calculated, and the drag force is applied as an external force to the center of mass of the particle. S32. Calculate the porosity and permeability coefficient based on the distribution of particles in the pore seepage grid cell.

7. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 6, characterized in that, In S31: S311. Determine the location of the particles; Suppose a certain particle i The coordinates are , Let be the total number of faces within the seepage mesh element, and let be the outward normal vector of a given face. Face center coordinates Calculate the direction vectors between the particle's center of mass and its face center. The expression is: ； The direction vector is quantized to obtain: ； In the formula, Represents the unit direction vector between the particle's centroid and face center; Calculate the outward normal direction vector and The dot product; if the dot product value of all surfaces is not less than 0, it indicates that the particles are located within the seepage grid cells; the expression for calculating the dot product is: ； In the formula, Represents the dot product of the surfaces; If the centroid of a particle is located inside a pore flow grid cell, a coupling relationship between the particle and the pore flow grid cell is established. The spatial position of the particle and the flow grid cell is quickly searched using the spatial lattice method. A background network is set with the maximum size of the pore flow grid cell as the lattice size, and the background network covers the region of the flow grid cell. The flow grid cell is mapped to the background network according to the centroid position of the flow grid cell. The background grid number where the particle is located is calculated according to the particle's centroid coordinates. With the background grid cell where the particle is located as the center grid, potential flow grid cells are searched from adjacent grids. S312. Calculate the drag force at the location of the particle; Suppose a certain particle Located in a certain seepage grid cell Inside, calculate the seepage grid cells. In particles The interpolated seepage velocity at the centroid is expressed as: ； In the formula, This represents the seepage velocity component at the interpolation point; Indicates the direction of velocity; Indicates the seepage grid cell in the first... Seepage velocity components at each node; Indicates the first Shape functions on each node; Calculate the flow velocity applied to the particles based on the seepage velocity at the interpolation point. The drag force on the surface is expressed as: ； In the formula, This indicates the total drag force. Indicates the drag coefficient; Indicates fluid density; This represents the sum of the flow velocities at the interpolation point; Indicates the cross-sectional area of ​​the particle; and satisfy: ； ； In the formula, This represents the seepage velocity component at interpolation point 1; This represents the seepage velocity component at interpolation point 2; This represents the seepage velocity component at interpolation point 3; Indicates the radius of the connecting key; The drag force components are calculated based on the resultant drag force acting at the particle's center of mass. The expression is as follows: ； In the formula, This represents the component of the drag force at the particle's center of mass. Adding the drag force component to the external force of the particle, we obtain the resultant force of the particle, expressed as: ； In the formula, This represents the new particle resultant force; This represents the original net force of the particles.

8. The microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 7, characterized in that, The specific process for calculating porosity and permeability coefficient in S32 is as follows: It has The particle is located in the pore flow grid cell. Inside the pore flow grid cell, the porosity is calculated. The expression is: ； In the formula, Indicates a certain particle Pore ​​seepage grid unit Internal volume; This represents the volume of a pore flow mesh cell; Calculate the permeability coefficient of the pore flow grid cell. The expression is: ； In the formula, Indicates the permeability coefficient; Indicates the characteristic size of the voids between particles. ; Indicates the dynamic viscosity of a fluid; By adjusting the density, the three-dimensional porosity of particles of equal size is converted into two-dimensional porosity, expressed as: ； ； In the formula, Indicates two-dimensional porosity; Indicates three-dimensional porosity; This indicates the relative density of sandy soil.

9. A microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 8, characterized in that, In S4: S41. Construct an integrated geometric model of the unpressurized pipeline, surrounding soil, and road surface structure based on the working conditions. S42. Divide the integrated geometric model into regions and clarify the structural layering of the surrounding soil and the road. The method for dividing the integrated geometric model into regions is as follows: the boundary between the pipeline and the integrated geometric model is established as a continuous element, the soil around the pipeline and the road surface are established as discrete particle elements, and the groundwater flow field is established as a finite flow element with pores. The finite flow element of the groundwater flow field is consistent with the region division of the discrete particle element of the soil around the pipeline, so that the groundwater flow field can effectively act on every particle in the soil around the pipeline. S43. Set the material parameters for the pavement layer, the surrounding soil and the unpressurized pipeline, including elastic modulus, Poisson's ratio, density, compaction degree, internal friction angle and void ratio of the soil. S44. Set the boundary conditions for the integrated geometric model; Set displacement boundary conditions and load boundary conditions according to actual working conditions, including vehicle loads on the road surface and lateral constraints on the soil. Initial conditions are set for the fluid component in the integrated geometric model, including the initial pressure and velocity of the fluid. S45. Solve and calculate the integrated geometric model of the unpressurized pipeline-peripheral soil-road structure to analyze the microscopic process of soil erosion and road subsidence.

10. A microscopic simulation method for soil erosion and road subsidence induced by leakage in an unpressurized pipeline according to claim 9, characterized in that, In S5: S51. Conduct a macroscopic scale test using a physical model device for pipe leakage erosion, and measure the volume of the collected sand-water mixture. and quality The expression for calculating soil seepage is: ； In the formula, This indicates the particle density of sand; Indicates the density of water; S52. Compare the physical model test results with the simulation results to evaluate the accuracy of the integrated geometric model.

Images