Landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling

Abstract

The invention discloses a landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling, and the method comprises the steps: building an SPH-DEM landslide surge calculation model according to a geological profile map, and setting boundary conditions; calculating flows used for water SPH calculation, landslide mass DEM calculation and rock-soil-water flow interaction SPH-DEM calculation are created in the GPU, and particle interaction force is asynchronously calculated in the corresponding calculating flows; asynchronously updating the acceleration, the position, the speed, the angle and the angular speed of the landslide mass DEM particles and the speed, the position and the density of the water body SPH particles in the GPU; and repeating the steps (2) and (3) to achieve a set termination condition, and outputting a calculation result. According to the method, on the basis of GPU asynchronous parallel calculation, the problem that consumed time is too long in SPH-DEM landslide surge calculation is solved.

Description

Technical Field

[0001] The invention relates to a landslide surge calculation method in a geotechnical engineering neighborhood, and in particular to a landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling. Technical Background

[0002] Landslide surge is a secondary disaster caused by landslides from mountains along the coast impacting water after the reservoir is built and filled with water. It is extremely harmful to the coast of the reservoir and is a difficult problem in the field of geotechnical engineering. Rapid and accurate simulation of landslide surge has important scientific significance and engineering value.

[0003] The landslide surge problem involves the deformation, fragmentation and movement of the landslide body, the flow of water and the propagation of surges, and the complex rock-soil-water flow interaction of the landslide body impacting the water body at high speed. The traditional grid method will cause grid distortion when calculating the large deformation of the landslide body and water body, which seriously reduces the calculation accuracy. The use of dynamic grid technology is prone to calculation divergence and greatly increases the amount of calculation, making the calculation of landslide surge difficult.

[0004] Smoothed Particle Hydrodynamics (SPH) is a Lagrangian meshless method for fluid dynamics simulation, which is widely used in the simulation of water flow and waves; Discrete Element Method (DEM) is a meshless method for simulating geotechnical materials, which is widely used in the field of geotechnical engineering. The SPH-DEM coupling method combines the advantages of the two methods to solve the large deformation and damage problems involved in geotechnical-water flow problems, but it also brings the disadvantage of large computational complexity, especially in high-resolution and large-scale simulations, the calculation time is too long.

[0005] Nowadays, GPU has become a more general parallel processor and is widely used in high-performance computing. Since the force and motion calculations of each particle in DEM and SPH are independent, they are naturally suitable for parallel computing. However, due to the different steps and methods of SPH calculation, DEM calculation and SPH-DEM calculation, it is impossible to perform parallel calculations synchronously. For example, it is generally necessary to complete DEM parallel calculation first, then perform SPH parallel calculation, and finally perform SPH-DEM parallel calculation. In the process of waiting for other parallel calculations to complete, GPU computing power will be wasted. Therefore, how to improve the computing efficiency of landslide surge based on GPU is a technical problem that needs to be solved urgently. Summary of the invention

[0006] Purpose of the invention: To solve the above technical problems, the present invention provides a landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling; by adopting the DEM method to simulate the movement of the landslide body, by using the SPH method to simulate the movement of the water flow, according to the SPH-DEM coupling method to simulate the interaction between the landslide body and the water flow, and based on GPU asynchronous parallelism, the calculation efficiency of the SPH-DEM landslide surge is improved, and the problem of long calculation time in large-scale SPH-DEM landslide surge calculation is solved.

[0007] Technical solution: The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling of the present invention includes the following steps:

[0008] (1) Establishing the SPH-DEM landslide surge calculation model based on the geological profile, including steps (1.1)-(1.4):

[0009] (1.1) Setting the boundary conditions of the DEM method: the sliding surface, ground and bottom of the landslide are set as boundaries to restrict the movement of the landslide body, and no boundary is set on the landslide surface;

[0010] (1.2) Setting DEM particle distribution: Filling DEM particles between the sliding surface and the landslide surface, discretizing the rock and soil materials into DEM particles, and realizing the simulation of the landslide body;

[0011] (1.3) Setting the boundary conditions of the SPH method: setting the bottom of the water, the ground and the sliding surface of the landslide as solid wall boundaries to restrict the movement of the water body;

[0012] (1.4) Set the SPH particle distribution: set the closed area enclosed by the water surface and the bottom of the water as the water body, and fill the water body with SPH particles.

[0013] (2) Create a landslide calculation flow in the GPU to calculate the normal and tangential forces between the DEM particles that make up the landslide body, a fluid calculation flow to calculate the pressure, density change rate and acceleration of the SPH particles that make up the water body, and a rock-soil-water flow interaction calculation flow to calculate the drag force and buoyancy of the SPH particles on the DEM particles in the rock-soil-water flow interaction, as well as the reaction force of the SPH particles on the DEM particles.

