Multi-block rock mass landslide surge risk assessment method and system

Abstract

The invention discloses a multi-block rock mass landslide surge risk assessment method and system, and the method comprises the following steps: 1, building a landslide and water body model, and carrying out the assignment of material parameters; 2, determining boundary conditions of the model, and performing discrete processing on the model according to precision requirements; 3, a block discrete element method is adopted, the contact acting force between multiple block rock masses is calculated, and the resultant acceleration of all the block masses is updated in real time; 4, simulating a water body and surge evolution process by adopting a smoothed particle dynamics method; 5, performing coupling calculation between the water body and the rock mass to obtain a coupling acceleration; and step 6, adding the calculated coupling acceleration to the step 3 and the step 4, carrying out speed and coordinate displacement updating on the whole device, and carrying out cyclic calculation until a calculation time step requirement is met. The device can accurately simulate the pressure distribution of the coupling interface, and can accurately predict the collision movement process of the rock mass and the generation and propagation process of the surge.

Description

Technical Field

[0001] The present invention belongs to the technical field of landslide-induced surge disaster simulation, and particularly relates to a method and device for risk assessment of multi-block rock mass landslide-induced surge. Background Art

[0002] In the field of geotechnical engineering disasters, landslide-induced surge is a very important research topic, and its disaster influence range even far exceeds the landslide body itself. Compared with soil landslides, rock mass landslides often show more obvious overall movement characteristics and contour characteristics during the process of instability, movement, and entry into water. When there are fewer joints inside the landslide body or no through areas are formed during its movement, in numerical simulation studies, the whole of the rock mass landslide-induced surge is usually regarded as a rigid body for equivalent treatment, that is, the local collision and fragmentation phenomenon between it and the bedrock during the movement process is ignored, which brings convenience to the selection of numerical simulation research tools. However, in actual engineering, the bedrock near the ground surface with strong weathering contains internal structural fissures due to long-term geological tectonic movements, and the core integrity is poor. The joints and fissures existing in the landslide body (taking the bedrock as an example) will cause different degrees of fragmentation during the landslide process, thus having different effects on the surge wave amplitude. Therefore, how to accurately simulate the whole process of multi-block landslide-induced surge is an urgent problem to be solved at present.

[0003] Due to its mesh-free dependence property, the Smoothed Particle Hydrodynamics (SPH) method has gradually been applied to simulate the propagation process of water surges in recent years, and shows advantages in dealing with complex free surface conditions such as wave impact and fragmentation. In the simulation of rigid body landslide-induced surge, existing research usually discretizes the rigid body boundary of the landslide body as the movement boundary of SPH fluid particles, and directly solves the coupling force between the rigid body and the water body through the force calculation of SPH boundary particles. This method does not require additional setting of boundary conditions during the calculation process, and also reduces the additional error introduced by the coefficient selection in the empirical formula method. However, it puts forward higher requirements for the simulation accuracy of boundary pressure fluctuations, and the interaction process between rock masses needs to be considered additionally in the overall calculation framework. Summary of the Invention

[0004] The present invention is proposed to solve the above deficiencies, and aims to provide a method and device for risk assessment of multi-block rock mass landslide-induced surge, which can consider the movement collision effect between multi-block rock masses and accurately simulate the induced propagation process of rock mass landslide-induced surge.

[0005] To achieve the above objectives, the present invention adopts the following solutions:

[0006] A method for risk assessment of multi-block rock mass landslide-induced surge, comprising the following steps:

[0007] Step 1: Determine the dimensions and sizes of the overall site of the landslide and surging wave to be simulated, establish the landslide and water body models, and assign parameters according to the material properties;

[0008] Step 2: Determine the boundary conditions of the model according to the boundary characteristics of the research site. Discretize the water body according to the analysis accuracy requirements, and perform discretization processing with the same accuracy on the boundary of the rock mass landslide body;

[0009] Step 3: The simulation calculation of multi-block rock masses adopts the block discrete element method. Select the discrete element contact model, calculate the contact force between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and update the resultant acceleration of each block in real time;

[0010] Step 4: The simulation calculation of the water body and surging wave adopts the smoothed particle hydrodynamics method. The fluid in the flow field is discretized into a series of fluid particles and calculated through the integral form of the Navier-Stokes equation;

[0011] Step 5: Conduct the coupling calculation between the water body and the rock mass to obtain the coupling acceleration;

[0012] Step 6: Add the calculated coupling acceleration to Steps 3 and 4, update the velocity and coordinate displacement of the overall device, and repeat the process of Steps 1 to 5 until the calculation time step requirements are met.

[0013] As a preferred implementation manner, in Step 2, the discretization process of the water body is simulated by using finite element mesh generation software. The mesh size is the target size of the particles, and the initialization of the water body discrete particles is carried out by exporting the node coordinates.

[0014] As a preferred implementation manner, Step 3 specifically includes the following steps:

[0015] Step 3.1: Discretize the boundary of the block discrete element into a series of particle discrete elements. Calculate the resultant force of the block discrete element by statistically analyzing the forces of the series of particle discrete elements. The particle size of the particle discrete element is significantly smaller than the block boundary size to accurately characterize the boundary contour of the block discrete element;

[0016] Step 3.2: The calculation of the boundary particle discrete element is carried out according to the force-displacement relationship, and its contact model is as follows:

[0017] F contact = F spring + F damp (1)

[0018] F n = k n δ n n + F ndamp (2)

[0019] Ft = max[(F t ) T-Δt + k t δ t + F tdamp , μ||F n ||] (3)

