Soil and rockfill dam breach development and warehouse-out flow numerical simulation prediction method and system considering overtopping dam break

Abstract

The invention discloses an earth and rockfill dam breach development and warehouse-out flow numerical simulation prediction method and system considering overtopping dam break, and the method comprises the steps: building grid information according to a topographic map and an earth and rockfill dam, carrying out the discretization of the grid information into a particle model needed by SPH, building a three-dimensional model of a soil body and a water body, building a soil body erosion-deformation model through an SPH method, and carrying out the simulation prediction of the breach development and warehouse-out flow of the earth and rockfill dam. A soil-water interaction coupling model is established by adopting two kinds of particles of a soil body and a water body respectively, and strength reduction, caused by water infiltration and fine particle scouring, of the soil body near a soil-water interface is emphatically considered. By researching the shear strength, erosion rate parameters and the like of the soil body, the breach development of overtopping dam break of the earth and rockfill dam and the prediction of the ex-warehouse flow are realized. According to the method, the coupling influence of soil erosion-deformation and soil-water interaction is considered at the same time, the established model and method can well reflect the soil-water interface evolution and landslide process of the earth and rockfill dam under the three-dimensional condition, and the method has the advantage that the breach development of overtopping dam break of the earth and rockfill dam can be predicted with high precision.

Description

Technical Field

[0001] The present invention belongs to the technical field of geotechnical engineering, and in particular relates to a numerical simulation prediction method and system for earth-rock dam breach development and outflow flow considering overtopping and dam breach. Background Art

[0002] Overtopping is one of the main failure modes of earth-rock dams. The prediction of the breach development and outflow of earth-rock dams is an important guarantee for evaluating the safety of the dam and downstream areas. During the breach process, a large amount of uncontrollable floods usually form, flooding the downstream areas. The evolution of the soil-water interface is related to the flood discharge in the downstream area. Its evolution law is a complex problem involving the soil mechanics properties and hydraulic principles of the dam body, and is affected by the collapse of the breach slope and the soil and water loss at its interface. At the same time, the three-dimensional movement process of the water body is also the key to the collapse of the soil slope. By accurately predicting the breach development and outflow of earth-rock dams considering overtopping, technical guidance can be provided for the safety protection design of downstream embankments.

[0003] Early studies focused more on field investigations, laboratory-scale experiments, centrifuge tests, and large-scale field experiments. Experimental studies can truly reflect the complex mechanical behavior of sediment during overtopping, but it is difficult to independently analyze the influence of key factors. Therefore, numerical models play an important role in the risk assessment of overtopping dam breaches. Previous studies have shown that "head" erosion (increase in downstream slope caused by overtopping flow) is the main mechanism affecting the evolution of the soil-water interface, but current numerical simulations are all based on some empirical formulas and physical parameters (such as earth-rock dam size, clay content, flow rate, tailwater depth, etc.) to reconstruct the evolution of the breach, and there is no sophisticated numerical scheme to reproduce this phenomenon. The evolution of the soil-water interface involves a critical soil-water interaction and large plastic deformation. Most models of existing numerical simulation technology recommend the use of non-Newtonian rheological models to simulate the movement of sediment, and recommend the use of additional shear failure criteria to describe the erosion behavior of sediment. A unified framework cannot be used to describe the problems of water flow and soil deformation. Therefore, it is difficult for the above methods to carry out detailed modeling of earth-rock dam overtopping and dam breach, and they cannot consider the different scouring and damage modes in the breach process caused by key factors, and they cannot effectively predict and correctly evaluate the outflow flow. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention proposes a numerical simulation prediction method and system for the breach development and outflow of an earth-rock dam considering overtopping and dam breach. Considering that smooth particle hydrodynamics (SPH) has advantages in simulating free surface flow and soil plastic flow, a smooth particle hydrodynamics method (SPH) is used to simulate the overtopping and dam breach process of the earth-rock dam in three dimensions. The overtopping and dam breach of the earth-rock dam is refinedly modeled through a soil erosion-deformation model and a coupling model of soil-water interaction, and the construction of inlet and outlet boundaries. The limitation of using an additional shear failure criterion to describe the erosion behavior of sediment is overcome, and a unified SPH framework is used to deal with the large deformation problem of water flow and soil. The influence of the overflow water velocity gradient, the viscosity of the fluid near the soil-water interface, and the shear strength of the earth-rock dam material on the soil slope collapse and the outflow is comprehensively considered.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A numerical simulation prediction method for earth-rock dam breach development and outflow flow considering overtopping dam breach includes the following steps:

[0007] Step 1: constructing a three-dimensional geometric model of the earth-rock dam soil body and a three-dimensional geometric model of the water body, and assigning parameters to the three-dimensional geometric model of the earth-rock dam soil body and the three-dimensional geometric model of the water body;

[0008] Step 2: Meshing the three-dimensional geometric model of the earth-rock dam soil and the three-dimensional geometric model of the water body, taking the grid center as the particle coordinate and the grid volume as the particle volume, to generate a soil particle discrete model and a water particle discrete model;

[0009] Step 3: Based on the soil particle discrete model and the water particle discrete model, initialize the soil and water information, set the inlet and outlet boundary conditions, and calculate the water flow velocity and water level at the inflow boundary; wherein setting the inlet and outlet boundary conditions includes: inflow boundary conditions and outflow boundary conditions;

[0010] Step 4: Based on the water velocity and water level of the inflow boundary, a soil erosion-deformation model and a coupled model of soil-water interaction are established. Based on the coupled model, the motion equations of soil particles and the motion equations of water particles are calculated by the SPH method;

[0011] Step 5: Based on the boundary conditions, the water level, the motion equations of the soil particles, the motion equations of the water particles, and the coupled model of soil-water interaction, measure the breach erosion width and soil erosion depth of the dam after the earth slope of the earth-rock dam collapses;

[0012] Step six: Based on the soil erosion-deformation model, the coupled model of soil-water interaction, the motion equations of the soil particles, the motion equations of the water particles and the outflow boundary conditions, measure the outflow flow and water level of the water body causing the overtopping dam break.

[0013] Preferably, the method of setting the inlet and outlet boundary conditions and calculating the water flow velocity and water level height at the inflow boundary includes:

[0014] The fluid velocity on the inflow boundary is obtained from the one-dimensional momentum conservation equation, expressed as:

[0015]

[0016] In the formula, Expressed as a time differential; It is expressed as the partial derivative in the x direction; h is the water level; g is the gravitational acceleration; v is the horizontal velocity; the equation is discretized using the forward difference method as follows: Among them l in is the length of the inflow area, h in is the initial water level of the reservoir in the upstream area; Δt is the time step; g is the gravitational acceleration; h se and v se is the actual water level and average flow rate of the reservoir in the upstream area; v in represents the average fluid velocity on the inflow boundary of the upstream region;

[0017] The water level at the inflow boundary is calculated as:

[0018]

[0019] in, Expressed as time differential; h r Indicates the water level in the reservoir; A r =dV r / dh r Represents the given storage capacity curve V r =h(h r ) or empirical formula at water level h r The reservoir area at that time; Q in is the inflow obtained based on given hydrological data; Q out It is the outflow flow including soil and water breach.

[0020] Preferably, based on the coupling model, the method of calculating the motion equations of soil particles and the motion equations of water particles by the SPH method includes:

[0021] The discrete mass conservation and momentum conservation equations of soil particle a are:

[0022]