[0014] (2.1.1) Traverse all DEM particles in parallel in the GPU, and use one thread to parallelly calculate the normal force of the surrounding DEM particles on the DEM particle:

[0015]

[0016] In the formula, F n is the normal force between particles, K n is the normal stiffness coefficient of the setting, η nis the normal damping coefficient, δ n is the normal overlap of the particles, is the normal relative velocity between particles.

[0017] (2.1.2) The tangential force between DEM particles is calculated using the following formula:

[0018]

[0019] In the formula, F t is the tangential force between particles, K t is the set tangential stiffness coefficient, η t is the tangential damping coefficient, δ t is the tangential overlap of the particles, is the tangential relative velocity between particles.

[0020] (2.1.3) Calculate the pressure of the SPH particles that make up the water body:

[0021]

[0022] Where p is the SPH particle pressure, γ is the input sound velocity coefficient, c0 is the input initial sound velocity, ρ is the current density of the particle, and ρ0 is the input initial density.

[0023] (2.1.4) Calculate the density change rate of the SPH particles that make up the water body:

[0024]

[0025] Where, subscript i is the particle number, subscript j is the number of the remaining particles in the i support domain, is the density change rate, ρ i is the density, u i is the velocity of particle i, u j is the velocity of particle j, is the kernel function gradient of the relative position of particles i and j, V j is the particle volume, δ is the coefficient, h is the smoothing length, c0 is the input initial sound velocity, ψ ij is the density dissipation coefficient:

[0026]

[0027] In the formula, r ji is the relative position vector between particles i and j, is the regularized density of particle i, is the regularized density of particle j, calculated using the following formula:

[0028]

[0029]

[0030] Where, L i is the regularization matrix of particle i, is the Kronecker product, is the kernel function gradient of the relative position of particles i and j;

[0031] (2.1.5) Calculate the acceleration of the SPH particles that make up the water body:

[0032]

[0033] Where a is acceleration, ρ is i is the density of particle i, p i is the pressure of particle i, p j is the pressure of particle j, g i is the acceleration of gravity, α is the input artificial viscosity coefficient, ρ0 is the input initial density, h is the smooth length, c0 is the input initial sound velocity, is the kernel function gradient value about the relative position of particles i and j, V j is the volume of particle j, π ij is the coefficient:

[0034]

[0035] In the formula, r ji is the relative position vector between particles i and j, u i is the velocity of particle i, u j is the velocity of particle j;

[0036] (2.1.6) The drag force of SPH particles on DEM particles is calculated from the relative velocity of DEM particles and SPH particles and the porosity of DEM particles:

[0037]

[0038] Where, subscript s represents DEM particles, subscript f represents SPH particles, and f d is the drag force on the DEM particles, μ is the set fluid dynamic viscosity, D is the diameter of the DEM particles, φ is the porosity of the DEM particles, and u f is the velocity of the SPH particle, u s is the velocity of the DEM particles, ρ f is the fluid density, α and β are the set drag coefficients.

[0039] Step (2.1.7) The buoyancy of the SPH particles on the DEM particles is calculated as follows:

[0040]

[0041] In the formula, f b is the buoyancy of the DEM particle, subscript a is the current DEM particle number, subscript i is the SPH particle number, and p i is the pressure of the SPH particle, p a is the pressure of DEM particles, V is the volume of SPH particles, is the kernel function gradient with respect to the relative position of the particle.

[0042] Step (3), calculate the acceleration, position, velocity, angle and angular velocity of the DEM particles that make up the landslide body, and the velocity, position and density of the SPH particles that make up the water body in the GPU:

[0043] (3.1) In the landslide calculation flow, the DEM particles are used as parallel calculation units, all DEM particles are traversed in parallel, and the acceleration of the DEM particles is calculated:

[0044] a=∑(F n +F t )+∑(f d +f b ) (14)

[0045] In the formula, F n is the normal force between DEM particles, F t is the tangential force between DEM particles, f d is the drag force of the SPH particles on the DEM particles, f b is the buoyancy of the SPH particles on the DEM particles.

[0046] The position, velocity, angle and angular velocity of DEM particles are calculated as follows:

[0047]

[0048] v i+1 =v i +δt(a i +a i+1 ) (16)

[0049] θ i+1 =θ i +δt(ω i +ω i+1 ) (17)

[0050] ω i+1 =ω i +δt(α i +α i+1 ) (18)

[0051] Wherein, subscript i represents the current time step, subscript i+1 represents the next time step, δt is the time step, r is the particle position vector, v is the velocity vector, a is the acceleration vector, θ is the angle, ω is the angular velocity, and α is the angular acceleration.