[0020] Where: F contact is the contact force, F spring is the contact elastic force, F damp is the contact damping force, n represents the common normal direction of the contact position, F n is the normal contact force, δ n is the normal overlap between particles, (F t ) T-Δt is the tangential contact force at this contact position at the end of the previous time step, δ t is the relative displacement between contact positions within a single time step, F ndamp and F tdamp represent the damping forces in the normal and tangential directions respectively, μ is the dynamic friction coefficient between blocks, taking the smaller value of the two dynamic friction coefficients, k n and k t are the normal stiffness and tangential stiffness of the contact interface respectively, and their calculation formulas are as follows:

[0021]

[0022] Where: G* = 0.5(G i + G j ), R* = 2R i R j / (R i + R j ), v * = 0.5(v i + v j ), i and j respectively represent the two block elements i and j to which the contact interface belongs; G is the shear modulus of the block element, R represents the radius of the sphere (disk) composed of the block element, v is the Poisson's ratio of the block element, s overlap represents the contact overlap distance. For the boundary, the radius of the sphere (disk) it composes is set to infinity, i.e., R* = 2R 块体 ;

[0023] Step 3.3: Calculate the resultant contact force and resultant moment:

[0024]

[0025] Where: F total，i is the resultant force on particle i, F contact,ij is the contact force between two particles i and j, M total,iis the resultant moment acting on particle i, N represents the total number of contact points with block i, represents the distance from contact point j to the centroid of block i, u i is the velocity vector of particle i, t is time, F gravity is the gravitational force, F coupling is the coupling force, and the calculation process of this force will be elaborated in step 5, I i is the moment of inertia of block i, M coupling is the coupling moment, I i is the moment of inertia of particle i.

[0026] As a preferred implementation manner, step 4 specifically includes the following steps:

[0027] Step 4.1: Convert the partial differential equation solution format of the fluid Navier-Stokes equation into a derivative integration format for the smooth kernel function. The smooth kernel function W ij = W(r i - r j , h), and its specific calculation format is:

[0028]

[0029] In the formula, q = ||r i - r j || / h, where r i and r j represent the displacements of particles i and j respectively, h is the kernel function smoothing length; α dim is a parameter related to the calculation dimension, taking 7 / (4πh 2 ) in two-dimensional conditions and 7 / (8πh 3 ) in three-dimensional conditions;

[0030] Step 4.2: Calculate the discrete fluid particle acceleration. The specific calculation format of the Navier-Stokes equation in integral format is as follows:

[0031]

[0032] In the formula: ρ i and ρ j are the densities of particles i and j respectively, t is time, u j and u i are the velocity vectors of particles j and i respectively, is the displacement partial differential form of the kernel function, W ij represents the value of the kernel function at the relative position of particles i and j, V j is the volume of particle j, m j is the mass of particle j, p i and p jThe pressures of particles i and j respectively, c 0 is the artificial sound speed. To ensure the weak compressibility assumption of the fluid, the density change gradient needs to be less than 1%. Therefore, c 0 should be taken not less than 10 times the maximum velocity of the flow field. ρ 0 is the reference density (at 20°C and one standard atmosphere, the density of water is taken as 1000 kg / m 3 ), γ = 7 is a constant, F vis is the viscous force, g is the acceleration due to gravity, with a value of -9.81 m / s 2 ;

[0033] Step 4.3, the viscous force F vis is calculated as shown in the following formula, and an artificial viscosity term is additionally added during the free surface calculation to avoid uneven particle distribution caused by surface tension instability:

[0034]

[0035] In the formula: Dim is the dimension (2 for two - dimensional and 3 for three - dimensional), μ is the fluid viscosity coefficient, ρ i and ρ j are the densities of particles i and j respectively, N is the set of acting particles, u i and u j are the velocity vectors of particles i and j respectively, r i and r j are the displacement vectors of particles i and j respectively, h i is the smoothing length of particle i, is the displacement partial differential form of the kernel function, W ij represents the value of the kernel function at the relative position of particles i and j, V j is the volume of particle j,, F vis,art is the artificial viscous force, α is a constant with a value of 0.02, c 0 is the artificial sound speed.

[0036] As a preferred implementation manner, in step 5, the coupled acceleration is calculated based on the following formula:

[0037]

[0038] In the formula: F Di represents the force exerted by the virtual particle i on the rock mass element D, that is, the resultant force of the fluid acting on a single virtual particle; m i and m j represent the masses of particles i and j respectively; p i and p j represent the pressures of particles i and j respectively; ρ i and ρ jrepresent the densities of particles i and j respectively; W is the smoothing kernel function, is the partial differential form of the kernel function with respect to the displacement between particles i and j; Dim is the computational dimension; the pressure p of the particle is calculated in step 4; j is the sum of all fluid particles within the computational domain of the boundary particle i; μ is the fluid viscosity coefficient; u i and u j represent the velocity vectors of particles i and j respectively; r i and r j represent the displacement vectors of particles i and j respectively; h i is the smoothing length of particle i, which is related to the form of the kernel function. In two-dimensional conditions, it is usually taken as 1.5 - 2 times the initial particle spacing, and in three-dimensional conditions, it is taken as 1.5 times the particle spacing; V is the particle volume; F coupling represents the coupling force; M coupling represents the coupling moment; Γ represents the set of all virtual particles acting on the same block element; r Di is the distance from the boundary particle i to the center of gravity of the block D.

[0039] As a preferred embodiment, step 6 specifically includes the following steps:

[0040] Step 6.1, perform the prediction step and correction step updates within a single time step:

[0041]

[0042] In the formula: is the velocity vector of particle i after the (n + 1 / 2) step, is the velocity vector of particle i after the n step, Δt is the time step, is the acceleration vector of particle i after the n step, is the angular velocity vector of particle i after the (n + 1 / 2) step, is the angular velocity vector of particle i after the n step, is the angular acceleration vector of particle i after the n step, is the rotation angle of particle i after the (n + 1 / 2) step, is the rotation angle of particle i after the n step, is the density of particle i after the (n + 1 / 2) step, is the density of particle i after the n step, is the density change gradient of particle i after the n step, p is the pressure;