[0023] Where a and b represent soil particles a and b respectively; m a ,u a , ρ a are the mass, acceleration, velocity and mass density associated with particle a respectively; represents the differential of the mass density of particle a with respect to time; represents the sum of all neighboring particles in the support domain of particle a; σ ab is the stress tensor at soil particle a including the artificial viscosity tensor; σ ba is the stress tensor at soil particle b including the artificial viscosity tensor; g ab It represents the discrete form of the derivative of the auxiliary kernel function of soil particle a to soil particle b. The specific formula is g ab =g(x a -x b ,l s ); Similarly, g ba It represents the discrete form of the derivative of the auxiliary kernel function of soil particle b to soil particle a. The specific formula is g ba =g(x b -x a ,l s ); g represents the auxiliary kernel function gradient; x a represents the coordinates of particle a; x b represents the coordinates of particle b; l s is the soil particle a supporting domain B m The kernel length of the size; V a and V b are the volumes of soil particles a and b respectively; δ ρ is a constant, generally taken as 0.1; c0 is the soil sound velocity; ψ ab represents the density difference between soil particles a and b; similarly, ψ ba represents the density difference between soil particle b and soil particle a;

[0024] where x ab is the coordinate difference between soil particles a and b, and the specific formula is x ab =x a -x b ; Similarly where x ba is the coordinate difference between soil particle b and particle a. The specific formula is x ba =x b -x a ; The supporting domain B is soil particle a mThe kernel length of size ;

[0025] The discrete mass conservation and momentum conservation equations of water particle A are:

[0026]

[0027] Where A and B represent water particles A and B respectively; m A ,u A , ρ A are the mass, acceleration, velocity, and mass density associated with particle A; represents the differential of the mass density of particle A with respect to time; represents the sum of all neighboring particles in the support domain of particle A; σ AB is the stress tensor at water particle A including the artificial viscosity tensor; σ BA is the stress tensor at water particle B including the artificial viscosity tensor; g AB It represents the discrete form of the derivative of the auxiliary kernel function of water particle A to water particle B. The specific formula is g AB =g(x A -x B ,l s ); Similarly, g BA It represents the discrete form of the derivative of the auxiliary kernel function of water particle B to water particle A. The specific formula is g BA =g(x B -x A ,l s ); g represents the auxiliary kernel function gradient; x A represents the coordinates of particle A; x B represents the coordinates of particle B; l s The water particle A supports domain B. s The kernel length of the size; V A and V B are the volumes of water particles A and B respectively; δ ρ is a constant, usually taken as 0.1; c0 is the speed of sound in water; ψ AB represents the density difference between water particle A and water particle B; similarly, ψ BA It represents the density difference between water particle B and water particle A; where x AB is the coordinate difference between water particle A and water particle B. The specific formula is x AB =x A -x B ; Similarly where x BA is the coordinate difference between water particle B and water particle A. The specific formula is x BA =x B -xA ; The water particle A supports the domain B. s The kernel length of size .

[0028] Preferably, the method for measuring the outflow rate and water level of the water body in the event of an overtopping dam breach comprises:

[0029]

[0030] Where C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflowing particle C.

[0031] The present invention also provides a numerical simulation prediction system for earth-rock dam breach development and outflow flow considering overtopping and dam breach, comprising: a division module, a discrete module, a first calculation module, a second calculation module, a first measurement module and a second measurement module;

[0032] The partitioning module is used to construct a three-dimensional geometric model of the earth body and a three-dimensional geometric model of the water body of the earth-rock dam, and to assign parameters to the three-dimensional geometric model of the earth body and the three-dimensional geometric model of the water body;

[0033] The discrete module is used to mesh the three-dimensional geometric model of the earth-rock dam soil and the three-dimensional geometric model of the water body, taking the grid center as the particle coordinate and the grid volume as the particle volume, to generate a soil particle discrete model and a water particle discrete model;

[0034] The first calculation module is used to initialize the information of the soil and water bodies based on the soil particle discrete model and the water particle discrete model, set the inlet and outlet boundary conditions, and calculate the water flow velocity and water level height at the inflow boundary; wherein setting the inlet and outlet boundary conditions includes: inflow boundary conditions and outflow boundary conditions;

[0035] The second calculation module is used to establish a soil erosion-deformation model and a coupled model of soil-water interaction based on the water flow velocity and water level of the inflow boundary, and calculate the motion equations of soil particles and the motion equations of water particles based on the coupled model by the SPH method;

[0036] The first measurement module is used to measure the breach erosion width and soil erosion depth of the earth-rock dam after the earth slope collapses based on the boundary conditions, the water level, the motion equations of the soil particles, the motion equations of the water particles, and the coupled model of soil-water interaction;

[0037] The second measurement module is used to measure the water outflow and water level height of the overtopping dam burst based on the soil erosion-deformation model, the soil-water interaction coupling model, the motion equations of the soil particles, the motion equations of the water particles and the outflow boundary conditions.

[0038] Preferably, the process of setting the inlet and outlet boundary conditions and calculating the water flow velocity and water level height at the inflow boundary includes:

[0039] The fluid velocity on the inflow boundary is obtained from the one-dimensional momentum conservation equation, expressed as:

[0040]

[0041] In the formula, Expressed as a time differential; It is expressed as the partial derivative in the x direction; h is the water level; g is the gravitational acceleration; v is the horizontal velocity; the equation is discretized using the forward difference method as follows: Among them l in is the length of the inflow area, h in is the initial water level of the reservoir in the upstream area; Δt is the time step; g is the gravitational acceleration; h se and v se is the actual water level and average flow rate of the reservoir in the upstream area; v in represents the average fluid velocity on the inflow boundary of the upstream region;

[0042] The water level at the inflow boundary is calculated as:

[0043]

[0044] in, Expressed as time differential; h r Indicates the water level in the reservoir; A r =dV r / dh r Represents the given storage capacity curve V r =h(h r ) or empirical formula at water level h r The reservoir area at that time; Q in is the inflow obtained based on given hydrological data; Q out It is the outflow flow including soil and water breach.

[0045] Preferably, based on the coupling model, the process of calculating the motion equations of soil particles and the motion equations of water particles by the SPH method includes:

[0046] The discrete mass conservation and momentum conservation equations of soil particle a are:

[0047]

[0048] Where a and b represent soil particles a and b respectively; m a ,u a , ρ a are the mass, acceleration, velocity and mass density associated with particle a respectively; represents the differential of the mass density of particle a with respect to time; represents the sum of all neighboring particles in the support domain of particle a; σ ab is the stress tensor at soil particle a including the artificial viscosity tensor; σ ba is the stress tensor at soil particle b including the artificial viscosity tensor; g ab It represents the discrete form of the derivative of the auxiliary kernel function of soil particle a to soil particle b. The specific formula is g ab =g(x a -x b ,l s ); Similarly, g ba It represents the discrete form of the derivative of the auxiliary kernel function of soil particle b to soil particle a. The specific formula is g ba =g(x b -x a ,l s ); g represents the auxiliary kernel function gradient; x a represents the coordinates of particle a; x b represents the coordinates of particle b; l s is the soil particle a supporting domain B m The kernel length of the size; V a and V b are the volumes of soil particles a and b respectively; δ ρ is a constant, generally taken as 0.1; c0 is the soil sound velocity; ψ ab represents the density difference between soil particles a and b; similarly, ψ ba represents the density difference between soil particle b and soil particle a;