[0052] (3.2) In the fluid calculation flow, for all SPH particles that make up the water body, each SPH particle is used as a parallel calculation unit, all SPH particles are traversed in parallel, and the velocity, position and density of the fluid SPH particles are calculated in parallel:

[0053]

[0054] v i+1 =v i +δt(a i +a i+1 ) (20)

[0055]

[0056] In the formula, subscript i represents the current time step, subscript i+1 represents the next time step, δt is the time step, r is the particle position vector, v is the velocity vector, a is the acceleration vector, ρ is the density, is the density change;

[0057] (3.3) If the current calculation time t = t out , where t out For the output interval time set in the calculation, the position, velocity, angle and angular velocity of the DEM particles of the current landslide body, and the position, velocity and pressure of the SPH particles that make up the water body are output;

[0058] (4) Set the current calculation time to t = t + δt, where δt is the time step. When t> t max Then the calculation ends, where t max is the set calculation time. When the calculation time is less than t≤t max Then return to step (3) to continue the loop calculation;

[0059] (5) Output the final calculated position, velocity, angle and angular velocity of the DEM particles that make up the landslide body, and the position, velocity and pressure of the SPH particles that make up the water body.

[0060] In step (1.3), the water surface is set as the free surface boundary.

[0061] In step (2), the fluid mechanics parameters are set including ρ0 initial density of the fluid, ν dynamic viscosity coefficient of the fluid; the basic parameters required for SPH calculation are set including c0 fluid sound velocity, α artificial viscosity coefficient, δ density dissipation coefficient, h smooth length, γ sound velocity coefficient. The landslide body parameters are set including μ friction coefficient, c cohesion, φ internal friction angle, ρ density; the parameters required for DEM calculation are set including E n Normal stiffness coefficient, E t Tangential stiffness coefficient, η n Normal damping coefficient, η t Tangential damping coefficient, η global damping coefficient, ξ restitution coefficient.

[0062] In step (2), set the calculation time: t out is the output interval time in the calculation, t max is the calculation time of coupling effect, and δt is the calculation time step.

[0063] During calculation, all flows such as landslide calculation flow, fluid calculation flow and rock-soil-water flow interaction calculation flow start calculation at the same time, and each flow performs internal parallel calculation steps asynchronously independently of other flows, and calculations of different flows are carried out simultaneously without interfering with each other.

[0064] In step (2.1.1), the DEM particles are used as parallel computing units, and all DEM particles are traversed in parallel in the GPU. Each DEM particle uses one thread to parallelly calculate the normal force of the surrounding particles on the DEM particle.

[0065] In step (2), the SPH particles are used as parallel computing units to traverse all SPH particles and calculate the pressure, density change rate and acceleration of the SPH particles that make up the water body.

[0066] (2.1.8) Traverse all the SPH particles that make up the water body, use the SPH particles as parallel computing units for parallel computing, and calculate the reaction force of the SPH particles that make up the water body on the DEM particles:

[0067]

[0068] Where, subscript i is the current SPH particle number, subscript b is the DEM particle number, and f is the i is the reaction force of the SPH particles on the DEM particles, m is the mass of the SPH particles, ρ is the density of the SPH particles, and f bb is the buoyancy of DEM particle b, f bd is the drag force on DEM particle b, W is the kernel function of the relative position of the particle, S b is the coefficient:

[0069]

[0070] In the formula, S b The coefficient is calculated from all SPH particles in the DEM particle influence domain, V b is the volume of DEM particles, W jb is the kernel function of the relative position between SPH particle j and DEM particle b.

[0071] In step (3.2), each SPH particle constituting the water body is taken as a parallel computing unit, all SPH particles are traversed, and each SPH particle uses one thread to parallelly calculate the sum of the forces exerted by the remaining SPH particles on the SPH particle.

[0072] Working principle: The present invention is based on the landslide surge calculation method of GPU asynchronous parallel SPH-DEM coupling, which uses the DEM method to simulate the movement of the landslide body, and the SPH method to simulate the movement of the water flow, and the SPH-DEM coupling method is used to simulate the interaction between the landslide body and the water flow. The present invention places the largest part of the landslide surge calculation in the GPU, uses the flow to perform asynchronous operations on the landslide surge calculation, and the fluid calculation flow, the landslide calculation flow, and the rock-soil-water flow interaction force calculation flow are calculated independently. When calculating the landslide surge, each flow performs internal parallel calculation steps independently of other flows, and the calculations of different flows do not interfere with each other. All flows are calculated asynchronously in the GPU at the same time.