[0043] Step 6.2, correction step:

[0044]

[0045] In the formula: is the velocity vector of particle i after the (n + 1 / 2)-th step, is the velocity vector of particle i after the n-th step, and Δt is the time step, is the acceleration vector of particle i after the (n + 1 / 2)-th step, is the angular velocity vector of particle i after the (n + 1 / 2)-th step, is the angular velocity vector of particle i after the n-th step, is the angular acceleration vector of particle i after the (n + 1 / 2)-th step, is the rotation angle of particle i after the (n + 1 / 2)-th step, is the rotation angle of particle i after the n-th step, is the density of particle i after the (n + 1 / 2)-th step, is the density of particle i after the n-th step, is the density change gradient of particle i after the (n + 1 / 2)-th step, is the displacement vector of particle i after the (n + 1 / 2)-th step, is the displacement vector of particle i after the n-th step;

[0046] Step 6.3, Parameter Update:

[0047]

[0048] In the formula: is the velocity vector of particle i after the (n + 1)-th step, is the velocity vector of particle i after the (n + 1 / 2)-th step, is the velocity vector of particle i after the n-th step, is the angular velocity vector of particle i after the (n + 1)-th step, is the angular velocity vector of particle i after the (n + 1 / 2)-th step, is the angular velocity vector of particle i after the n-th step, is the rotation angle of particle i after the (n + 1)-th step, is the rotation angle of particle i after the (n + 1 / 2)-th step, is the rotation angle of particle i after the n-th step, is the density of particle i after the (n + 1)-th step, is the density of particle i after the (n + 1 / 2)-th step, is the density of particle i after the n-th step, is the displacement vector of particle i after the (n + 1)-th step, is the displacement vector of particle i after the (n + 1 / 2)-th step, is the displacement vector of particle i after the n-th step;

[0049] In the DEM-SPH framework, the time steps of the two modules are kept consistent and determined according to the following formula (22):

[0050]

[0051] where: Δt is the single-step duration, m i is the mass of the discrete element i of the block, k ni is the maximum contact stiffness of the discrete element i, h i is the smoothing length of the SPH particle i, v is the kinematic viscosity coefficient of the fluid, a i is the acceleration of particle i, c 0 is the artificial sound speed, u max is the maximum flow velocity of the flow field.

[0052] The present invention also provides a multi-block rock mass landslide surge risk assessment system, which is characterized by including:

[0053] A modeling module that determines the dimensions and sizes of the overall site of the landslide and surge to be simulated, establishes a landslide and water body model, and assigns parameters according to the material properties;

[0054] A preprocessing module that determines the boundary conditions of the model according to the boundary characteristics of the research site, discretizes the water body according to the analysis accuracy requirements, and performs discretization processing on the boundary of the rock mass landslide body with the same accuracy;

[0055] A rock mass landslide simulation module that uses the block discrete element method for the simulation calculation of multi-block rock masses, selects a discrete element contact model, calculates the contact forces between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and updates the resultant acceleration of each block in real time;

[0056] A surge simulation module that uses the smoothed particle hydrodynamics method for the simulation calculation of water bodies and surges, discretizes the flow field water body into a series of fluid particles, and calculates through the Navier-Stokes equation in integral form;

[0057] A time step update module that adds the calculated coupled acceleration to the rock mass landslide simulation module and the surge simulation module, updates the velocity and coordinate displacement of the overall device, and repeats the above module update process until the calculation time step requirements are met;

[0058] A control module that is communicatively connected to the modeling module, the preprocessing module, the rock mass landslide simulation module, the surge simulation module, and the time step update module, and controls their operations.

[0059] As a preferred implementation manner, the multi-block rock mass landslide surge risk assessment system further includes:

[0060] The regional division module divides the research area of the actual project, generates a numerical model in STL format, imports it into the pre-processing software for discretization processing, and generates node coordinates for export to subsequent modules.

[0061] As a preferred embodiment, the multi-block rock mass landslide and surge risk assessment system further includes:

[0062] The post-processing module is used to output the calculation results as data files, directly process them using post-processing software, and obtain visual calculation results.

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

[0064] The multi-block rock mass landslide and surge risk assessment device provided by the present invention solves the problem that it is difficult to simultaneously simulate the collision between multi-block rock masses and the fluid-structure coupling effect with water bodies. It can accurately simulate the pressure distribution at the coupling interface, accurately predict the rock mass collision movement process, the generation and propagation process of surges, and can also accurately invert the local impact flow pattern. The provided multi-block rock mass landslide and surge risk assessment device can simulate various working conditions such as single-block and multi-block rock mass landslides and surges, providing a reliable numerical calculation method for rock mass landslide and surge disaster assessment and prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 is a flowchart of the multi-block rock mass landslide and surge risk assessment method according to the present invention;

[0066] Figure 2 is a schematic diagram of a landslide body, water body and terrain according to an embodiment of the present invention;

[0067] Figure 3 is a schematic diagram of a particle retrieval method according to an embodiment of the present invention;

[0068] Figure 4 is a flowchart of particle retrieval according to an embodiment of the present invention;

[0069] Figure 5 is a diagram of the local flow pattern simulation results according to an embodiment of the present invention, where (a) is the result diagram of the physical model test, and (b) is the numerical simulation result diagram calculated using the present invention;

[0070] Figure 6 is a diagram of the change results of the surge height at the wave gauge in three groups of implementation working conditions according to an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0071] The following will describe in detail the specific implementation schemes of the multi-block rock mass landslide and surge risk assessment method and device according to the present invention with reference to the accompanying drawings.

[0072] A method for risk assessment of multi - block rock mass landslide surges according to the present invention includes the following steps:

[0073] Step 1: Determine the dimensions and size of the overall site of the landslide and surge to be simulated, establish a landslide and water body model, and assign parameters according to material properties;