[0049] where x ab is the coordinate difference between soil particles a and b, and the specific formula is x ab =x a -x b ; Similarly where x ba is the coordinate difference between soil particle b and particle a. The specific formula is x ba =x b -x a ; The supporting domain B is soil particle a mThe kernel length of size ;

[0050] The discrete mass conservation and momentum conservation equations of water particle A are:

[0051]

[0052] Where A and B represent water particles A and B respectively; m A ,u A , ρ A are the mass, acceleration, velocity, and mass density associated with particle A; represents the differential of the mass density of particle A with respect to time; represents the sum of all neighboring particles in the support domain of particle A; σ AB is the stress tensor at water particle A including the artificial viscosity tensor; σ BA is the stress tensor at water particle B including the artificial viscosity tensor; g AB It represents the discrete form of the derivative of the auxiliary kernel function of water particle A to water particle B. The specific formula is g AB =g(x A -x B ,l s ); Similarly, g BA It represents the discrete form of the derivative of the auxiliary kernel function of water particle B to water particle A. The specific formula is g BA =g(x B -x A ,l s ); g represents the auxiliary kernel function gradient; x A represents the coordinates of particle A; x B represents the coordinates of particle B; l s The water particle A supports domain B. s The kernel length of the size; V A and V B are the volumes of water particles A and B respectively; δ ρ is a constant, usually taken as 0.1; c0 is the speed of sound in water; ψ AB represents the density difference between water particle A and water particle B; similarly, ψ BA It represents the density difference between water particle B and water particle A; where x AB is the coordinate difference between water particle A and water particle B. The specific formula is x AB =x A -x B ; Similarly where x BA is the coordinate difference between water particle B and water particle A. The specific formula is x BA =x B -xA ; The water particle A supports the domain B. s The kernel length of size .

[0053] Preferably, the process of measuring the outflow flow and water level of the water body causing the dam to overflow and fail includes:

[0054]

[0055] Where C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflowing particle C.

[0056] Compared with the prior art, the present invention has the following beneficial effects:

[0057] The present invention adopts the SPH method to realize the numerical simulation technology of the breach development and outflow of earth-rock dam considering overtopping dam breach, and establishes a refined numerical scheme containing a soil-water interaction coupling model, which can describe the evolution of the soil-water interface under a unified framework. At the same time, the soil erosion-deformation model is considered, and the strength-reduced soil erosion model is used to describe the influence of water infiltration and fine particle scouring effect on the soil. The change of shear strength and friction angle related to the erosion rate parameter affects the soil strength. The influence of shear strength and erosion rate parameters on the breach erosion width, soil erosion depth, water outflow and water level during the overtopping dam breach of the earth-rock dam can be calculated, thereby realizing the simulation of the whole process of overtopping dam breach.

[0058] The present invention not only overcomes the inconvenience of large-scale field experiments in early studies and the difficulty in independently analyzing key factors, but also makes up for the defect that a unified framework cannot be used to describe water flow and soil deformation problems. The established numerical simulation method for earth-rock dam breach development and outflow considering overtopping dam breach can use a unified SPH framework to deal with water flow and soil deformation problems, and conducts refined modeling of earth-rock dam overtopping dam breach. It can effectively measure the breach erosion width, soil erosion depth, water outflow and water level height, and can also perform a large number of numerical simulations at low cost by changing the shear strength and erosion rate parameters. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0060] Figure 1 A flow chart of a numerical simulation prediction method for earth-rock dam breach development and outflow flow considering overtopping dam breach according to an embodiment of the present invention;

[0061] Figure 2 The three-dimensional geometric model of the soil body and the three-dimensional geometric model of the water body of the earth-rock dam according to the embodiment of the present invention;

[0062] Figure 3 A schematic diagram of an inflow boundary condition and an outflow boundary condition of an embodiment of the present invention;

[0063] Figure 4 A schematic diagram of a soil-water interaction coupling model according to an embodiment of the present invention;

[0064] Figure 5 A schematic diagram of a size model of an earth-rock dam and a reservoir according to an embodiment of the present invention;

[0065] Figure 6 A schematic diagram of an elevation-storage capacity curve of an earth-rock dam model reservoir according to an embodiment of the present invention;

[0066] Figure 7 Schematic diagram of the breach development form of the dam axis section of the earth-rock dam according to the embodiment of the present invention; wherein (a) is a schematic diagram of a low-viscosity earth-rock dam, (b) is a schematic diagram of a medium-viscosity earth-rock dam, and (c) is a schematic diagram of a high-viscosity earth-rock dam;

[0067] Figure 8 Schematic diagram of the overtopping and dam breaching morphology and outflow flow of an earth-rock dam according to an embodiment of the present invention; wherein (a) is a schematic diagram of a low-viscosity earth-rock dam, (b) is a schematic diagram of a medium-viscosity earth-rock dam, and (c) is a schematic diagram of a high-viscosity earth-rock dam. DETAILED DESCRIPTION

[0068] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0069] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0070] Embodiment 1

[0071] The purpose of the present invention is to overcome the defects of the prior art, and adopt the SPH method to realize the numerical simulation technology of the breach development and outflow flow of the earth-rock dam considering the overtopping and dam breach, and to carry out fine modeling of the whole process of the earth-rock dam overtopping and dam breach. The evolution of the soil-water interface is described by the coupling model of the soil erosion-deformation model and the soil-water interaction, breaking through the limitation of the traditional use of empirical formulas to describe the evolution of the soil-water interface, and comprehensively considering the influence of different shear strength and erosion rate parameters on the collapse of the soil slope, the breach development and the outflow flow.

[0072] like Figure 1 The flowchart shown is a flowchart of a numerical simulation prediction method for earth-rock dam breach development and outflow considering overtopping dam breach proposed in an embodiment of the present invention. The method comprises the following steps:

[0073] Step 1: Obtain the elevation of the earth-rock dam and the surrounding mountains according to the terrain elevation map and the earth-rock dam design drawing. Construct a 3D geometric model of the earth-rock dam body and the surrounding mountains based on the elevation information, and construct a 3D geometric model of the water body of the reservoir upstream of the earth-rock dam according to the reservoir water level, and assign material parameters to the 3D geometric model of the earth-rock dam soil, the 3D geometric model of the mountain, and the 3D geometric model of the water body; create a 3D geometric model of the earth-rock dam soil according to the actual situation, such as Figure 2 As shown in the figure, the earth-rock dam has a trapezoidal cross-section, and the dam body gradually narrows from the bottom to the top. The height, width, slope, and preset breach shape of the dam are set according to the design requirements. A three-dimensional geometric model of the water body is created according to the actual reservoir conditions. The shape corresponds to the shape of the dam body. The depth and width of the water body depend on the design and water storage capacity of the reservoir. The soil body is assigned a density ρ, as well as material constitutive parameters: elastic modulus E, Poisson's ratio v, internal friction angle Shear strength c, dilatancy angle ψ and erosion rate parameter k e etc.; Set the water density to ρ′;

[0074] Step 2: Based on the above geometric model, the three-dimensional models of soil and water are divided by using hexahedral unit grids. The unit division should be as uniform as possible, and the soil and water bodies should use the same discrete spacing, and the size of the discrete spacing should be less than 1 / 100 of the maximum side length of the solid geometric model to meet the calculation accuracy requirements. Assuming that the particles are three-dimensional spheres, the center of the grid is the particle coordinate, and the grid volume is the particle volume V, a soil particle discrete model and a water particle discrete model are generated (the discrete model refers to the discretization of the original model into a large number of particles);