[0073] The landslide surge calculation process of the present invention is as follows: (1) inputting a landslide surge geological profile, establishing an SPH-DEM landslide surge calculation model, and inputting calculation parameters; (2) creating a calculation flow in a GPU for water body SPH calculation, landslide body DEM calculation, and rock-soil-water flow interaction SPH-DEM calculation, and synchronously calculating the particle interaction force in the corresponding calculation flow; (3) updating the acceleration, position, velocity, angle, and angular velocity of the landslide body DEM particles, and the velocity, position, and density of the water body SPH particles in the GPU; (4) repeating the calculation steps (2) and (3), reaching the set termination condition, and outputting the position, velocity, angle, and angular velocity of the DEM particles constituting the landslide body, and the position, velocity, and pressure of the SPH particles constituting the water body. The present invention adopts GPU asynchronous parallel calculation, gives full play to the advantages of GPU in large-scale calculation, and significantly improves the calculation efficiency of rock-soil-water flow interaction.

[0074] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0075] (1) The present invention adopts an asynchronous parallel method to avoid the problem of long calculation time in serial calculation and insufficient utilization of GPU computing resources in non-asynchronous parallel calculation, thereby improving the efficiency of parallel calculation of SPH-DEM landslide surge.

[0076] (2) The present invention shortens the time consumption of landslide surge calculation based on SPH-DEM by using GPU, laying a foundation for the engineering application of SPH-DEM landslide surge model. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 This is a flow chart of the landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling of the present invention;

[0078] Figure 2 A cross-sectional view of a landslide surge geological map input in an embodiment of the present invention;

[0079] Figure 3 This is a diagram of the SPH-DEM landslide surge calculation model according to an embodiment of the present invention;

[0080] Figure 4 A diagram showing the formation process of a landslide surge according to an embodiment of the present invention;

[0081] in, Figure 4 (a) is the result diagram when the calculation reaches 25s;

[0082] Figure 4 (b) is the result diagram when the calculation reaches 40s;

[0083] Figure 4 (c) is the result diagram when the calculation reaches 69s. DETAILED DESCRIPTION

[0084] like Figure 1 As shown, the landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling of the present invention includes the following steps:

[0085] (1) Establishing the SPH-DEM landslide surge calculation model based on the geological profile, including steps (1.1)-(1.4):

[0086] (1.1) Setting the boundary conditions of the DEM method: the sliding surface, ground and bottom of the landslide are set as boundaries to restrict the movement of the landslide body, and no boundary is set on the landslide surface;

[0087] (1.2) Setting the DEM particle distribution, filling DEM particles between the sliding surface and the landslide surface, discretizing the rock and soil materials into DEM particles, and realizing the simulation of the landslide body;

[0088] (1.3) Set the boundary conditions of the SPH method: the bottom of the water, the ground and the sliding surface of the landslide are regarded as solid wall boundaries to restrict the movement of the water body, and the water surface is regarded as a free surface boundary;

[0089] (1.4) Setting the SPH particle distribution: The closed area enclosed by the water surface and the bottom is set as the water body, and the SPH particles are filled in the water body;

[0090] (2) In the GPU, a landslide calculation flow is created to calculate the normal and tangential forces between the DEM particles that make up the landslide body, a fluid calculation flow is created to calculate the pressure, density change rate and acceleration of the SPH particles that make up the water body, and a rock-soil-water flow interaction calculation flow is created to calculate the drag force and buoyancy of the SPH particles on the DEM particles in the rock-soil-water flow interaction, as well as the reaction force of the SPH particles on the DEM particles.

[0091] Set the fluid mechanics parameters including ρ0 initial density of the fluid, ν dynamic viscosity coefficient of the fluid; set the basic parameters required for SPH calculation including c0 fluid sound velocity, α artificial viscosity coefficient, δ density dissipation coefficient, h smooth length, γ sound velocity coefficient. Set the landslide parameters including μ friction coefficient, c cohesion, φ internal friction angle, ρ density; set the parameters required for DEM calculation including E n Normal stiffness coefficient, E t Tangential stiffness coefficient, η n Normal damping coefficient, η t Tangential damping coefficient, η global damping coefficient, ξ restitution coefficient.

[0092] Set calculation time: t out The output interval of the running file, t max Calculation time of coupling effect, δt is the calculation time step.

[0093] During calculation, all streams start calculating at the same time, and each stream performs internal parallel calculation steps asynchronously independent of other streams. The calculations of different streams do not interfere with each other and are performed simultaneously. The process is:

[0094] (2.1.1) Traverse all DEM particles in parallel in the GPU, and use one thread to parallelly calculate the normal force of the surrounding DEM particles on the DEM particle:

[0095]