[0074] Step 2: Determine the boundary conditions of the model according to the boundary characteristics of the research site. According to the analysis accuracy requirements, discretize the water body and perform the same - precision discretization treatment on the boundary of the rock mass landslide body. The discretization process of the water body is simulated using finite - element mesh generation software. The mesh size is the particle target size, and the initialization of the water - body discrete particles is carried out by exporting the node coordinates.

[0075] Step 3: The simulation calculation of multi - block rock mass adopts the block - discrete element method. Select a discrete - element contact model, calculate the contact force between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and update the resultant acceleration of each block in real time. It specifically includes the following steps:

[0076] Step 3.1: Discretize the block - discrete element boundary into a series of particle - discrete elements. Calculate the resultant force of the block - discrete element by statistically analyzing the forces of the series of particle - discrete elements. The particle size of the particle - discrete element is significantly smaller than the block - boundary size to accurately represent the block - discrete element boundary profile;

[0077] Step 3.2: The calculation of the boundary particle - discrete element is carried out according to the force - displacement relationship, and its contact model is as follows:

[0078] F contact =F spring +F damp (1)

[0079] F n =k n δ n n+F ndamp (2)

[0080] F t =max[(F t ) T-Δt +k t δ t +F tdamp ,μ||F n ||] (3)

[0081] In the formula: F contact is the contact force, F spring is the contact elastic force, F damp is the contact damping force, n represents the common - normal direction of the contact position, F n is the normal contact force, δn is the normal overlap between particles, (F t ) T-Δt is the tangential contact force at this contact position at the end of the previous time step, δ t is the relative displacement between contact positions within a single time step, F ndamp and F tdamp represent the normal and tangential damping forces respectively, μ is the dynamic friction coefficient between the blocks, taking the smaller value of the two dynamic friction coefficients, k n and k t are the normal stiffness and tangential stiffness of the contact interface respectively, and their calculation formulas are as follows:

[0082]

[0083] In the formula: G* = 0.5(G i + G j ), R* = 2R i R j / (R i + R j ), v * = 0.5(v i + v j ), i and j respectively represent the two block elements i and j to which the contact interface belongs; G is the shear modulus of the block element, R represents the radius of the sphere (disc) composed of the block element, v is the Poisson's ratio of the block element, s overlap represents the contact coincidence distance. For the boundary, the radius of the sphere (disc) it composes is set to infinity, that is, R* = 2R 块体 ;

[0084] Step 3.3: Calculate the resultant contact force and resultant moment:

[0085]

[0086] In the formula: F total，i is the resultant force on particle i, F contact,ij is the contact force between the two particles i and j, M total,i is the resultant moment on particle i, N represents the total number of contact points with block i, r Gij represents the distance from contact point j to the center of gravity of block i, u i is the velocity vector of particle i, t is time, F gravity is the gravitational force, F coupling is the coupling force, and the calculation process of this force will be elaborated in Step 5, I i is the moment of inertia of block i, M coupling is the coupling moment, I i is the moment of inertia of particle i.

[0087] Step 4. The simulation calculation of water body and surging waves adopts the smoothed particle hydrodynamics method. The water body in the flow field is discretized into a series of fluid particles and calculated through the integral form of the Navier-Stokes equation, which specifically includes the following steps:

[0088] Step 4.1. The partial differential equation solution format of the fluid Navier-Stokes equation is transformed into a derivative integral format of the smoothing kernel function. The smoothing kernel function W ij = W(r i - r j , h), and its specific calculation format is:

[0089]

[0090] In the formula, q = ||r i - r j || / h, where r i and r j represent the displacements of particles i and j respectively, h is the smoothing length of the kernel function; α dim is a parameter related to the calculation dimension, taking 7 / (4πh 2 ) under two-dimensional conditions and 7 / (8πh 3 ) under three-dimensional conditions;

[0091] Step 4.2. Calculate the acceleration of discrete fluid particles. The specific calculation format of the integral form of the Navier-Stokes equation is as follows:

[0092]

[0093] In the formula: ρ i and ρ j are the densities of particles i and j respectively, t is the time, u j and u i are the velocity vectors of particles j and i respectively, is the displacement partial differential form of the kernel function, W ij represents the value of the kernel function at the relative position of particles i and j, V j is the volume of particle j, m j is the mass of particle j, p i and p j are the pressures of particles i and j respectively, c 0 is the artificial sound speed. To ensure the weak compressibility assumption of the fluid, the density change gradient needs to be less than 1%. Therefore, the value of c 0 should not be less than 10 times the maximum velocity of the flow field. ρ 0 is the reference density (at 20°C and one standard atmosphere, the density of water is taken as 1000 kg / m 3 ), γ = 7 is a constant, F visis the viscous force, and g is the acceleration due to gravity, with a value of -9.81 m / s 2 ;

[0094] Step 4.3, Viscous force F vis The calculation formula is as shown below, and an artificial viscosity term is additionally added during the free surface calculation to avoid uneven particle distribution caused by surface tension instability:

[0095]

[0096] In the formula: Dim is the dimension (2 for two-dimensional and 3 for three-dimensional), μ is the fluid viscosity coefficient, ρ i and ρ j are the densities of particles i and j respectively, N is the set of acting particles, u i and u j are the velocity vectors of particles i and j respectively, r i and r j are the displacement vectors of particles i and j respectively, h i is the smoothing length of particle i, is the displacement partial differential form of the kernel function, W ij represents the value of the kernel function at the relative position of particles i and j, V j is the volume of particle j,, F vis,art is the artificial viscous force, α is a constant with a value of 0.02, c 0 is the artificial sound speed.

[0097] Step 5. Perform the coupled calculation between the water body and the rock mass. The discrete particles at the boundary of the rock mass landslide body participate in the fluid calculation process in Step 4, and the coupled acceleration is calculated based on the following formula:;