[0075] like Figure 4 As shown in the figure, the soil and water of the earth-rock dam are discretized into a set of a large number of particles with a certain mass, and the soil particles and water particles are assigned a mass m = ρV. The soil particle set is set as B m and water particle set B s, the collection of soil particles and water particles corresponds to the mass of the dam body, mountain body and water body respectively. The material constitutive and density parameters of the soil and water body assigned in step 1 are input into the corresponding particle discrete model, so that the collection of particle discrete models is equivalent to the original soil and water body models of the earth-rock dam.

[0076] Step 3: Initialize the soil and water information, set the inlet and outlet boundary conditions, and calculate the water flow velocity and water level at the inflow boundary.

[0077] The area where the initial boundary conditions of soil and water are applied is divided into blocks separately, and the structural parts of soil and water are numbered in blocks, and each unit block is defined as a smooth particle hydrodynamic particle. The displacement and velocity boundary conditions can be given by directly assigning the particle position information, and the force boundary conditions are all converted into force density conditions. Set the water body storage capacity curve, that is, give the water body height as h r The volume of water in the reservoir at that time; the water inflow rate is set according to the hydrological data; the three rigid body degrees of freedom (three-directional velocities) of the soil and water bodies are set.

[0078] like Figure 3 As shown in the figure, the inlet and outlet boundary conditions include inflow boundary conditions and outflow boundary conditions. The inflow boundary condition is to pre-set a water inflow buffer zone, the starting position of the buffer zone coordinate is x0, and there is a certain distance from the earth-rock dam. The inflow boundary is set at a width of l in The hollow circle represents the water body virtual particle in the box, and it keeps a constant distance Δx from the nearest water body particle. The water body virtual particle is initially inactive. When the virtual particle is located in the inflow buffer zone, the water body virtual particle will be activated and join the water body particle system at the set inflow speed. The outflow boundary condition is to pre-set a soil and water particle outflow buffer zone, such as Figure 3 As shown, the outflow boundary is set to a width of l out The number and velocity of particles passing through the outflow boundary box are recorded according to the displacement of soil particles and water particles. When soil or water particles leave the outflow buffer zone, the soil or water particles are no longer calculated, and the outflow flow of soil particles and water particles is calculated by the following formula: In the formula, Q out represents the outflow rate of soil particles and water particles; C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflow particle C. The water velocity and water level at the inflow boundary are calculated by the given water body storage capacity curve and the inflow flow. The fluid velocity v on the inflow boundary can be obtained from the one-dimensional momentum conservation equation as follows: Where v is the average horizontal velocity; Expressed as a time differential; It is expressed as the partial derivative in the x direction; h is the water level; g is the gravitational acceleration. The equation can be discretized using the forward difference method as Among them l in is the length of the inflow area, h in is the initial water level of the reservoir in the upstream area; Δt is the time step; g is the gravitational acceleration; h se and v se is the actual water level and average flow rate of the reservoir in the upstream area; v in It represents the average horizontal velocity of the fluid on the inflow boundary of the upstream region. The water level height h at the inflow boundary r It can be calculated by the following formula:

[0079] in Expressed as time differential; h r Indicates the water level in the reservoir; A r =dV r / dh r Represents the given storage capacity curve V r =h(h r ) or empirical formula at water level h r The reservoir area at that time; Q in is the inflow obtained based on given hydrological data; Q out It is the outflow flow including soil and water breach. In the specific calculation process, water particles in the upstream area exceeding the water level are not included in the calculation.

[0080] Step 4: Use the soil erosion-deformation model and the coupled model of soil-water interaction to jointly describe the water flow and soil deformation problems. In the coupled model, the soil particles use the strength-reduced erosion soil model to describe the influence of water infiltration and fine particle scouring on the soil, and the water particles are regarded as weakly compressible viscous fluids. Based on the coupled model, the SPH method is used to calculate the motion equations of soil particles and water particles.