[0096] In the formula, F n is the normal force between particles, K n is the normal stiffness coefficient of the setting, η n is the normal damping coefficient, δ n is the normal overlap of the particles, is the normal relative velocity between particles.

[0097] (2.1.2) The tangential force between DEM particles is calculated using the following formula:

[0098]

[0099] In the formula, F t is the tangential force between particles, Kt is the set tangential stiffness coefficient, η t is the tangential damping coefficient, δ t is the tangential overlap of the particles, is the tangential relative velocity between particles.

[0100] (2.1.3) Calculate the pressure of the SPH particles that make up the water body:

[0101]

[0102] Where p is the SPH particle pressure, γ is the input sound velocity coefficient, c0 is the input initial sound velocity, ρ is the current density of the particle, and ρ0 is the input initial density.

[0103] (2.1.4) Calculate the density change rate of the SPH particles that make up the water body:

[0104]

[0105] Where, subscript i is the particle number, subscript j is the number of the remaining particles in the i support domain, is the density change rate, ρ i is the density, u i is the velocity of particle i, u j is the velocity of particle j, is the kernel function gradient of the relative position of the particle, V j is the particle volume, δ is the coefficient, h is the smoothing length, c0 is the input initial sound velocity, ψ ij is the density dissipation coefficient:

[0106]

[0107] In the formula, r ji is the relative position vector between particles i and j, is the regularized density of particle i, is the regularized density of particle j, calculated using the following formula:

[0108]

[0109]

[0110] Where, L i is the regularization matrix of particle i, is the Kronecker product, is the kernel function gradient of the relative position of particles i and j;

[0111] (2.1.5) Calculate the acceleration of the SPH particles that make up the water body:

[0112]

[0113] Where a is acceleration, ρ is i is the density of particle i, p i is the pressure of particle i, p j is the pressure of particle j, g i is the acceleration of gravity, α is the input artificial viscosity coefficient, ρ0 is the input initial density, h is the smooth length, c0 is the input initial sound velocity, is the kernel function gradient value about the relative position of particles i and j, V j is the volume of particle j, π ij is the coefficient:

[0114]

[0115] In the formula, r ji is the relative position vector between particles i and j, u i is the velocity of particle i, u j is the velocity of particle j.

[0116] (2.1.6) Traverse all the DEM particles that make up the landslide body and perform parallel calculations with DEM particles as parallel calculation units. The drag force of SPH particles on DEM particles is calculated by the relative velocity between DEM particles and SPH particles and the porosity of DEM particles:

[0117]

[0118] Where, subscript s represents DEM particles, subscript f represents SPH particles, and f d is the drag force on the DEM particles, μ is the set fluid dynamic viscosity, D is the diameter of the DEM particles, φ is the porosity of the DEM particles, and u f is the velocity of the SPH particle, u s is the velocity of the DEM particles, ρ f is the fluid density, α and β are the set drag coefficients.

[0119] (2.1.7) The buoyancy of SPH particles on DEM particles is calculated as follows:

[0120]

[0121] In the formula, f b is the buoyancy of the DEM particle, subscript a is the current DEM particle number, subscript i is the SPH particle number, and p i is the pressure of the SPH particle, p a is the pressure of DEM particles, V is the volume of SPH particles, is the kernel function gradient with respect to the relative position of the particle.

[0122] (2.1.8) Traverse all the SPH particles that make up the water body, use the SPH particles as parallel computing units for parallel computing, and calculate the reaction force of the SPH particles that make up the water body on the DEM particles:

[0123]

[0124] Where, subscript i is the current SPH particle number, subscript b is the DEM particle number, and f is the i is the reaction force of the SPH particles on the DEM particles, m is the mass of the SPH particles, ρ is the density of the SPH particles, and f bb is the buoyancy of DEM particle b, f bd is the drag force on DEM particle b, W is the kernel function of the relative position of the particle, S b is the coefficient:

[0125]

[0126] In the formula, S b The coefficient is calculated from all SPH particles in the DEM particle influence domain, V b is the volume of DEM particles, W jb is the kernel function of the relative position between SPH particle j and DEM particle b.

[0127] (3) The position, velocity, angle and angular velocity of the DEM particles that make up the landslide body and the velocity, position and density of the SPH particles that make up the water body are calculated in the GPU:

[0128] (3.1) In the landslide calculation flow, the DEM particles are used as parallel calculation units, all DEM particles are traversed in parallel, and the acceleration of the DEM particles is calculated:

[0129] a=Σ(F n +F t )+∑(f d +f b ) (14)

[0130] In the formula, F n is the normal force between DEM particles, F t is the tangential force between DEM particles, f d is the drag force of the SPH particles on the DEM particles, f b is the buoyancy of the SPH particles on the DEM particles.