[0098]

[0099] In the formula: F Di represents the force exerted by the virtual particle i on the rock mass element D, that is, the resultant force of the nearby fluid acting on a single virtual particle; m i and m j represent the masses of particles i and j respectively; p i and p j represent the pressures of particles i and j respectively; ρ i and ρ j represent the densities of particles i and j respectively; W is the smoothing kernel function, is the displacement partial differential form of the kernel function between particles i and j; Dim is the calculation dimension; the pressure p of the particle is calculated in Step 4; j is the sum of all fluid particles within the calculation domain of the boundary particle i; μ is the fluid viscosity coefficient; u i and u j represent the velocity vectors of particles i and j respectively; ri and r j represent the displacement vectors of particles i and j respectively; h i is the smoothing length of particle i, which is related to the form of the kernel function. In two-dimensional conditions, it is usually taken as 1.5 - 2 times the initial particle spacing, and in three-dimensional conditions, it is taken as 1.5 times the particle spacing; V is the particle volume; F coupling represents the coupling force; M coupling represents the coupling moment; Γ represents the set of all virtual particles acting on the same block element; r Di is the distance from the boundary particle i to the center of gravity of the block D.

[0100] Step 6: Add the calculated coupling acceleration to Steps 3 and 4, update the velocity and coordinate displacement of the overall device, and repeat the process of Steps 1 - 5 until the calculation time step requirement is met, which specifically includes the following steps:

[0101] Step 6.1, perform the prediction step and correction step updates within a single time step:

[0102]

[0103] In the formula: is the velocity vector of particle i after the (n + 1 / 2) step, is the velocity vector of particle i after the n step, Δt is the time step size, is the acceleration vector of particle i after the n step, is the angular velocity vector of particle i after the (n + 1 / 2) step, is the angular velocity vector of particle i after the n step, is the angular acceleration vector of particle i after the n step, is the rotation angle of particle i after the (n + 1 / 2) step, is the rotation angle of particle i after the n step, is the density of particle i after the (n + 1 / 2) step, is the density of particle i after the n step, is the density change gradient of particle i after the n step, p is the pressure;

[0104] Step 6.2, correction step:

[0105]

[0106] In the formula: is the velocity vector of particle i after the (n + 1 / 2) step, is the velocity vector of particle i after the n step, Δt is the time step size, is the acceleration vector of particle i after the (n + 1 / 2) step, is the angular velocity vector of particle i after the (n + 1 / 2)-th step, is the angular velocity vector of particle i after the n-th step, is the angular acceleration vector of particle i after the (n + 1 / 2)-th step, is the rotation angle of particle i after the (n + 1 / 2)-th step, is the rotation angle of particle i after the n-th step, is the density of particle i after the (n + 1 / 2)-th step, is the density of particle i after the n-th step, is the density change gradient of particle i after the (n + 1 / 2)-th step, is the displacement vector of particle i after the (n + 1 / 2)-th step, is the displacement vector of particle i after the n-th step;

[0107] Step 6.3, parameter update:

[0108]

[0109] p n+1 = f(ρ n+1 )(21)

[0110] where: is the velocity vector of particle i after the (n + 1)-th step, is the velocity vector of particle i after the (n + 1 / 2)-th step, is the velocity vector of particle i after the n-th step, is the angular velocity vector of particle i after the (n + 1)-th step, is the angular velocity vector of particle i after the (n + 1 / 2)-th step, is the angular velocity vector of particle i after the n-th step, is the rotation angle of particle i after the (n + 1)-th step, is the rotation angle of particle i after the (n + 1 / 2)-th step, is the rotation angle of particle i after the n-th step, is the density of particle i after the (n + 1)-th step, is the density of particle i after the (n + 1 / 2)-th step, is the density of particle i after the n-th step, is the displacement vector of particle i after the (n + 1)-th step, is the displacement vector of particle i after the (n + 1 / 2)-th step, is the displacement vector of particle i after the n-th step;

[0111] In the DEM-SPH framework, the time steps of the two modules are kept consistent and determined according to the following formula (22):

[0112]

[0113] where: Δt is the single-step duration, m i is the mass of the block discrete element i, k ni is the maximum contact stiffness of the discrete element i, h i is the smoothing length of the SPH particle i, v is the kinematic viscosity coefficient of the fluid, a i is the acceleration of the particle i, c 0 is the artificial sound speed, u max is the maximum flow velocity of the flow field.

[0114] The present invention also provides a multi-block rock mass landslide surge risk assessment system, which is characterized by comprising:

[0115] A modeling module, which determines the dimensions and sizes of the overall site of the landslide and surge to be simulated, establishes a landslide and water body model, and assigns parameters according to the material properties;

[0116] A preprocessing module, which determines the boundary conditions of the model according to the boundary characteristics of the research site, discretizes the water body according to the analysis accuracy requirements, and performs the same-precision discretization processing on the boundary of the rock mass landslide body;

[0117] A rock mass landslide simulation module, which uses the block discrete element method for the simulation calculation of multi-block rock masses, selects a discrete element contact model, calculates the contact forces between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and updates the resultant acceleration of each block in real time;

[0118] A surge simulation module, which uses the smoothed particle hydrodynamics method for the simulation calculation of water bodies and surges, discretizes the flow field water body into a series of fluid particles, and calculates through the Navier-Stokes equation in integral form;

[0119] A time step update module, which adds the calculated coupled acceleration to the rock mass landslide simulation module and the surge simulation module, updates the velocity and coordinate displacement of the overall device, and repeats the above module update process until the calculation time step requirements are met;

[0120] A control module, which is communicatively connected to the modeling module, the preprocessing module, the rock mass landslide simulation module, the surge simulation module, and the time step update module, and controls their operations;

[0121] A regional division module, which divides the research area of the actual project, generates a numerical model in STL format, imports it into the preprocessing software for discretization processing, and generates node coordinates for export to subsequent modules;