[0081] In the process of overtopping and dam failure of earth-rock dam, the soil erosion-deformation model is adopted. That is, under the scouring of water flow, by considering the erosion effect of water on soil, the erosion rate parameter is introduced to describe the change of soil strength, and based on the Drucker-Prager yield criterion (including elastic modulus E, Poisson's ratio v, shear strength c and internal friction angle), the soil erosion rate parameter is introduced to describe the change of soil strength. The ultimate stability state of the slope of the earth-rock dam is calculated by using the parameters such as pressure, pressure, and pressure (such as pressure, pressure, and pressure), and the deformation state of the soil of the earth-rock dam is obtained. The migration of fine particles in the soil plays a vital role in the soil-water interface during overflow-induced failure. During the overtopping and dam-breaking process, fine particles may be washed into the shallow soil by surface water. The redistribution of particle size leads to changes in the soil structure, resulting in a decrease in soil strength. Therefore, during the erosion process, the effect of water infiltration on the strength of the soil in the saturated period is reflected through the softening of the soil strength. The description of the softening of soil strength can be expressed as the strength reduction of the shear strength and internal friction angle of the soil, and the calculation formula is:

[0082]

[0083] In the formula, c is the cohesive shear strength; c* represents the shear strength of the soil under soil strength softening; is the friction angle; It indicates the internal friction angle of soil under soil strength softening; is the tangent of the friction angle; tan -1 represents the cotangent function; F e It is an erosion parameter that depends on a number of physical parameters, including water velocity, sediment transport, particle size, etc.

[0084]

[0085] In the formula, k e Represents the erosion rate parameter involving some specific physical parameters; t0 represents the time when the erosion process starts.

[0086] During the overtopping and dam failure of an earth-rock dam, the coupled model of soil-water interaction is considered to calculate the water flow and soil deformation. It is assumed that there is no material exchange process between the soil and the water, that is, the two media are discretized into soil particles and water particles respectively. The support domain of the soil particles only contains the soil, and the support domain of the water particles only contains the water. However, both the soil and the water involve large deformation problems, so the interaction between the soil and the water is considered. The coupled model of soil-water interaction is as follows: Figure 4 As shown in the figure, when calculating the motion of soil particles near the soil-water interface, the support domain size is centered on the particle, l s , and the water particles in the support domain are included in the calculation. Similarly, when calculating the motion of water particles near the soil-water interface, the support domain size is centered on the particle, l s A sphere with a radius of , and the soil particles in the support domain are included in the calculation. In summary, the process of the soil-water interaction coupling model is considered, that is, the calculation of water particles includes soil particles, and soil particles include water particles.

[0087] The calculation method for establishing the motion equation of soil particles is:

[0088] According to the symbolic conventions in tensor theory, based on the virtual work principle of the FPM framework and the deformation control equation form of the soil, the discrete mass conservation and momentum conservation equations of each soil particle a are as follows:

[0089]

[0090] Where a and b represent soil particles a and b respectively; m a ,u a , ρ a are the mass, acceleration, velocity and mass density associated with particle a respectively; represents the differential of the mass density of particle a with respect to time; represents the sum of all neighboring particles in the support domain of particle a; σ ab is the stress tensor at soil particle a including the artificial viscosity tensor; σ ba is the stress tensor at soil particle b including the artificial viscosity tensor; g ab It represents the discrete form of the derivative of the auxiliary kernel function of soil particle a to soil particle b. The specific formula is g ab =g(x a -x b ,l s ); Similarly, g ba It represents the discrete form of the derivative of the auxiliary kernel function of soil particle b to soil particle a. The specific formula is g ba =g(x b -x a ,l s ); g represents the auxiliary kernel function gradient; x a represents the coordinates of particle a; x b represents the coordinates of particle b; l s is the soil particle a supporting domain B m The kernel length of the size; V a and V b are the volumes of soil particles a and b respectively; δ ρ is a constant, generally taken as 0.1; c0 is the soil sound velocity; ψ ab represents the density difference between soil particle a and soil particle b; similarly, ψ ba represents the density difference between soil particle b and soil particle a;

[0091] where x ab is the coordinate difference between soil particles a and b, and the specific formula is x ab =x a -x b ; Similarly where x bais the coordinate difference between soil particle b and particle a. The specific formula is x ba =x b -x a ; The supporting domain B is soil particle a m The kernel length of size .

[0092] Since the calculation results of the FPM method often show unstable numerical oscillations, artificial viscosity can be used to increase internal constraints. The stress tensor including the artificial viscosity tensor can be expressed as:

[0093] σ ab =σ a +∏ ab I

[0094] In the formula, σ a is the Cauchy stress of soil particle a; I is the second-order unit tensor; in order to better reflect the influence of soil-water interface viscosity, the artificial viscosity proposed by Morris et al. can produce a better numerical solution, as shown in the following formula:

[0095]

[0096] In the formula, μ a and μ b is the physical viscosity coefficient associated with the corresponding soil particles; l s is the soil particle a supporting domain B m The kernel length of size v ab =v a -v b is the velocity difference between soil particles a and b; x ab =x a -x b is the coordinate difference between soil particles a and b.

[0097] In order to obtain more accurate prediction results, the FPM (finite particle method) kernel function h ab =h(x a -x b ,l s ) and its derivative g ab =g(x a -x b ,l s ) is based on the coordinate information of the surrounding particles (x ba =[x b -x a ,y b -y a ,z b -z a ] T , where x a、x b Respectively represent the x-axis coordinates of soil particles a and b, and y a ,y b Respectively represent the y-axis coordinates of soil particles a and b, and z a 、z b Respectively represent the z-axis coordinates of soil particles a and b. ) is solved and obtained as follows:

[0098]

[0099] Among them, h ab is the auxiliary kernel function; g ab =[g x,ab g y,ab g z,ab ] T They represent the discrete forms of the derivative of the auxiliary kernel function in the x, y, and z directions respectively; a and b represent soil particles a and b respectively; b is the support domain B of soil particle a m The number of neighboring particles within V b represents the particle volume associated with soil particle b; W ab , is the discrete form of the original kernel function and the three directional derivatives, and the Gaussian function is used as the SPH kernel function. The specific expression is:

[0100]

[0101] In the formula, l s is the soil particle a supporting domain B m The kernel length of the size of x a and x b are the coordinates of soil particles a and b, respectively.

[0102] Assuming that the constitutive relation of the soil conforms to the Drucker-Prager constitutive model, according to the objective strain relationship and elastic-plastic theory, the constitutive relation expression is as follows:

[0103]

[0104] In the formula, is the objective strain rate of the Cauchy stress; D represents the fourth-order elastic tensor; m represents the predefined plastic potential function stress gradient; represents the gradient operator; Indicates speed; It represents the plastic multiplier in the elastic-plastic theory and is determined by the consistency condition of the yield function f, as follows:

[0105] f≤0,

[0106] The yield function and plastic potential function expressions of the Drucker-Prager constitutive model are:

[0107]

[0108] In the formula, c* represents the shear strength of soil under soil strength softening; represents the internal friction angle of soil under soil strength softening; p = -tr(σ) / 3 represents the hydrostatic pressure; s = σ-pI represents the deviatoric stress; tr(·) represents the trace of the tensor; ψ represents the birefringence angle, and tanψ represents its tangent. For plastic flow in large deformation problems, ψ is usually set to 0.

[0109] From the consistency condition, the plastic strain factor It is expressed as:

[0110]

[0111] Where K represents bulk modulus; G represents shear modulus; represents the internal friction angle of soil under soil strength softening; s = σ-pI represents the deviatoric stress; represents the strain rate; tr(·) represents the trace of the tensor; ψ represents the birefringence angle, and tanψ represents its tangent. For plastic flow in large deformation problems, ψ is usually set to 0.

[0112] The calculation method for establishing the motion equation of water particles is:

[0113] According to the symbolic conventions in tensor theory, based on the virtual work principle of the FPM framework and the deformation control equation form of the water body, the discrete mass conservation and momentum conservation equations of each water particle A are as follows:

[0114]

[0115] Where A and B represent water particles A and B respectively; m A ,u A , ρ A are the mass, acceleration, velocity, and mass density associated with particle A; represents the differential of the mass density of particle A with respect to time; represents the sum of all neighboring particles in the support domain of particle A; σ AB is the stress tensor at water particle A including the artificial viscosity tensor; σ BA is the stress tensor at water particle B including the artificial viscosity tensor; g AB It represents the discrete form of the derivative of the auxiliary kernel function of water particle A to water particle B. The specific formula is g AB =g(xA -x B ,l s ); Similarly, g BA It represents the discrete form of the derivative of the auxiliary kernel function of water particle B to water particle A. The specific formula is g BA =g(x B -x A ,l s ); g represents the auxiliary kernel function gradient; x A represents the coordinates of particle A; x B represents the coordinates of particle B; l s The water particle A supports domain B. s The kernel length of the size; V A and V B are the volumes of water particles A and B respectively; δ ρ is a constant, usually taken as 0.1; c0 is the speed of sound in water; ψ AB represents the density difference between water particle A and water particle B; similarly, ψ BA It represents the density difference between water particle B and water particle A; where x AB is the coordinate difference between water particle A and water particle B. The specific formula is x AB =x A -x B ; Similarly where x BA is the coordinate difference between water particle B and water particle A. The specific formula is x BA =x B -x A ; The water particle A supports the domain B. s The kernel length of size .

[0116] Since the calculation results of the FPM method often show unstable numerical oscillations, artificial viscosity can be used to increase internal constraints. The stress tensor including the artificial viscosity tensor can be expressed as:

[0117] σ AB =p A I+∏ AB I

[0118] In the formula, σ A is the Cauchy stress of water particle A; I is the second-order unit tensor; in order to better reflect the influence of soil-water interface viscosity, the artificial viscosity proposed by Morris et al. can produce a better numerical solution, as shown in the following formula:

[0119]

[0120] In the formula, μ A and μB is the physical viscosity coefficient associated with the corresponding water particles; l s The water particle A supports domain B. s The kernel length of size v AB =v A -v B is the velocity difference between water particle A and water particle B; x AB =x A -x B is the coordinate difference between water particle A and particle B.

[0121] In order to obtain more accurate prediction results, the FPM (finite particle method) kernel function h AB =h(x A -x B ,l s ) and its derivative g AB =g(x A -x B ,l s ) is based on the coordinate information of the surrounding particles (x BA =[x B -x A ,y B -y A ,z B -z A ] T , where x A 、x B Respectively represent the x-axis coordinates of water particles A and B, and y-axis coordinates of water particles A and B. A ,y B Respectively represent the y-axis coordinates of water particles A and B, and z A 、z B Respectively represent the z-axis coordinates of water particles A and B. ) is solved, and the specific formula is as follows:

[0122]

[0123] Among them, h AB is the auxiliary kernel function; g AB =[g x,AB g y,AB g z,AB ] T Respectively represent the discrete forms of the derivative of the auxiliary kernel function in the x, y, and z directions; A and B represent water particles A and B respectively; B is the support domain B of water particle A s The number of neighboring particles within V B W represents the particle volume associated with water particle B; AB , is the discrete form of the original kernel function and the three directional derivatives, and the Gaussian function is used as the SPH kernel function. The specific expression is:

[0124]

[0125] In the formula, l s The water particle A supports domain B. s The kernel length of the size of x A and x B are the coordinates of water particles A and B respectively.

[0126] The water body is considered as a weakly compressible viscous fluid characterized by the following constitutive relations:

[0127]

[0128] Where p is the hydrostatic pressure; β is an empirical parameter, usually taken as 7; ρ0 is the initial density of water particles; ρ is the current density of water particles and c is the speed of sound, which is defined as Among them, g is the acceleration due to gravity and h is the initial water level.

[0129] Step 5: Based on the boundary conditions, water level, motion equations of soil particles, motion equations of water particles and soil-water coupling model, record the breach erosion width and soil erosion depth of the dam after the earth slope of the earth-rock dam collapses.

[0130] Based on the set water velocity and water level, soil erosion-deformation model and coupled model of soil-water interaction, the motion equation of soil particles and the motion equation of water particles, the breach development process can be obtained as the soil particles at the top of the earth-rock dam move, thereby obtaining the process of soil-water interface evolution. Standard post-processing software is used to visualize the output of soil slope collapse process and breach development process, and the section in the middle of the earth-rock dam is selected to observe the final breach damage form. The length measurement tool in the post-processing software is used to measure the erosion width of the breach and the soil erosion depth after the collapse of the earth slope of the earth-rock dam.

[0131] Step 6: Based on the soil erosion-deformation model and the coupled model of soil-water interaction, the motion equations of soil particles, the motion equations of water particles and the outflow boundary conditions, record the outflow flow and water level of the water body that will cause the dam to overflow and break.

[0132] Based on the soil erosion-deformation model and the coupled model of soil-water interaction, the motion equations of soil particles, the motion equations of water particles and the outflow boundary conditions, the displacement coordinates of soil particles and water particles at the outflow boundary box are output in the result file to obtain the height of the water level after the overtopping dam burst. The calculation formula for the water outflow flow of the overtopping dam burst is:

[0133]

[0134] Where C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflowing particle C.

[0135] Step 7: Calculate according to the control equations in steps 5, 6, and 7 and the initial conditions in step 4. Given a time step, when the soil slope collapse stops, it is considered that the breach development and discharge process of the earth-rock dam considering overtopping dam breach is completed. Here is a test example to show the effect of the method proposed in this paper. Figure 5 As shown, a simplified diagram of the geometric dimensions of the earth-rock dam and water body is shown, in which the width of the earth-rock dam is 2.6m, the initial height is 0.6m, and the width of the dam crest is 2.5m; there is a gap in the middle of the dam crest with a height and width of 0.05m×0.5m. Figure 6 The storage capacity curve of the upstream reservoir is shown. The entire earth-rock dam is a homogeneous dam, and the water inflow adopts a constant flow rate Q in =0.0245m 3 / s; soil density is 2000kg / m; fluid density is 1000kg / m; soil particles and water particles adopt uniform discrete model, particle spacing and smooth length are both 0.01m, a total of about 4658150 particles. In order to compare the breach development of homogeneous dams under different viscosities, three material parameters with different viscosities are used, and the detailed parameters are shown in Table 1. The calculation results are shown in Figure 7 and Figure 8 shown.

[0136] Table 1 Material parameters

[0137]

[0138] This method simulates the entire process of breach development and outflow rate of an earth-rock dam considering overtopping and dam failure, and comprehensively considers the main factors of breach development such as soil strength softening, soil-water interaction, shear strength, and erosion rate parameters. This method can realize the refined modeling of overtopping and dam failure of earth-rock dams, and can effectively predict the breach development process, outflow rate and slope instability of earth-rock dams. It can then independently analyze the influencing parameters of overtopping and dam failure of earth-rock dams and make a correct assessment of dam safety.

[0139] Embodiment 2

[0140] The present invention also provides a numerical simulation prediction system for earth-rock dam breach development and outflow flow considering overtopping and dam breach, comprising: a division module, a discrete module, a first calculation module, a second calculation module, a first measurement module and a second measurement module;

[0141] The partitioning module is used to construct a three-dimensional geometric model of the earth-rock dam soil body and a three-dimensional geometric model of the water body, and to assign parameters to the three-dimensional geometric model of the earth-rock dam soil body and the three-dimensional geometric model of the water body, as well as to perform mesh partitioning;

[0142] The discrete module is used to generate soil particle discrete models and water particle discrete models with the grid center as the particle coordinate and the grid volume as the particle volume;

[0143] The first calculation module is used to initialize the information of the soil and water bodies based on the soil particle discrete model and the water particle discrete model, set the inlet and outlet boundary conditions, and calculate the water flow velocity and water level height at the inflow boundary; wherein setting the inlet and outlet boundary conditions includes: inflow boundary conditions and outflow boundary conditions;

[0144] The second calculation module is used to establish a soil erosion-deformation model and a coupled model of soil-water interaction based on the water flow velocity and water level of the inflow boundary, and calculate the motion equations of soil particles and the motion equations of water particles based on the coupled model by the SPH method;

[0145] The first measurement module is used to measure the breach erosion width and soil erosion depth of the earth-rock dam after the earth slope collapses based on the boundary conditions, the water level, the motion equations of the soil particles, the motion equations of the water particles, and the coupled model of soil-water interaction;

[0146] The second measurement module is used to measure the outflow rate and water level of the water body causing the overtopping dam break based on the soil erosion-deformation model, the coupled model of soil-water interaction, the motion equations of the soil particles, the motion equations of the water particles and the outflow boundary conditions.

[0147] In this embodiment, the process of setting the inlet and outlet boundary conditions and calculating the water flow velocity and water level height at the inflow boundary includes:

[0148] The fluid velocity on the inflow boundary is obtained from the one-dimensional momentum conservation equation, expressed as:

[0149]

[0150] In the formula, Expressed as a time differential; It is expressed as the partial derivative in the x direction; h is the water level; g is the gravitational acceleration; v is the horizontal velocity; the equation is discretized using the forward difference method as follows: Among them l in is the length of the inflow area, h inis the initial water level of the reservoir in the upstream area; Δt is the time step; g is the gravitational acceleration; h se and v se is the actual water level and average flow rate of the reservoir in the upstream area; v in represents the average fluid velocity on the inflow boundary of the upstream region;

[0151] The water level at the inflow boundary is calculated as:

[0152]

[0153] in, Expressed as time differential; h r Indicates the water level in the reservoir; A r =dV r / dh r Represents the given storage capacity curve V r =h(h r ) or empirical formula at water level h r The reservoir area at that time; Q in is the inflow obtained based on given hydrological data; Q out It is the outflow flow including soil and water breach.

[0154] In this embodiment, based on the coupling model, the process of calculating the motion equations of soil particles and the motion equations of water particles by the SPH method includes:

[0155] The discrete mass conservation and momentum conservation equations of soil particle a are:

[0156]

[0157] Where a and b represent soil particles a and b respectively; m a ,u a , ρ a are the mass, acceleration, velocity and mass density associated with particle a respectively; represents the differential of the mass density of particle a with respect to time; represents the sum of all neighboring particles in the support domain of particle a; σ ab is the stress tensor including the artificial viscosity tensor; g ab It represents the discrete form of the derivative of the auxiliary kernel function of soil particle a to soil particle b. The specific formula is g ab =g(x a -x b ,l s ); g represents the auxiliary kernel function gradient; x a represents the coordinates of particle a; x brepresents the coordinates of particle b; l s is the soil particle a supporting domain B m The kernel length of the size; V a and V b are the volumes of soil particles a and b respectively; δ ρ is a constant, generally taken as 0.1; c0 is the soil sound velocity; ψ ab represents the density difference between soil particles a and b; where x ab is the coordinate difference between soil particles a and b, and the specific formula is x ab =x a -x b ; The supporting domain B is soil particle a m The kernel length of size ;

[0158] The discrete mass conservation and momentum conservation equations of water particle A are:

[0159]

[0160] Where A and B represent water particles A and B respectively; m A ,u A , ρ A are the mass, acceleration, velocity, and mass density associated with particle A; represents the differential of the mass density of particle A with respect to time; represents the sum of all neighboring particles in the support domain of particle A; σ AB is the stress tensor including the artificial viscosity tensor; g AB It represents the discrete form of the derivative of the auxiliary kernel function of water particle A to water particle B. The specific formula is g AB =g(x A -x B ,l s ); g represents the auxiliary kernel function gradient; x A represents the coordinates of particle A; x B represents the coordinates of particle B; l s The water particle A supports domain B. s The kernel length of the size; V A and V B are the volumes of water particles A and B respectively; δ ρ is a constant, usually taken as 0.1; c0 is the speed of sound in water; ψ AB It represents the density difference between water particle A and water particle B; where x AB is the coordinate difference between water particle A and particle B. The specific formula is x AB =x A-x B ; The water particle A supports the domain B. s The kernel length of size .

[0161] In this embodiment, the process of measuring the outflow rate of water body and the height of water level of dam overtopping and dam failure includes:

[0162]

[0163] Where C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflowing particle C.

[0164] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A numerical simulation prediction method for the breach development and outflow of an earth-rock dam considering overtopping dam breach, characterized in that: The following steps are involved: Step 1: constructing a three-dimensional geometric model of the earth-rock dam soil body and a three-dimensional geometric model of the water body, and assigning parameters to the three-dimensional geometric model of the earth-rock dam soil body and the three-dimensional geometric model of the water body; Step 2: Meshing the three-dimensional geometric model of the earth-rock dam soil and the three-dimensional geometric model of the water body, taking the grid center as the particle coordinate and the grid volume as the particle volume, to generate a soil particle discrete model and a water particle discrete model; Step 3: Based on the soil particle discrete model and the water particle discrete model, initialize the soil and water information, set the inlet and outlet boundary conditions, and calculate the water flow velocity and water level at the inflow boundary; wherein setting the inlet and outlet boundary conditions includes: inflow boundary conditions and outflow boundary conditions; Step 4: Based on the water velocity and water level of the inflow boundary, a soil erosion-deformation model and a coupled model of soil-water interaction are established. Based on the coupled model, the motion equations of soil particles and the motion equations of water particles are calculated by the SPH method; Step 5: Based on the boundary conditions, the water level, the motion equations of the soil particles, the motion equations of the water particles, and the coupled model of soil-water interaction, measure the breach erosion width and soil erosion depth of the dam after the earth slope of the earth-rock dam collapses; Step six: Based on the soil erosion-deformation model, the coupled model of soil-water interaction, the motion equations of the soil particles, the motion equations of the water particles and the outflow boundary conditions, measure the outflow flow and water level of the water body causing the overtopping dam break.

2. The numerical simulation prediction method for earth-rock dam breach development and outflow considering overtopping dam breach according to claim 1 is characterized in that: The methods for setting inlet and outlet boundary conditions and calculating the water flow velocity and water level at the inflow boundary include: The fluid velocity on the inflow boundary is obtained from the one-dimensional momentum conservation equation, expressed as: In the formula, Expressed as a time differential; It is expressed as the partial derivative in the x direction; h is the water level; g is the gravitational acceleration; v is the horizontal velocity; the equation is discretized using the forward difference method as follows: Among them l in is the length of the inflow area, h in is the initial water level of the reservoir in the upstream area; Δt is the time step; g is the gravitational acceleration; h se and v se is the actual water level and average flow rate of the reservoir in the upstream area; v in represents the average fluid velocity on the inflow boundary of the upstream region; The water level at the inflow boundary is calculated as: in, Expressed as time differential; h r Indicates the water level in the reservoir; A r =dV r / dh r Represents the given storage capacity curve V r =h(h r ) or empirical formula at water level h r The reservoir area at that time; Q in is the inflow obtained based on given hydrological data; Q out It is the outflow flow including soil and water breach.

3. The numerical simulation prediction method for earth-rock dam breach development and outflow considering overtopping dam breach according to claim 1 is characterized in that: Based on the coupling model, the method of calculating the motion equations of soil particles and the motion equations of water particles by the SPH method includes: The discrete mass conservation and momentum conservation equations of soil particle a are: Where a and b represent soil particles a and b respectively; m a , ρ a are the mass, acceleration, velocity and mass density associated with particle a respectively; represents the differential of the mass density of particle a with respect to time; represents the sum of all neighboring particles in the support domain of particle a; σ ab is the stress tensor at soil particle a including the artificial viscosity tensor; σ ba is the stress tensor at soil particle b including the artificial viscosity tensor; g ab It represents the discrete form of the derivative of the auxiliary kernel function of soil particle a to soil particle b. The specific formula is g ab =g(x a -x b ,l s ); Similarly, g ba It represents the discrete form of the derivative of the auxiliary kernel function of soil particle b to soil particle a. The specific formula is g ba =g(x b -x a ,l s ); g represents the auxiliary kernel function gradient; x a represents the coordinates of particle a; x b represents the coordinates of particle b; l s is the soil particle a supporting domain B m The kernel length of the size; V a and V b are the volumes of soil particles a and b respectively; δ ρ is a constant, generally taken as 0.1; c0 is the soil sound velocity; ψ ab represents the density difference between soil particles a and b; similarly, ψ ba represents the density difference between soil particle b and soil particle a; where x ab is the coordinate difference between soil particles a and b, and the specific formula is x ab =x a -x b ; Similarly where x ba is the coordinate difference between soil particle b and particle a. The specific formula is x ba =x b -x a ; The supporting domain B is soil particle a m The kernel length of size ; The discrete mass conservation and momentum conservation equations of water particle A are: Where A and B represent water particles A and B respectively; m A , ρ A are the mass, acceleration, velocity, and mass density associated with particle A; represents the differential of the mass density of particle A with respect to time; represents the sum of all neighboring particles in the support domain of particle A; σ AB is the stress tensor at water particle A including the artificial viscosity tensor; σ BA is the stress tensor at water particle B including the artificial viscosity tensor; g AB It represents the discrete form of the derivative of the auxiliary kernel function of water particle A to water particle B. The specific formula is g AB =g(x A -x B ,l s ); Similarly, g BA It represents the discrete form of the derivative of the auxiliary kernel function of water particle B to water particle A. The specific formula is g BA =g(x B -x A ,l s ); g represents the auxiliary kernel function gradient; x A represents the coordinates of particle A; x B represents the coordinates of particle B; l s The water particle A supports domain B. s The kernel length of the size; V A and V B are the volumes of water particles A and B respectively; δ ρ is a constant, usually taken as 0.1; c0 is the speed of sound in water; ψ AB represents the density difference between water particle A and water particle B; similarly, ψ BA It represents the density difference between water particle B and water particle A; where x AB is the coordinate difference between water particle A and water particle B. The specific formula is x AB =x A -x B ; Similarly where x BA is the coordinate difference between water particle B and water particle A. The specific formula is x BA =x B -x A ; The water particle A supports the domain B. s The kernel length of size .

4. The numerical simulation prediction method for earth-rock dam breach development and outflow considering overtopping dam breach according to claim 1 is characterized in that: Methods for measuring the outflow and water level of a dam overtopping event include: Where C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflowing particle C.

5. A numerical simulation prediction system for earth-rock dam breach development and outflow considering overtopping dam breach, characterized in that: include: A division module, a discrete module, a first calculation module, a second calculation module, a first measurement module, and a second measurement module; The partitioning module is used to construct a three-dimensional geometric model of the earth body and a three-dimensional geometric model of the water body of the earth-rock dam, and to assign parameters to the three-dimensional geometric model of the earth body and the three-dimensional geometric model of the water body; The discrete module is used to mesh the three-dimensional geometric model of the earth-rock dam soil and the three-dimensional geometric model of the water body, taking the grid center as the particle coordinate and the grid volume as the particle volume, to generate a soil particle discrete model and a water particle discrete model; The first calculation module is used to initialize the information of the soil and water bodies based on the soil particle discrete model and the water particle discrete model, set the inlet and outlet boundary conditions, and calculate the water flow velocity and water level height at the inflow boundary; wherein setting the inlet and outlet boundary conditions includes: inflow boundary conditions and outflow boundary conditions; The second calculation module is used to establish a soil erosion-deformation model and a coupled model of soil-water interaction based on the water flow velocity and water level of the inflow boundary, and calculate the motion equations of soil particles and the motion equations of water particles based on the coupled model by the SPH method; The first measurement module is used to measure the breach erosion width and soil erosion depth of the earth-rock dam after the earth slope collapses based on the boundary conditions, the water level, the motion equations of the soil particles, the motion equations of the water particles, and the coupled model of soil-water interaction; The second measurement module is used to measure the water outflow and water level height of the overtopping dam burst based on the soil erosion-deformation model, the soil-water interaction coupling model, the motion equations of the soil particles, the motion equations of the water particles and the outflow boundary conditions.

6. The numerical simulation prediction system for earth-rock dam breach development and outflow considering overtopping and dam breach according to claim 5 is characterized in that: The process of setting the inlet and outlet boundary conditions and calculating the water flow velocity and water level at the inflow boundary includes: The fluid velocity on the inflow boundary is obtained from the one-dimensional momentum conservation equation, expressed as: In the formula, Expressed as a time differential; It is expressed as the partial derivative in the x direction; h is the water level; g is the gravitational acceleration; v is the horizontal velocity; the equation is discretized using the forward difference method as follows: Among them l in is the length of the inflow area, h in is the initial water level of the reservoir in the upstream area; Δt is the time step; g is the gravitational acceleration; h se and v se is the actual water level and average flow rate of the reservoir in the upstream area; v in represents the average fluid velocity on the inflow boundary of the upstream region; The water level at the inflow boundary is calculated as: in, Expressed as time differential; h r Indicates the water level in the reservoir; A r =dV r / dh r Represents the given storage capacity curve V r =h(h r ) or empirical formula at water level h r The reservoir area at that time; Q in is the inflow obtained based on given hydrological data; Q out It is the outflow flow including soil and water breach.

7. The numerical simulation prediction system for earth-rock dam breach development and outflow considering overtopping and dam breach according to claim 5 is characterized in that: Based on the coupling model, the process of calculating the motion equations of soil particles and the motion equations of water particles by the SPH method includes: The discrete mass conservation and momentum conservation equations of soil particle a are: Where a and b represent soil particles a and b respectively; m a , ρ a are the mass, acceleration, velocity and mass density associated with particle a respectively; represents the differential of the mass density of particle a with respect to time; represents the sum of all neighboring particles in the support domain of particle a; σ ab is the stress tensor at soil particle a including the artificial viscosity tensor; σ ba is the stress tensor at soil particle b including the artificial viscosity tensor; g ab It represents the discrete form of the derivative of the auxiliary kernel function of soil particle a to soil particle b. The specific formula is g ab =g(x a -x b ,l s ); Similarly, g ba It represents the discrete form of the derivative of the auxiliary kernel function of soil particle b to soil particle a. The specific formula is g ba =g(x b -x a ,l s ); g represents the auxiliary kernel function gradient; x a represents the coordinates of particle a; x b represents the coordinates of particle b; l s is the soil particle a supporting domain B m The kernel length of the size; V a and V b are the volumes of soil particles a and b respectively; δ ρ is a constant, generally taken as 0.1; c0 is the soil sound velocity; ψ ab represents the density difference between soil particles a and b; similarly, ψ ba represents the density difference between soil particle b and soil particle a; where x ab is the coordinate difference between soil particles a and b, and the specific formula is x ab =x a -x b ; Similarly where x ba is the coordinate difference between soil particle b and particle a. The specific formula is x ba =x b -x a ; The supporting domain B is soil particle a m The kernel length of size ; The discrete mass conservation and momentum conservation equations of water particle A are: Where A and B represent water particles A and B respectively; m A , ρ A are the mass, acceleration, velocity, and mass density associated with particle A; represents the differential of the mass density of particle A with respect to time; represents the sum of all neighboring particles in the support domain of particle A; σ AB is the stress tensor at water particle A including the artificial viscosity tensor; σ BA is the stress tensor at water particle B including the artificial viscosity tensor; g AB It represents the discrete form of the derivative of the auxiliary kernel function of water particle A to water particle B. The specific formula is g AB =g(x A -x B ,l s ); Similarly, g BA It represents the discrete form of the derivative of the auxiliary kernel function of water particle B to water particle A. The specific formula is g BA =g(x B -x A ,l s ); g represents the auxiliary kernel function gradient; x A represents the coordinates of particle A; x B represents the coordinates of particle B; l s The water particle A supports domain B. s The kernel length of the size; V A and V B are the volumes of water particles A and B respectively; δ ρ is a constant, usually taken as 0.1; c0 is the speed of sound in water; ψ AB represents the density difference between water particle A and water particle B; similarly, ψ BA It represents the density difference between water particle B and water particle A; where x AB is the coordinate difference between water particle A and water particle B. The specific formula is x AB =x A -x B ; Similarly where x BA is the coordinate difference between water particle B and water particle A. The specific formula is x BA =x B -x A ; The water particle A supports the domain B. s The kernel length of size .

8. The numerical simulation prediction system for earth-rock dam breach development and outflow considering overtopping and dam breach according to claim 5 is characterized in that: The process of measuring the outflow and water level of the overtopping dam includes: Where C represents soil or water particles; nt is the number of particles flowing out of the buffer zone; l out is the width of the outflow buffer; v x,c is the x-direction velocity of the outflowing particle C; V c is the volume of the outflowing particle C.

Images