[0131] The position, velocity, angle and angular velocity of DEM particles are calculated as follows:

[0132]

[0133] v i+1 =v i +δt(a i +a i+1 ) (16)

[0134] θ i+1 =θ i +δt(ω i +ω i+1 ) (17)

[0135] ω i+1 =ω i +δt(α i +α i+1 ) (18)

[0136] Wherein, subscript i represents the current time step, subscript i+1 represents the next time step, δt is the time step, r is the particle position vector, v is the velocity vector, a is the acceleration vector, θ is the angle, ω is the angular velocity, and α is the angular acceleration.

[0137] (3.2) In the fluid calculation flow, for all SPH particles that make up the water body, each SPH particle is used as a parallel calculation unit, all SPH particles are traversed in parallel, and the velocity, position and density of the fluid SPH particles are calculated in parallel:

[0138]

[0139] v i+1 =v i +δt(a i +a i+1 ) (20)

[0140]

[0141] In the formula, subscript i represents the current time step, subscript i+1 represents the next time step, δt is the time step, r is the particle position vector, v is the velocity vector, a is the acceleration vector, ρ is the density, is the density change;

[0142] (3.3) If the current calculation time t = t out , where t out For the output time in the set calculation, the position, velocity, angle and angular velocity of the DEM particles of the current landslide body, and the position, velocity and pressure of the SPH particles that make up the water body are output;

[0143] (4) Set the current calculation time to t = t + δt, where δt is the time step. When t> t max Then the calculation ends, where t maxis the set calculation time, the current calculation time is less than t≤t max Then return to step (3) to continue the loop calculation;

[0144] (5) Output the final calculated position, velocity, angle and angular velocity of the DEM particles that make up the landslide body, and the position, velocity and pressure of the SPH particles that make up the water body.

[0145] Example:

[0146] Taking the calculation of rock-soil-water interaction of a landslide-generated surge as an example, the landslide surge calculation process based on GPU asynchronous parallel SPH-DEM coupling in the present invention is as follows:

[0147] 1. Input the landslide geological profile, such as Figure 2 , input calculation parameters. In the GPU, create a landslide calculation flow to calculate the normal force and tangential force between the DEM particles that make up the landslide body, a fluid calculation flow to calculate the pressure, density change rate and acceleration of the SPH particles that make up the water body, and a rock-soil-water flow interaction calculation flow to calculate the drag force and buoyancy of the SPH particles on the DEM particles in the rock-soil-water flow interaction, as well as the reaction force of the SPH particles on the DEM particles.

[0148] 2. Set the model boundary conditions: Set the landslide sliding surface, ground and bottom of the water as DEM boundaries, and do not set boundaries on the landslide surface. Set the bottom of the water, ground and landslide sliding surface as SPH solid wall boundaries, and the water surface as a free surface boundary.

[0149] 3. Set particle distribution: The closed area surrounded by the water surface and the bottom of the water is the water body, and the SPH particles are filled in the water body. DEM particles are filled between the landslide sliding surface and the landslide surface, and the rock and soil materials are discretized into DEM particles, such as Figure 3 shown.

[0150] 4. Perform landslide surge calculations, calculate the interaction force between the DEM particles of the landslide body, and obtain the movement process of the landslide body; calculate the interaction between the SPH particles of the water body, and obtain the process of water flow and surge propagation; calculate the interaction force between the SPH particles and the DEM particles, and obtain the process of the landslide body impacting the water body and generating surges. The calculation is performed in the corresponding flow, all flows start calculations at the same time, and each flow performs internal parallel calculation steps asynchronously independent of other flows, and the calculations of different flows do not interfere with each other and are performed simultaneously.

[0151] 5. Calculate the position, speed, angle and angular velocity of the DEM particles that make up the landslide body, and the speed, position and density of the SPH particles that make up the water body in the GPU. The operation time set in this embodiment is 69s, and the current calculation results are output at 25s and 40s. According to the calculation results, the DEM particles that make up the landslide body and the SPH particles that make up the water body are plotted on the coordinate axis. The colors of the DEM particles and the SPH particles are set according to the speed, and the landslide surge formation process diagram is obtained as shown below: Figure 4 shown.

[0152] 6. From Figure 4 The positions of the DEM particles that make up the landslide and the SPH particles that make up the water body can be seen in the figure. The lighter the color, the faster the particle speed. The DEM particles are impacting the water on the left, and the SPH particles are on the right. The arrow in the water body represents the direction of movement of the SPH particles.