[0122] The post-processing module is used to output the calculation results in a data file and directly process them using post-processing software to obtain visual calculation results.

[0123] Embodiment:

[0124] As Figure 1 shown, the multi-block rock mass landslide surge risk assessment method and device adopted in this embodiment include the following implementation steps:

[0125] Step 1: Determine that the research scale is two-dimensional, import the model file, and perform visual processing of the pre-processing using open-source software such as Solidworks and Hypermesh to ensure the rationality of the calculation domain settings and determine the parameter values of the calculation model.

[0126] Step 2: As Figure 2 shown, determine the simulation conditions, select an appropriate target particle size for mesh generation, extract the mesh nodes and export them, name them Input.dat, generate discrete element particles and fluid particles, assign the attribute number 1 to the discrete element particles, assign the attribute number 2 to the fluid particles, and assign particle attribute parameters according to the landslide body parameters and water body parameters to complete the pre-processing stage. The particle parameter assignment is shown in Table 1.

[0127] Table 1 Parameter values for multi-block landslide surge calculation

[0128]

[0129] Traverse the overall particle group, retrieve the particles with attribute number 1 for block discrete element calculation, and retrieve the particles with attribute number 1 and substitute them into the smooth particle hydrodynamics equation for calculation. The parameter information of the two types of particles is stored separately.

[0130] Retrieve the interacting particle pairs. The retrieval method uses the grid retrieval method. The principle and process schematic diagram are as Figure 3 、 Figure 4 shown. Set a background grid in the calculation domain and number the grids where each particle is located. After processing in this way, under two-dimensional conditions, only need to retrieve whether there are particles in the 8 grids near the target grid, and then further retrieve whether there are interacting particles with the target unit in the nearby range. The interacting particles in the discrete element calculation refer to the contacting particles, and the interacting particles in the smooth particle hydrodynamics method calculation refer to the particles in the kernel function influence domain.

[0131] Step 3: Calculate the acting forces between the rock mass landslides according to the discrete element retrieved particle pairs. The calculation process refers to Formulas (1) to (9), and update the resultant force according to the force-displacement relationship.

[0132] Step 4: Calculate the water body force based on the fluid particle pairs. The calculation process refers to Formulas (10) to (13), and update the acceleration at each location in the flow field.

[0133] Step 5: Conduct the coupling calculation between the water body and the rock mass. The calculation process refers to Formulas (14) to (16), and add the forces and moments obtained from the coupling calculation to the rock mass and fluid calculation modules.

[0134] Step 6: Update the motion parameters such as the acceleration, velocity, and displacement of the overall calculation domain. During the above calculation process, the time step parameters of the landslide and surge calculation modules are updated simultaneously, and their integration formats are both the prediction-correction method.

[0135] After the above calculation is completed, export the parameters such as the coordinates, velocities, and pressures of the updated particles in the.dat file format, and then import them into the post-processing software for visualization processing, such as Figure 5 、 Figure 6 It can be seen that the accuracy of the proposed method and device can be verified by analyzing the flow state, landslide body trajectory, and surge height.

[0136] In summary, the present invention proposes a multi-block rock mass landslide and surge risk assessment method and device that can simultaneously consider the collision effects between multiple rock masses and the fluid-structure coupling effects with the water body, accurately simulate the pressure distribution at the coupling interface, accurately predict the rock mass collision movement process, the generation and propagation process of surges, and can also accurately invert the local impact flow state; moreover, the parameters selected in the present invention can be obtained through methods such as on-site material measurement and indoor tests, with high credibility and easy expansion, and can provide reliable numerical calculation means for the assessment and prediction of rock mass landslide and surge disasters.

[0137] In this embodiment, a multi-block rock mass landslide and surge risk assessment device that can automatically implement the above method of the present invention is also provided. The device includes a modeling module, a preprocessing module, a rock mass landslide simulation module, a surge simulation module, a landslide and surge coupling calculation module, a time step update module, and a control module.

[0138] The modeling module executes the content described in Step 1 above, determines the dimensions and sizes of the overall site of the landslide and surge to be simulated, establishes the landslide and water body models, and assigns parameters according to the material properties.

[0139] The preprocessing module executes the content described in Step 2 above, determines the boundary conditions of the model according to the boundary characteristics of the research site, discretizes the water body according to the analysis accuracy requirements, and performs the same-precision discretization processing on the boundary of the rock mass landslide body.

[0140] The rock mass landslide simulation module executes the content described in step 3 above. The simulation calculation of multiple-block rock masses adopts the block discrete element method, selects the discrete element contact model, calculates the contact force between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and updates the resultant acceleration of each block in real time.

[0141] The surge simulation module executes the content described in step 4 above. The simulation calculation of water bodies and surges adopts the smoothed particle hydrodynamics method. The flow field water body is discretized into a series of fluid particles and calculated through the integral form of the Navier-Stokes equation.

[0142] The landslide-surge coupling calculation module executes the content described in step 5 above, conducts the coupling calculation between the water body and the rock mass, and the discrete particles at the boundary of the rock mass landslide body participate in the fluid calculation process in step 4.

[0143] The time step update module executes the content described in step 6 above, adds the calculated coupling acceleration to the rock mass landslide simulation module and the surge simulation module, updates the velocity and coordinate displacement of the overall device, and repeats the above module update process until the requirements of the calculation time step are met.

[0144] The control module is communicatively connected to the modeling module, the preprocessing module, the rock mass landslide simulation module, the surge simulation module, and the time step update module, and controls their operations.

[0145] The regional division module divides the research area of the actual project, generates a numerical model in STL format, imports it into the preprocessing software for discretization processing, and generates node coordinates for export to subsequent modules.

[0146] The post-processing module is used to output the calculation results in a data file, directly process them using post-processing software, and obtain visual calculation results.