[0153] 7. If Figure 4 (a), at t = 25s, the sliding body enters the water and the first wave of the surge begins to form. Since the sliding body is wedge-shaped when entering the water and the sliding speed is slow, the pressure on the water body is small and no sharp crest is formed. Figure 4 (b), at t = 40s, the surge propagates to the center axis of the river channel, and the surge height is 34.5m. At this time, the water surface on the opposite bank begins to rise and enters the climbing stage. Figure 4 (c) At t = 69s, the calculation ends, and the climbing height on the opposite bank reaches the maximum, which is 85.05m.

Claims

1. A landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling, characterized by: The following steps are involved: (1) Establishing the SPH-DEM landslide surge calculation model based on the geological profile, including steps (1.1)-(1.4): (1.1) Setting the boundary conditions of the DEM method: the landslide sliding surface, ground and bottom are set as boundaries; (1.2) Setting DEM particle distribution: filling DEM particles between the sliding surface and the landslide surface, discretizing the geotechnical materials into DEM particles; (1.3) Set the boundary conditions of the SPH method: set the bottom of the water, the ground and the sliding surface of the landslide as solid wall boundaries; (1.4) Setting the SPH particle distribution: The closed area enclosed by the water surface and the bottom is set as the water body, and the SPH particles are filled in the water body; (2) Create a landslide calculation flow in the GPU to calculate the normal and tangential forces between the DEM particles that make up the landslide body, a fluid calculation flow to calculate the pressure, density change rate and acceleration of the SPH particles that make up the water body, and a rock-soil-water flow interaction calculation flow to calculate the drag force and buoyancy of the SPH particles on the DEM particles in the rock-soil-water flow interaction, as well as the reaction force of the SPH particles on the DEM particles: (2.1.1) Traverse all DEM particles in parallel in the GPU, and use one thread to parallelly calculate the normal force of the surrounding DEM particles on the DEM particle: In the formula, F n is the normal force between particles, K n is the normal stiffness coefficient of the setting, η n is the normal damping coefficient, δ n is the normal overlap of the particles, is the normal relative velocity between particles; (2.1.2) The tangential force between DEM particles is calculated using the following formula: In the formula, F t is the tangential force between particles, K t is the set tangential stiffness coefficient, η t is the tangential damping coefficient, δ t is the tangential overlap of the particles, is the tangential relative velocity between particles; (2.1.3) Calculate the pressure of the SPH particles that make up the water body: Where p is the SPH particle pressure, γ is the input sound velocity coefficient, c0 is the input initial sound velocity, ρ is the current density of the particle, and ρ0 is the input initial density; (2.1.4) Calculate the density change rate of the SPH particles that make up the water body: Where, subscript i is the particle number, subscript j is the number of the remaining particles in the i support domain, is the density change rate, ρ i is the density, u i is the velocity of particle i, u j is the velocity of particle j, is the kernel function gradient of the relative position of particles i and j, V j is the volume of particle j, δ is the coefficient, h is the smoothing length, c0 is the initial sound velocity input, ψ ij is the density dissipation coefficient: In the formula, r ji is the relative position vector between particles i and j, is the regularized density of particle i, is the regularized density of particle j, calculated using the following formula: Where, L i is the regularization matrix of particle i, is the Kronecker product, is the kernel function gradient of the relative position of particles i and j; (2.1.5) Calculate the acceleration of the SPH particles that make up the water body: Where a is acceleration, ρ is i is the density of particle i, p i is the pressure of particle i, p j is the pressure of particle j, g i is the acceleration of gravity, α is the input artificial viscosity coefficient, h is the smooth length, c0 is the input initial sound velocity, ρ0 is the input initial density, is the kernel function gradient value about the relative position of particles i and j, V j is the volume of particle j, π ij is the coefficient: In the formula, r ji is the relative position vector between particles i and j, u i is the velocity of particle i, u j is the velocity of particle j; (2.1.6) The drag force of SPH particles on DEM particles is calculated from the relative velocity of DEM particles and SPH particles and the porosity of DEM particles: Where, subscript s represents DEM particles, subscript f represents SPH particles, and f d is the drag force on the DEM particles, μ is the set fluid dynamic viscosity, D is the diameter of the DEM particles, φ is the porosity of the DEM particles, and u f is the velocity of the SPH particle, u s is the velocity of the DEM particles, ρ f is the fluid density, α and β are the set drag coefficients; (2.1.7) The buoyancy of SPH particles on DEM particles is calculated as follows: In the formula, f b is the buoyancy of the DEM particle, subscript a is the current DEM particle number, subscript i is the SPH particle number, and p i is the pressure of the SPH particle, p a is the pressure of the DEM particle, V is the volume of the SPH particle, ▽W is the kernel function gradient for the relative position of the particle; (3) Calculate the acceleration, position, velocity, angle and angular velocity of the DEM particles that make up the landslide body, and the velocity, position and density of the SPH particles that make up the water body in the GPU: (3.1) In the landslide calculation flow, calculate the acceleration of the DEM particles: a=∑(F n +F t )+∑(f d +f b ) (14) In the formula, F n is the normal force between DEM particles, F t is the tangential force between DEM particles, f d is the drag force of the SPH particles on the DEM particles, f b is the buoyancy of the SPH particles on the DEM particles; The position, velocity, angle and angular velocity of DEM particles are calculated as follows: v i+1 =v i +δt(a i +a i+1 ) (16) i i+1 =θ i +δt(ω i +oh i+1 ) (17) oh i+1 =ω i +δt(a i +a i+1 ) (18) Where, subscript i represents the current time step, subscript i+1 represents the next time step, δt is the time step, r is the particle position vector, v is the velocity vector, a is the acceleration vector, θ is the angle, ω is the angular velocity, and α is the angular acceleration; (3.2) In the fluid calculation flow, each SPH particle is used as a parallel calculation unit, all SPH particles are traversed in parallel, and the speed, position and density of the SPH particles are calculated in parallel: v i+1 =v i +δt(a i +a i+1 ) (20) In the formula, subscript i represents the current time step, subscript i+1 represents the next time step, δt is the time step, r is the particle position vector, v is the velocity vector, a is the acceleration vector, ρ is the density, is the density change; (3.3) If the current calculation time t = t out , where t out For the output interval time set in the calculation, the position, velocity, angle and angular velocity of the DEM particles of the current landslide body, and the position, velocity and pressure of the SPH particles that make up the water body are output; (4) Set the current calculation time to t = t + δt, where δt is the calculation time step. When t> t max Then the calculation ends, where t max is the set calculation time. When the calculation time is less than t≤t max Then return to step (3) to continue the loop calculation; (5) Output the final calculated position, velocity, angle and angular velocity of the DEM particles that make up the landslide body, and the position, velocity and pressure of the SPH particles that make up the water body.

2. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized by: In step (1.3), the water surface is set as the free surface boundary.

3. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized by: In step (2), the fluid mechanics parameters, the parameters required for SPH calculation, the landslide body parameters and the parameters required for DEM calculation are set.

4. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized in that: In step (2), set the calculation time: t out is the output interval time in the calculation, t max is the calculation time of coupling effect, and δt is the calculation time step.

5. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized in that: In step (2.1.1), the DEM particles are used as parallel computing units, and all DEM particles are traversed in parallel in the GPU. Each DEM particle uses one thread to parallelly calculate the normal force of the surrounding particles on the DEM particle.

6. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized by: In step (2), the SPH particles are used as parallel computing units to traverse all SPH particles and calculate the pressure, density change rate and acceleration of the SPH particles that make up the water body.

7. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized by: In step (2), set the fluid mechanics parameters: ρ0 initial density of the fluid, ν dynamic viscosity coefficient of the fluid; set the parameters required for SPH calculation: c0 fluid sound velocity, α artificial viscosity coefficient, δ density dissipation coefficient, h smooth length, γ sound velocity coefficient; set the landslide body parameters: μ friction coefficient, c cohesion, φ internal friction angle, ρ density; set the parameters required for DEM calculation: E n Normal stiffness coefficient, E t Tangential stiffness coefficient, η n Normal damping coefficient, η t Tangential damping coefficient, η global damping coefficient, ξ restitution coefficient.

8. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized by: In step (2.1.8), all SPH particles constituting the water body are traversed, and parallel calculation is performed with the SPH particles as parallel calculation units to calculate the reaction force of the SPH particles on the DEM particles.

9. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 8 is characterized by: In step (2.1.8), the reaction force of the SPH particles constituting the water body on the DEM particles is calculated: Where, subscript i is the current SPH particle number, subscript b is the DEM particle number, and f is the i is the reaction force of the SPH particles on the DEM particles, m is the mass of the SPH particles, ρ is the density of the SPH particles, and f bb is the buoyancy of DEM particle b, f bd is the drag force on DEM particle b, W is the kernel function of the relative position of the particle, S b is the coefficient: In the formula, S b The coefficient is calculated from all SPH particles in the DEM particle influence domain, V b is the volume of DEM particles, W jb is the kernel function of the relative position between SPH particle j and DEM particle b.

10. The landslide surge calculation method based on GPU asynchronous parallel SPH-DEM coupling according to claim 1 is characterized by: In step (3.1), in the landslide calculation flow, DEM particles are used as parallel calculation units, all DEM particles are traversed in parallel, and the acceleration of the DEM particles is calculated.

Images