[0147] The above embodiments are merely illustrative examples of the technical solutions of the present invention. The multi-block rock mass landslide surge risk assessment method and device involved in the present invention are not limited solely to the content described in the above embodiments, but are subject to the scope defined by the claims. Any modification, supplement, or equivalent replacement made by those skilled in the art to the present invention based on this embodiment falls within the scope protected by the claims of the present invention.

Claims

1. A multi-block rock mass landslide surge risk assessment method, characterized in that: The following steps are involved: Step 1: Determine the dimensions and size of the overall site of the landslide and surge to be simulated, establish landslide and water body models, and assign parameters according to material properties; Step 2: According to the boundary characteristics of the research site, the boundary conditions of the model are determined. According to the analysis accuracy requirements, the water body is discretized, and the boundary of the rock mass and landslide is discretized with the same accuracy; Step 3: The simulation calculation of multi-block rock mass adopts the block discrete element method, selects the discrete element contact model, calculates the contact force between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and updates the resultant acceleration of each block in real time; Step 4: The simulation calculation of water body and surge wave adopts the smooth particle dynamics method. The flow field water body is discretized into a series of fluid particles and calculated by the integral format Navier-Stokes equation; Step 5: perform coupling calculation between the water body and the rock mass to obtain coupling acceleration; Step 6: Add the calculated coupled acceleration to steps 3 and 4, update the velocity and coordinate displacement of the entire device, and repeat steps 1 to 5 until the calculation time step requirements are met.

2. The multi-block rock mass landslide surge risk assessment method according to claim 1, characterized in that: In step 2, the water body discretization process is simulated by using finite element meshing software, the mesh size is the particle target size, and the discrete particles of the water body are initialized by deriving the node coordinates.

3. The multi-block rock mass landslide surge risk assessment method according to claim 1, characterized in that: The step 3 specifically comprises the following steps: Step 3.1, the boundary of the block discrete element is discretized into a series of particle discrete elements, and the resultant force of the block discrete element is calculated by calculating the forces of the series of particle discrete elements. The particle size of the particle discrete element is obviously smaller than the block boundary size, so as to accurately characterize the boundary contour of the block discrete element; Step 3.2: The calculation of the boundary particle discrete element is carried out according to the force-displacement relationship, and its contact model is as follows: F contact =F spring +F damp (1) F n =k n d n n+F ndamp (2) F t =max[(F t ) T-Δt +k t δ t +F tdamp ,μ||F n ||] (3) Where: F contact is the contact force, F spring is the contact elastic force, F damp is the contact damping force, n represents the direction of the common normal line of the contact position, F n is the normal contact force, δ n is the normal overlap between particles, (F t ) T-Δt is the tangential contact force at the contact position at the end of the previous time step, δ t is the relative displacement between contact positions in a single time step, F ndamp With F tdamp They represent the normal and tangential damping forces respectively, μ is the dynamic friction coefficient between the blocks, and the smaller value of the dynamic friction coefficients is taken, k n With k t are the normal stiffness and tangential stiffness of the contact interface, respectively, and their calculation formulas are: Where: G*=0.5(G i +G j ), R*=2R i R j / (R i +R j ), v * =0.5(v i +v j ), i and j represent the two block elements i and j to which the contact interface belongs respectively; G is the shear modulus of the block element, R represents the radius of the sphere composed of the block element, v is the Poisson's ratio of the block element, s overlap Represents the contact overlap distance. For the boundary, the radius of the sphere that makes up the boundary is set to infinity, that is, R*=2R 块体 ; Step 3.3: Calculate the contact force resultant force and moment: Where: F total，i is the resultant force on particle i, F contact,ij is the contact force between two particles i and j, M total,i is the total moment of particle i, N is the total number of contact points with block i, r Gij represents the distance from contact point j to the center of mass of block i, u i is the velocity vector of particle i, t is the time, F gravity is gravity, F coupling is the coupling force, I i is the moment of inertia of block i coupling is the coupling torque, I i is the moment of inertia of particle i.

4. The multi-block rock mass landslide surge risk assessment method according to claim 1, characterized in that: The step 4 specifically comprises the following steps: Step 4.1: The partial differential equation solution format of the fluid Navier-Stokes equation is converted into the derivative integral format of the smooth kernel function. The smooth kernel function W ij =W(r i -r j ,h), the specific calculation format is: In the formula, q = || r i -r j || / h, where r i and r j Respectively represent the displacement of particles i and j, h is the smooth length of the kernel function; α dim It is a parameter related to the calculation dimension. In two-dimensional conditions, it is 7 / (4πh 2 ), take 7 / (8πh for three-dimensional working conditions 3 ); Step 4.2: Calculate the acceleration of discrete fluid particles. The specific calculation format of the integral format Navier-Stokes equation is as follows: Where: i and ρ j are the densities of particles i and j respectively, t is the time, u j and u i are the velocity vectors of particles j and i, respectively, i W ij is the displacement partial differential form of the kernel function, W ij represents the value of the kernel function at the relative position of particles i and j, V j is the volume of particle j, m j is the mass of particle j, p i and p j are the pressures of particles i and j, c0 is the artificial sound speed, ρ0 is the reference density, γ=7 is a constant, F vis is the viscous force, g is the acceleration due to gravity, and its value is -9.81m / s 2 ; Step 4.3: Viscous force F vis The calculation format is shown below, and an artificial viscosity term is added during free surface calculation: Where: Dim is the dimension (2 for two-dimensional and 3 for three-dimensional), μ is the fluid viscosity coefficient, ρ i and ρ j are the densities of particles i and j respectively, N is the set of interacting particles, u i and u j are the velocity vectors of particles i and j, r i and r j are the displacement vectors of particles i and j, h i is the smooth length of particle i, ▽ i W ij is the displacement partial differential form of the kernel function, W ij represents the value of the kernel function at the relative position of particles i and j, V j is the volume of particle j, F vis,art is the artificial viscous force, α is a constant with a value of 0.02, and c0 is the artificial sound speed.

5. The multi-block rock mass landslide surge risk assessment method according to claim 1, characterized in that: In step 5, the coupled acceleration is calculated based on the following formula: Where: F Di represents the force exerted by virtual particle i on rock mass unit D, i.e., the resultant force exerted by the nearby fluid on a single virtual particle; m i and m j denote the masses of particles i and j respectively; p i and p j denote the pressure of particles i and j respectively; ρ i and ρ j denote the density of particles i and j respectively; W is the smooth kernel function, ▽ i W ij is the partial differential form of the kernel function for the displacement between particles i and j; Dim is the calculation dimension; p is the pressure of the particle; j is the sum of all fluid particles in the calculation domain of boundary particle i; μ is the fluid viscosity coefficient; u i and u j Represent the velocity vectors of particles i and j respectively; r i and r j denote the displacement vectors of particles i and j respectively; h i is the smooth length of particle i; V is the volume of particle; F coupling Represents coupling force; M coupling represents the coupling torque; Γ represents the set of all virtual particles acting on the same block unit; r Di is the distance from the boundary particle i to the center of gravity of the block D.

6. The multi-block rock mass landslide surge risk assessment method according to claim 1, characterized in that: The step 6 specifically comprises the following steps: Step 6.1, update the prediction step and correction step within a single time step: Where: is the velocity vector of particle i after the (n+1 / 2)th step, is the velocity vector of particle i after the (nth) step, Δt is the time step, is the acceleration vector of particle i after the (nth) step, is the angular velocity vector of particle i after the (n+1 / 2)th step, is the angular velocity vector of particle i after the (nth) step, is the angular acceleration vector of particle i after the (nth) step, is the turning angle of particle i after the (n+1 / 2) step, is the turning angle of particle i after the (nth) step, is the density of particle i after the (n+1 / 2) step, is the density of particle i after the (nth) step, is the density change gradient of particle i after the (nth) step, and p is the pressure; Step 6.2, correction step: Where: is the velocity vector of particle i after the (n+1 / 2)th step, is the velocity vector of particle i after the (nth) step, Δt is the time step, is the acceleration vector of particle i after the (n+1 / 2)th step, is the angular velocity vector of particle i after the (n+1 / 2)th step, is the angular velocity vector of particle i after the (nth) step, is the angular acceleration vector of particle i after the (n+1 / 2)th step, is the turning angle of particle i after the (n+1 / 2) step, is the turning angle of particle i after the (nth) step, is the density of particle i after the (n+1 / 2) step, is the density of particle i after the (nth) step, is the density change gradient of particle i after the (n+1 / 2) step, is the displacement vector of particle i after the (n+1 / 2)th step, is the displacement vector of particle i after the (nth) step; Step 6.3, parameter update: Where: is the velocity vector of particle i after the (n+1)th step, is the velocity vector of particle i after the (n+1 / 2)th step, is the velocity vector of particle i after the (nth) step, is the angular velocity vector of particle i after the (n+1)th step, is the angular velocity vector of particle i after the (n+1 / 2)th step, is the angular velocity vector of particle i after the (nth) step, is the turning angle of particle i after the (n+1)th step, is the turning angle of particle i after the (n+1 / 2) step, is the turning angle of particle i after the (nth) step, is the density of particle i after the (n+1)th step, is the density of particle i after the (n+1 / 2) step, is the density of particle i after the (nth) step, is the displacement vector of particle i after the (n+1)th step, is the displacement vector of particle i after the (n+1 / 2)th step, is the displacement vector of particle i after the (nth) step; In the DEM-SPH framework, the time steps of the two modules are kept consistent and determined according to the following formula (22): Where: Δt is the single step duration, m i is the mass of the block discrete element i, k ni is the maximum contact stiffness of discrete element i, h i is the smooth length of SPH particle i, v is the fluid kinematic viscosity coefficient, a i is the acceleration of particle i, c0 is the artificial sound speed, u max is the maximum flow velocity in the flow field.

7. A multi-block rock mass landslide surge risk assessment system, characterized in that: include: Modeling module, which determines the dimensions and size of the overall site of the landslide and surge to be simulated, establishes landslide and water body models, and assigns parameters according to material properties; The pre-processing module determines the boundary conditions of the model according to the boundary characteristics of the research site, discretizes the water body according to the analysis accuracy requirements, and discretizes the boundaries of the rock mass and landslide body with the same accuracy; Rock mass landslide simulation module: The simulation calculation of multi-block rock mass adopts the block discrete element method, selects the discrete element contact model, calculates the contact force between rock masses according to the material parameters and motion characteristics of the rock mass landslide, and updates the resultant acceleration of each block in real time; Surge simulation module: The simulation calculation of water body and surge wave adopts smooth particle dynamics method. The flow field water body is discretized into a series of fluid particles and calculated by Navier-Stokes equation in integral format. The time step update module adds the calculated coupled acceleration to the rock mass landslide simulation module and the surge simulation module, updates the velocity and coordinate displacement of the entire device, and repeats the above module update process until the calculation time step requirement is met; The control module is connected to the modeling module, the pre-processing module, the rock mass landslide simulation module, the surge simulation module and the time step updating module to control their operation.

8. The multi-block rock mass landslide surge risk assessment system according to claim 7, characterized in that: Also includes: The area division module divides the research area of ​​the actual project, generates a numerical model in STL format, imports it into the pre-processing software, performs discretization, and generates node coordinates for export to subsequent modules.

9. The multi-block rock mass landslide surge risk assessment system according to claim 8, characterized in that: Also includes: The post-processing module is used to output the calculation results in the form of data files, and directly process them using post-processing software to obtain visual calculation results.

Images