Simulation method of rock slope failure process and related equipment

Abstract

The present disclosure provides a simulation method for failure process of rock slope and related equipment. Specifically, the simulation method comprises: constructing a slope model of a target slope; wherein the slope model comprises a plurality of particles respectively assigned with different parameter values, and the plurality of particles belong to at least two types of particles; according to a preset time step and the parameter values, cyclically calculating a damage state and a target parameter value of the plurality of particles until a predetermined condition; and outputting the damage state and the target parameter value of the plurality of particles at least once. The simulation method of the present disclosure can quickly solve and evaluate the cracking, sliding and large deformation of the rock slope.

Description

Technical Field

[0001] This disclosure relates to the field of simulation and prediction technology, and in particular to a simulation method and related equipment for the failure process of rock slopes. Background Technology

[0002] The failure process of rock slopes includes crack initiation, propagation, and contact slip. Crack initiation, propagation, and contact slip are particularly important for rock slope stability assessment and slope protection engineering. The cracking process in rock mass is controlled by the physical and mechanical properties of the rock matrix and pre-existing discontinuities, as well as their spatial distribution. Therefore, research on the formation laws and sliding characteristics of cracks in discontinuous rock slopes is of great significance for rock slope stability assessment and slope protection engineering. Summary of the Invention

[0003] In view of this, the purpose of this disclosure is to propose a simulation method and related equipment for the failure process of rock slopes.

[0004] To achieve the above objectives, this disclosure provides a method for simulating the failure process of rock slopes, including:

[0005] Construct a slope model of the target slope; wherein the slope model includes multiple particles that are assigned different parameter values, and the multiple particles belong to at least two types of particles;

[0006] Based on the preset time step and the parameter values, the damage state and target parameter values ​​of the multiple particles are calculated iteratively until the predetermined conditions are met.

[0007] Output at least once the damage state of the plurality of particles and the target parameter value; wherein,

[0008] The calculation of the damage state and target parameter values ​​of the plurality of particles includes:

[0009] Determine the target particle from among the plurality of particles;

[0010] In response to determining that the target particle's historical damage state is that it is an undamaged particle, the nearest neighbor particle of the target particle is determined based on the previous parameter value and a preset compact support region.

[0011] Based on the Mohr-Coulomb theory, the damage state of the target particle is determined.

[0012] Based on historical parameter values, the current damage state of the target particle, and the historical damage state of the nearest neighboring particle, determine the corresponding preset target parameter control equation and calculate the target parameter change rate of the current target particle;

[0013] Based on a preset time step, the previous parameter value, and the target parameter change rate, the current target parameter value of the target particle is calculated; wherein, the current target parameter value is used to update the previous parameter value.

[0014] Based on the same inventive concept, this disclosure also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the simulation method as described in any of the foregoing embodiments.

[0015] Based on the same inventive concept, this disclosure also provides a non-transitory computer-readable storage medium that stores computer instructions for causing a computer to perform any of the simulation methods described above.

[0016] As can be seen from the above, the simulation method and related equipment for the failure process of rock slopes provided in this disclosure fully restore the characteristics of the target slope by constructing slope models with different types of particles. Using the Mohr-Coulomb theory, the damage state of the particles is determined, thereby enabling the determination of the corresponding target parameter control equations based on the damage state and the calculation of the target parameter values. Finally, the cracking, sliding, and large deformation of the rock slope can be quickly solved and evaluated. Furthermore, the target parameter control equations incorporate the contact force of the particles to characterize the discontinuities caused by the damaged particles, making the simulation results more reliable. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in this disclosure or related technologies, the accompanying drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the accompanying drawings described below are only embodiments of this disclosure. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 This diagram illustrates the principle of a rock slope fracture contact frame according to an embodiment of the present disclosure.

[0019] Figure 2 Show Figure 1 Enlarged view of the middle section structure;

[0020] Figure 3 This diagram illustrates a DSPH model of a bedding rock slope according to an embodiment of the present disclosure.

[0021] Figure 4(a) shows a schematic diagram of the failure process of a bedding rock slope provided in an embodiment of the present disclosure;

[0022] Figure 4(b) shows a schematic diagram of displacement deformation corresponding to Figure 4(a);

[0023] Figure 4(c) shows a schematic diagram of the evolution of the maximum principal stress corresponding to Figure 4(a);

[0024] Figure 4(d) shows a schematic diagram of the experimental results of a bedding rock slope provided in an embodiment of this disclosure;

[0025] Figure 4(e) shows a schematic diagram of the FDEM simulation results of a bedding rock slope provided in an embodiment of the present disclosure;

[0026] Figure 5 This diagram illustrates a DSPH model of a reverse-dip layered rock slope provided in an embodiment of the present disclosure.

[0027] Figure 6(a) shows a schematic diagram of the failure process of a reverse-dip layered rock slope provided in an embodiment of the present disclosure;

[0028] Figure 6(b) shows a schematic diagram of displacement deformation corresponding to Figure 6(a);

[0029] Figure 6(c) shows a schematic diagram of the evolution of the maximum principal stress corresponding to Figure 6(a);

[0030] Figure 6(d) shows a schematic diagram of experimental results for a reverse-dip layered rock slope provided in an embodiment of this disclosure;

[0031] Figure 6(e) shows a schematic diagram of the simulation results of a reverse-dip layered rock slope using the general discrete element method according to an embodiment of the present disclosure;

[0032] Figure 7 This diagram illustrates a DSPH model of a discontinuous jointed rock slope according to an embodiment of the present disclosure.

[0033] Figure 8(a) shows a schematic diagram of the failure process of a discontinuous jointed rock slope provided in an embodiment of the present disclosure;

[0034] Figure 8(b) shows a schematic diagram of displacement deformation corresponding to Figure 8(a);

[0035] Figure 8(c) shows a schematic diagram of the evolution of the maximum principal stress corresponding to Figure 8(a);

[0036] Figure 9 This diagram illustrates a partial flow chart of a simulation method for the failure process of a rock slope provided in an embodiment of this disclosure.

[0037] Figure 10 This diagram illustrates the structure of an electronic device provided in an embodiment of the present disclosure. Detailed Implementation

[0038] To make the objectives, technical solutions, and advantages of this disclosure clearer, the following detailed description is provided in conjunction with specific embodiments and the accompanying drawings.

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

[0040] To facilitate understanding of the technical solutions disclosed herein, some technical terms involved in this disclosure will be introduced below.

[0041] Smoothed Particle Hydrodynamics (SPH) is a typical meshless Lagrangian particle method.

[0042] Discontinuous SPH (DSPH) is an improved method of SPH for interface problems. Its essence is that when estimating the derivative value of a particle at a certain interface, it is less affected by particles on the opposite side of the interface.

[0043] The Discrete Element Method (DEM) is a method for analyzing particulate discrete materials. It was first applied to the analysis of rock mechanics problems and has since been gradually applied to the fields of bulk materials and powder engineering.

[0044] Discontinuous Deformation Analysis (DDA) is a numerical simulation method applicable to discontinuous media such as rock masses.

[0045] The hybrid discrete finite element method (FDEM) is an advanced finite element-discrete element coupled method that can effectively simulate large deformation and fracture processes of materials.

[0046] As mentioned in the background section, the study of the formation law and sliding characteristics of cracks in discontinuous rock slopes is of great significance for rock slope stability evaluation and slope protection engineering. Numerical simulation is a useful tool for revealing the sliding mechanism of rock slopes. However, traditional finite element methods and finite difference methods, due to mesh limitations, can usually only simulate small deformations of rock slopes and cannot realistically simulate the entire process of rock mass cracking and failure. For example, discrete element methods and discontinuous deformation analysis methods have many micro-parameters that are not practically meaningful, require a lot of manual calibration, and have relatively low computational efficiency, making them difficult to apply in engineering practice. Therefore, there is an urgent need to develop a new numerical simulation method for the crack propagation and contact sliding process of slopes to accurately reflect the characteristics of crack propagation and contact sliding in rock slopes.

[0047] In view of this, the present disclosure provides a simulation method and related equipment for the failure process of rock slopes. The simulation method is based on a discontinuous smooth particle hydrodynamics (DSPH) framework. DSPH has unique advantages in simulating large deformations and discontinuous problems in materials: First, compared with continuous methods such as the finite element method and the finite difference method, DSPH has no mesh constraints and can simulate large deformations of rock slopes; second, compared with discontinuous methods such as the discrete element method and discontinuous deformation analysis methods, the mechanical parameters of DSPH are calculated based on partial differential equations in a continuous domain, with very clear physical meaning and no need for complex microscopic parameter calibration; finally, compared with the traditional SPH method, DSPH effectively solves the problem that the fracture behavior of rock slopes and the frictional contact of the fracture surface cannot be characterized.

[0048] To make the technical solution of this disclosure clearer and easier to understand, the principle of the rock slope fracture contact frame provided in the embodiments of this disclosure will be introduced below with reference to the accompanying drawings.

[0049] Figure 1 This diagram illustrates the principle of a rock slope fracture contact frame according to an embodiment of the present disclosure. Figure 1 As shown, in the initial state of the rock slope model, multiple particles form a continuous domain through virtual bonds. Under tensile and shear failure, some virtual bonds break to form fracture bonds, causing particles located within these fracture bonds to break and become damaged particles. Next, cracks form along the fracture bonds, and the damaged particles detach from the remaining particles, becoming discontinuous particles. Then, these discontinuous particles establish contact with particles within a preset critical distance (including both discontinuous and undamaged particles), and are subjected to contact forces (especially tangential forces), causing the initial continuous domain to become discontinuous. Here, the criterion for determining whether virtual bonds break to form fracture bonds is based on the pseudo-spring method using the Mohr-Coulomb and maximum tensile stress failure criteria.

[0050] Next, we will provide an exemplary illustration of the various particle parameter control equations for the slope model at different fracture stages.

[0051] In some embodiments, each particle in the model is an undamaged particle, and its corresponding control equation is shown below. It should be noted that equation (1) is constructed based on the mass conservation discrete construction of DSPH.

[0052]

[0053] Where ρ is the particle density; m is the particle mass; i is the index of the target particle; j is the index of the nearest neighbor particle; v is the particle velocity; x is the particle position coordinates; N is the total number of nearest neighbor particles; α and β represent the x, y, and z directions in the Cartesian coordinate system, respectively, i.e., the Einstein summation convention, applied to repeated indexing; W ij The kernel function is a smooth kernel; σ represents the Cauchy stress tensor; δ αβ For Dirac function; Π ij Artificial viscosity; T ij For artificial stress; f i α It could be gravity, and the corresponding acceleration could be, for example, 9.81 m / s². 2 .

[0054] Optionally, the nearest neighbor particles in the compactly supported region of the target particle can be determined using a linked list search method or a tree search method. Here, using a linked list search method or a tree search method is beneficial for improving search efficiency, increasing computation speed, and saving computation time. It should be noted that the linked list search method and the tree search method are merely examples, and this disclosure does not specifically limit their application.

[0055] The following embodiments of this disclosure illustrate the calculation methods for smooth kernel function, Cauchy stress, artificial viscosity, and artificial stress.

[0056] In some embodiments, the smoothing kernel function employs a B-spline function (Monaghan, 1994) to reduce instabilities in particle simulations of compressed states, and its form is:

[0057]

[0058] In the formula, R = r ij / h,r ij =x i -x j , representing the distance between target particle i and its nearest neighbor particle j, where λ is the scaling factor of the smooth kernel function space dimension, set to in two-dimensional space. h is the radius of the smooth kernel.

[0059] Here, the compactly supported region can be the radius of a smooth core, for example, 0.25m.

[0060] In some embodiments, the formula for calculating Cauchy stress may be:

[0061] σ αβ =-pδ αβ +S αβ (3)

[0062] In the formula, P is the isotropic pressure and S is the deviatoric stress.

[0063] Optionally,

[0064] In equation (4), ρ0 is the initial density of the target particle, ρ is the current density of the target particle, and E is the elastic modulus.

[0065] In solid mechanics, the basic model generally allows stress to be a function of strain and strain rate. For anisotropic shear stress, if a small displacement is assumed, the stress rate and strain rate are proportional. The proportionality constant is the shear modulus G = E / (2 × (1 + u)), where E is the elastic modulus and u is Poisson's ratio. Based on this, the equation for the stress rate tensor is:

[0066]

[0067] In equation (5), ε γγ =ε xx +ε yy +ε zz , ε αβ The strain rate tensor is defined as:

[0068]

[0069] To align the material information with strain, a continuity equation is derived by introducing the Jaumann rate of change:

[0070]

[0071] γ represents the x, y, and z axes in the Cartesian coordinate system, i.e., the Einstein summation convention, applied to repeated indexing; ω βγ and ω αγ These are torsion rate tensors, which are defined as:

[0072]

[0073] In some embodiments, Π ij Artificial viscosity, its calculation formula is:

[0074]

[0075] In the formula, c0 is the numerical sound velocity of the solid particle, the magnitude of which can be determined by the elastic modulus E and the initial density, for example... v ijIt is speed; β1 and β2 are control parameters predetermined according to the specific problem, which can be found in Examples 1 to 3; This represents the average particle density.

[0076] In some embodiments, T ij Artificial stress, used to improve the numerical stability of DSPH, can be calculated using the following formula:

[0077] T ij =f n ×R ij (10)

[0078]

[0079] In the formula, ΔP is the initial particle spacing. and It is related to the stress state of the rock and soil materials.

[0080] Optionally, and It can be calculated and determined by the following formula:

[0081]

[0082]

[0083] in, and θ t It is calculated by the following formula:

[0084]

[0085]

[0086] Where, σ t xx and σ t yy It is σ in equation (3) αβ The superscripts x and y are different.

[0087] It should be understood that the above description is merely exemplary, and those skilled in the art can select appropriate functions to calculate the parameters in equation (1) according to their needs, and this disclosure does not limit this.

[0088] Next, as Figure 1 As shown, under the action of tensile failure and shear failure, some virtual bonds of the undamaged particles in the model break to form broken bonds, which in turn causes the particles located in the broken bonds to break and form damaged particles.

[0089] For example, the damage to particles within a slope can be determined using the Mohr-Coulomb criterion with tensile truncation. Specifically, the criterion is as follows:

[0090] σ n =1 / 2(σ1+σ3)-1 / 2(σ1-σ3)sinφ≥σ t (16)

[0091] τ n =1 / 2(σ1-σ3)×cosφ≥τ t =c+σ n tanφ (17)

[0092] In the formula: σ n and τ n For tensile stress and shear stress, σ t and τ t For tensile strength and shear strength; σ1 is the maximum principal stress, σ3 is the minimum principal stress, c is the cohesion of the rock, and φ is the angle of internal friction. Where σ xx and σ yy From σ in equation (3) αβ The superscripts are used to represent different x and y values.

[0093] It should be noted that σ n τ n σ t and τ t If the conditions of equations (16) and (17) are satisfied, then particle breakage can be determined.

[0094] It should be understood that the parameter control equation of equation (1) cannot be applied to the presence of broken particles, whether the broken particles are present in the compactly supported domain of the undamaged particles or the particles themselves are broken particles. In view of this, the exemplary embodiments of this disclosure further describe the parameter control equations for the presence of broken particles in the compactly supported domain of the undamaged particles and for the particles themselves being broken particles.

[0095] By introducing a parameter f to improve equation (1), when a local particle satisfies the above destruction criterion, it is considered completely damaged, and f = 0; when a local particle does not satisfy the above destruction criterion, it is considered intact, and f = 1. The governing equations for the improved undamaged particles and completely damaged particles are as follows:

[0096]

[0097] In the formula: U represents an undamaged particle, ζ represents a completely damaged particle; W ij (i.e. W(x) i -x j ,h)) and (Right now ) are the smooth kernel function and its gradient for the undamaged particle (IP); and , and , respectively, the smooth kernel function and its gradient for a completely broken particle; where, Where f is the damage coefficient.

[0098] refer to Figure 1 Each completely damaged particle corresponds to a discontinuous particle (DP) (i.e., a broken particle). Since the fracture is instantaneous, the fracture surface becomes a contact surface, and based on the contact, contact forces are generated between the particles.

[0099] In cases where there are broken particles in the compactly supported domain of undamaged particles, when the fracture surface is transformed into a contact surface, the discontinuous particles in the compactly supported domain of the undamaged particles (i.e., the aforementioned undamaged particles) are no longer applicable to the smooth kernel function. Instead, the contact force of the discontinuous particles that are in contact with the target particle is used.

[0100] In some embodiments, a method is provided for searching for broken particles in contact with unbroken particles. Specifically, contact occurs when the distance between the unbroken and broken particles is less than a critical distance. Optionally, the critical distance can be defined as twice the radius of the smooth kernel. Optionally, using a linked list method for the search helps improve search efficiency.

[0101] Figure 2 Show Figure 1 Enlarged view of the middle section structure. (Reference) Figure 2 The contact force is the normal contact force F i Cn and tangential contact force F i Cs composition.

[0102] In some embodiments, the normal contact force F i Cn The calculation formula is:

[0103] F i Cn =k n U n n (19)

[0104] U n =R i +R d -d id (20)

[0105]

[0106] Where, k n For normal stiffness, U nR is the normal overlap, where n is the unit outward normal vector at the contact point; i R d Let R be the radius of the two particles. i The initial interparticle spacing of IP is set to ΔP / 2, and the radius of DP particles is typically the radius of a smooth nucleus; d id E represents the distance between the IP and DP particles. i and E d These are the elastic moduli of IP particles and DP particles, respectively.

[0107] In some embodiments, the formula for calculating the tangential contact force involves:

[0108] ΔU s =Δv i -Δv i ·n (22)

[0109] F i Cs ={F i Cs} update +k s ΔU s (twenty three)

[0110]

[0111] In the formula, k s For tangential stiffness, λ i ΔU is the ratio of tangential to normal stiffness. s Δv is the shear increment, n is the unit outward normal vector at the contact point; i The velocity difference between the two IP and DP particles; d id The distance between the IP and DP particles; {F i Cs} update This represents the tangential contact force from the previous step. When IP and DP particles are in simulated contact, due to the adherence to Coulomb's law of friction, their final tangential contact force... The comparison between static friction and sliding friction is shown in the following equation:

[0112]

[0113] Wherein, the static friction force is F i Cs The sliding friction force is When the object is at rest, i.e., ||F i Cs ||<μ||F i CnLet the static friction be the tangential contact force at this point; when the object slides, the maximum static friction is greater than the sliding friction, i.e., ||F. i Cs ||≥μ||F i Cn ||, take the sliding friction force as the tangential contact force at this time.

[0114] Therefore, when there are broken particles in the compactly supported domain of unbroken particles, the parameter governing equation for those particles changes from equation (18) to:

[0115]

[0116] In the formula, f i α It's an external force; here we can take gravity as the force. The value of can be calculated according to formula (25).

[0117] When the target particle is a broken particle, once the fracture surface is transformed into a contact surface, the broken particle is no longer subject to the smooth kernel function. Instead, the contact force of the particles in contact with the target particle takes its place.

[0118] In some embodiments, the method described above for searching for damaged particles in contact with non-damaged particles can be used to search for contacting particles of damaged particles, which will not be repeated here. The calculation of the contact force follows Newton's second law, and can be referred to equations (19) to (25), which will not be repeated here.

[0119] Therefore, the parameter governing equation for the damaged particles changes from equation (18) to:

[0120]

[0121] Finally, based on the damage status of the particles in the model, by selecting an appropriate equation (1), equation (26), or equation (27), the density change rate, acceleration, and velocity of each particle in the model can be calculated. Combined with a preset time step, the density change rate, acceleration, and velocity of a single time step can be integrated over time to obtain the particle density, velocity, and displacement of each particle at each time step.

[0122] In some embodiments, a frog-leapfrog algorithm is used for time integration, recursively obtaining the position, velocity, and density of each base particle at each step. An exemplary integration formula is shown below:

[0123]

[0124] In the formula: t and t0 are the calculation time and the initial time, respectively; Δt is the calculation time step; ρ is the particle density; v is the particle velocity; and x is the particle position coordinates. It should be noted that the density, velocity, and position coordinates corresponding to t0 can be values ​​obtained from the previous calculation or values ​​assigned by the model, while the density, velocity, and position coordinates corresponding to t are obtained from the current calculation.

[0125] The frog-leap algorithm, which uses time integration to solve for the stress, velocity, and displacement changes of materials, is more efficient than the discrete element method and discontinuous deformation analysis method for simulating slope sliding.

[0126] Below, we will describe in detail the simulation method for the failure process of rock slopes provided by the exemplary embodiments of this disclosure, based on the principle of the rock slope fracture contact frame and the above-mentioned parameter control equations applicable to different particles.

[0127] It should be noted that different rock slopes typically exhibit different failure processes, and the characteristic parameters of a rock slope directly determine its failure process. Here, characteristic parameters are parameters used to characterize the properties of a material or phenomenon. Specifically for rock slopes, characteristic parameters can include geometric parameters, structural parameters, matrix particle physical and mechanical parameters, and structural particle physical and mechanical parameters. Geometric parameters describe the shape and size of the rock slope, such as length, width, height, and dip angle. Structural parameters describe structural surfaces, including strike, dip angle, and dip direction. It should be noted that structural surfaces refer to geological interfaces or zones with a certain extension direction and length, and relatively small thickness, formed within the rock mass during geological history; examples include bedding planes, unconformities, joints, faults, foliation surfaces, and cleavage. Physical and mechanical parameters describe the physical and mechanical characteristics of particles, such as density, elastic modulus, Poisson's ratio, cohesion, and angle of internal friction. It should be noted that particles located at structural planes (i.e., structural particles) have different physical and mechanical parameters than matrix particles (i.e., particles not located at structural planes). Those skilled in the art will understand that different rock slopes can have different structural planes, such as bedding planes or joint planes; the same rock slope can also have multiple structural planes, such as simultaneously having bedding planes and joint planes. This disclosure does not specifically limit the structural planes of rock slopes.

[0128] Therefore, in the simulation method for the failure process of rock slopes provided in this disclosure, before constructing the slope model, the characteristic data of the target slope are first obtained. Here, the target slope refers to the slope to be simulated; the characteristic data refers to the characteristic parameter data.

[0129] Next, an initial model is constructed based on the geometric parameters of the target slope. This initial model has a certain height, length, and slope, and can include multiple matrix particles within its region, with preset parameters such as the initial spacing between the matrix particles and the radius of the smoothing kernel. Figure 3 This diagram illustrates a DSPH model of a bedding rock slope according to an embodiment of the present disclosure. Figure 3 It can be seen that the slope height is 32m, the length is 48m, and the slope is 70°.

[0130] Then, based on the structural parameters of the target slope, at least one matrix particle is determined as the structural particle to facilitate the simulation of the slope's bedding, joints, and other structures.

[0131] The method of simulating structural features is illustrated as follows: Based on the structural parameters of the target slope, a preset search algorithm is used to determine the structural particles from at least some of the matrix particles. For example, if the structure is bedding, then the structural particles are bedding particles; if the structure is joint, then the structural particles are joint particles.

[0132] Taking the region search algorithm as an example, the steps for identifying some matrix particles as structural particles are explained in detail. In some embodiments, based on the structural parameters of the target slope, structural feature surfaces can be formed within the region of the initial model. For example, the structural parameters of the target slope include bedding planes inclined at 45°, with a distance of approximately 3m between bedding planes. Based on this, multiple layers corresponding to the structural parameters (such as...) are drawn within the initial model. Figure 3 (As shown). Here, S can be used in multiple layers. i (i = 1, 2, ..., n). Next, search particles are introduced and set on the structural feature surfaces. It should be noted that the number of search particles can be set as needed; for example, one or more search particles can be set for each structural feature surface. When multiple search particles are set for each structural feature surface, the search particles are evenly distributed on the structural feature surface. Here, the search particles can be M. i (i = 1, 2, ..., q) represents the search. It should be understood that each search particle has a pre-defined search radius R. s This allows for accurate positioning of the matrix particles to be transformed. Search radius R s The resulting search domain can be represented by Ω. i Indicated. Optionally, the search radius R s It can be a smooth core with a radius of h. Finally, the search particle M is moved along the structural feature surface. i Search domain Ω i Covering any matrix particle determines that the matrix particle as a structural particle. This sampling method enables rapid simulation of the structural surfaces of a target slope, allowing the construction of rock slope models with structural surfaces.

[0133] Next, corresponding physical and mechanical parameters are assigned to the matrix particles and structural particles, and corresponding boundary conditions are also assigned to the slope model. It should be noted that the boundary conditions can be pre-set, for example... Figure 3 As shown, the bottom boundary is fixed in all directions by the bottom three layers of real particles. Figure 3 (Middle triangle), with three layers of real particles fixed horizontally on each of the left and right boundaries. Figure 3 (Medium-circular shape). By setting boundary conditions, unexpected failure of rock slopes can be effectively avoided, ensuring that the failure process closely resembles the actual failure process. Those skilled in the art can set the boundary conditions reasonably, and this disclosure does not limit this.

[0134] In some embodiments, the process of assigning corresponding physical and mechanical parameters can be carried out in steps. For example, firstly, based on physical and mechanical characteristic data, elastic modulus, Poisson's ratio, density, and gravity are applied to restore the internal stress characteristics of the slope; as the in-situ stress gradually reaches equilibrium, the matrix particles and layered particles are given actual strengths, including but not limited to cohesion, internal friction angle, tensile strength, etc.

[0135] Then, the accumulated displacement generated during the stress equilibrium stage is cleared, and the calculation time is set to 0 seconds, completing the slope model construction. Here, the slope model can be considered as a second model.

[0136] Next, using the slope model and the previously constructed parameter control equations and integral equations, the failure process of the rock slope can be simulated. It should be noted that at any time step, the target parameter values ​​for each particle, such as velocity, position, and density, are calculated separately.

[0137] Figure 9 This diagram illustrates a partial flow chart of a simulation method for the failure process of a rock slope according to an embodiment of this disclosure. Specifically, Figure 9 The process of the k-th calculation is shown. Here, taking the k-th loop as an example, the calculation steps for the target particle are explained as follows: k is an integer greater than 1 and less than or equal to the preset number of calculations.

[0138] In the initial calculation, both matrix particles and structural particles are in an undamaged state in their historical damage state, so step 901 can be omitted. The current damage state 902 of the target particle is determined according to equations (16) and (17). If the target particle is undamaged, the target parameter change rate 9083 is calculated according to equation (1), and the target parameter value 909 is obtained by integration according to equation (28). If the target particle is damaged, the series of steps marked by dashed box A are executed, and the target parameter change rate 9081 is calculated according to equation (27). The target parameter value 909 is obtained by integration according to equation (28). Here, in each calculation process, multiple particles of the slope model are selected as target particles one by one. After all particles have been calculated, the target parameters of each particle are updated. The current damage state is distinguished from the historical damage state. The current damage state is the damage state obtained in this calculation; the historical damage state refers to the damage state determined from the first calculation to the previous calculation.

[0139] Next, the damage state and target parameters of multiple particles are calculated in a loop.

[0140] Determine the historical damage state of the target particle 901; It should be noted that the particle will not recover after it is damaged. Therefore, if any calculation determines that the particle is destroyed, the particle becomes a damaged particle. The damage state of the particle is recorded by marking the damaged particles. Step 901 can directly query the mark to determine the historical damage state of the target particle.

[0141] If the target particle is damaged, the series of steps marked by the dashed box A are executed, including: determining at least one contacting particle 903 of the target particle according to the preset contact distance. Here, the preset contact distance can be referred to the aforementioned contact distance, and will not be repeated here; next, determining the relative motion between the two. It should be noted that the relative motion can be determined by calculating the change in distance between the two. For example, if the distance between the two at k-1 times is different from the distance at k-2 times, it indicates that there is relative motion. The contact force 905 is calculated based on the sliding friction force. Otherwise, there is no relative motion. The contact force 905 is calculated based on the static friction force; then the target parameter change rate 9081 is calculated using formula (27).

[0142] If the history of the target particle is intact, then the current damaged state 902 of the target particle is determined. Here, as mentioned before, the method for determining the damaged state can be calculated using the Mohr-Coulomb theory, which will not be elaborated further.

[0143] If it is determined that the target particle is currently damaged, then execute dashed box A and step 9081.

[0144] If it is determined that the current state of the target particle is undamaged, then the nearest neighbor particle is determined, and its historical damage state 906 is determined; here, the method for determining the nearest neighbor particle is as described above, and will not be repeated.

[0145] If the nearest neighboring particle has a broken particle, then execute step 907 of dashed box A, and then use equation (26) to calculate the target parameter change rate 9082.

[0146] If there are no damaged particles in the nearest neighboring particles, the target parameter change rate 9083 is calculated using equation (1).

[0147] Finally, based on the target parameter change rate and preset time step obtained from 9081, 9082, and 9083, the target parameter is calculated using equation (28).

[0148] The calculation is repeated multiple times until a preset condition is met, at which point the calculation can stop. It should be noted that the preset condition can be a preset number of calculations, or the target parameter values ​​of some particles meeting a preset rule. Here, the preset rule could be that the particle's position changes by a preset distance relative to time 0, such as 5 meters.

[0149] The simulation method for the failure process of rock slopes provided in this disclosure will be described below with reference to specific embodiments.

[0150] Example 1: Bean-parallel rock slope

[0151] Based on the characteristic parameters of bedding rock slopes, construct as follows: Figure 3 The model shown is of a slope. This model contains a set of bedding planes inclined at 45° to the slope surface, with a spacing of approximately 3m between the bedding planes. The slope height is 32m, the length is 48m, and the slope angle is 70°. The initial model consists of a total of 32,000 DSPH particles with an initial particle spacing of dx = 0.2m and a computation time step of Δt = 1.5 × 10⁻⁶. -5 Using a region search method, 27,400 rock matrix particles and 4,600 rock bedding particles were identified.

[0152] Table 1 shows some of the physical and mechanical properties of the rock matrix particles and layering particles. Based on the parameters in Table 1, the particles are assigned corresponding parameters. Specific steps may include: stress equilibrium, determining location, and setting a zero time point. Simultaneously, corresponding boundary conditions are assigned to the layered slope model. The bottom boundary is fixed in all directions by the bottom three layers of solid particles, while the left and right boundaries are each fixed in the horizontal direction by three layers of solid particles.

[0153] Table 1 DSPH parameters of rock slopes

[0154]

[0155]

[0156] To accelerate slope instability, gravitational acceleration is continuously increased to induce rock slope fracture. The artificial viscosity control parameters for the rock mass are β1 = 4 and β2 = 1. The radius of the smooth core is 0.25 m. The maximum time step is set to 40,000 steps.

[0157] Next, the target parameters are calculated for each time step of the model particles until the preset time step is reached, which is 40,000 steps in this case, at which point the slope calculation terminates. The specific calculation method can be found above and will not be repeated here.

[0158] Based on the calculation results from the above steps, a drawing can be performed. Figures 4(a) to 4(c) Figure 4(a) shows a schematic diagram of the failure process of a bedding rock slope provided in an embodiment of this disclosure. Figure 4(a) illustrates the crack initiation and propagation process of the bedding slope. Figure 4(b) shows a schematic diagram of displacement deformation corresponding to Figure 4(a); Figure 4(c) shows a schematic diagram of the evolution of the maximum principal stress corresponding to Figure 4(a).

[0159] Figure 4(d) shows a schematic diagram of the experimental results of a bedding rock slope provided in an embodiment of this disclosure. The experiment here is a physical simulation experiment. Figure 4(e) shows a schematic diagram of the FDEM simulation results of a bedding rock slope provided in an embodiment of this disclosure. Comparison Figures 4(a) to 4(e) It can be seen that the numerical simulation method provided in the embodiments of this disclosure is effective and feasible.

[0160] Example 2: Anti-dip layered rock slope

[0161] Based on the characteristic parameters of the anti-dip rock slope, construct as follows: Figure 5 The slope model shown. Figure 5 This diagram illustrates a DSPH model of a reverse-dipping bedding rock slope according to an embodiment of this disclosure. The model includes a set of bedding planes inclined inwards at 70°, with a spacing of approximately 2m between the bedding planes. The slope height is 32m, the length is 48m, and the slope is 70°. The initial model consists of a total of 32,000 DSPH particles, with an initial particle spacing of dx = 0.2m and a calculation time step of Δt = 1.5 × 10⁻⁵ s. Using a region search method, 30,000 rock matrix particles and 2,000 rock bedding particles were identified.

[0162] Table 2 shows some of the physical and mechanical properties of the rock matrix particles and layering particles. Based on the parameters in Table 2, the particles are assigned corresponding parameters. Specific steps may include: stress equilibrium, determining location, and setting a zero time point. Simultaneously, corresponding boundary conditions are assigned to the layered slope model. The bottom boundary is fixed in all directions by the bottom three layers of solid particles, while the left and right boundaries are each fixed in the horizontal direction by three layers of solid particles.

[0163] Table 2 DSPH parameters for rock slopes

[0164]

[0165]

[0166] To accelerate slope instability, gravitational acceleration is continuously increased to induce rock slope fracture. The artificial viscosity control parameters for the rock mass are β1 = 4 and β2 = 1. The radius of the smooth core is 0.25 m. The maximum time step is set to 80,000 steps.

[0167] Next, the target parameters are calculated for each time step of the model particles until the preset time step is reached, which is 80,000 steps in this case, at which point the slope calculation terminates. The specific calculation method can be found above and will not be repeated here.

[0168] Based on the calculation results from the above steps, a drawing can be performed. Figures 6(a) to 6(c) Figure 6(a) shows a schematic diagram of the failure process of a reverse-dip rock slope provided in an embodiment of this disclosure. Figure 6(a) illustrates the crack initiation and propagation process of the reverse-dip slope. Figure 6(b) shows a schematic diagram of displacement deformation corresponding to Figure 6(a); Figure 6(c) shows a schematic diagram of the evolution of the maximum principal stress corresponding to Figure 6(a).

[0169] Figure 6(d) shows a schematic diagram of the experimental results of a reverse-dip rock slope provided in an embodiment of this disclosure. The experiment here is a physical simulation experiment. Figure 6(e) shows a schematic diagram of the simulation results of a reverse-dip rock slope using the general discrete element method provided in an embodiment of this disclosure. Comparison Figures 6(a) to 6(e) It can be seen that the numerical simulation method provided in the embodiments of this disclosure is effective and feasible.

[0170] Example 3: Discontinuous Jointed Rock Slope

[0171] Based on the characteristic parameters of discontinuous jointed rock slopes, construct as follows: Figure 7 The slope model shown. Figure 7 This diagram illustrates a DSPH model of a discontinuous jointed rock slope according to an embodiment of this disclosure. The model includes a set of joint surfaces inclined at approximately 45° to the slope face. The joints are filled with DP particles, representing pre-existing joints. The joint length is 4m, the parallel spacing between adjacent joints is 3m, the joint discontinuity is 66%, and the salt bridge length is 2m. The slope height is 32m, the length is 48m, and the slope is 70°. The initial model consists of a total of 32,000 DSPH particles with an initial particle spacing of dx = 0.2m and a calculation time step of Δt = 1.5 × 10⁻⁵ s. Using a region search method, 31,000 rock matrix particles and 1,000 rock joint particles were identified.

[0172] Table 3 shows some of the physical and mechanical properties of the rock matrix particles and layering particles. Based on the parameters in Table 3, the particles are assigned corresponding parameters. Specific steps may include: stress equilibrium, determining location, and setting a zero time point. Simultaneously, corresponding boundary conditions are assigned to the layered slope model. The bottom boundary is fixed in all directions by the bottom three layers of solid particles, while the left and right boundaries are each fixed in the horizontal direction by three layers of solid particles.

[0173] Table 3 DSPH parameters for rock slopes

[0174]

[0175] To accelerate slope instability, gravitational acceleration is continuously increased to induce rock slope fracture. The artificial viscosity control parameters for the rock mass are β1 = 4 and β2 = 1. The radius of the smooth core is 0.25 m. The maximum time step is set to 60,000 steps.

[0176] Next, the target parameters are calculated for each time step of the model particles until the preset time step is reached, which is 60,000 steps in this case, at which point the slope calculation terminates. The specific calculation method can be found above and will not be repeated here.

[0177] Based on the calculation results from the above steps, a drawing can be performed. Figures 8(a) to 8(c) Figure 8(a) shows a schematic diagram of the failure process of a discontinuous jointed rock slope provided in an embodiment of this disclosure. Figure 8(a) illustrates the crack initiation and propagation process of a reverse-dip slope. Figure 8(b) shows a schematic diagram of displacement deformation corresponding to Figure 8(a); Figure 8(c) shows a schematic diagram of the evolution of the maximum principal stress corresponding to Figure 8(a).

[0178] It should be noted that the normal stiffness in each embodiment can be calculated according to equation (21).

[0179] This exemplary embodiment provides a method for simulating the failure process of a rock slope, including:

[0180] Construct a slope model of the target slope; wherein the slope model includes multiple particles, each assigned different parameter values, and the multiple particles belong to at least two types of particles; here, the at least two types of particles may include matrix particles, bedding particles, joint particles, etc.; the parameters of the particles may be physical and mechanical parameters, including but not limited to density, elastic modulus, Poisson's ratio, cohesion, tensile strength, internal friction angle, normal stiffness, friction coefficient, position, velocity, smooth core radius, etc.

[0181] Based on a preset time step and the parameter values, the damage state and target parameter values ​​of the plurality of particles are calculated iteratively until a predetermined condition is met; here, the target parameters may include at least one of velocity, acceleration and density; optionally, the predetermined condition is selected from a predetermined number of times or the target parameter values ​​of some of the particles satisfying a preset rule, the preset rule can be referred to the foregoing and will not be repeated here;

[0182] Output at least once the damage state of the plurality of particles and the target parameter value; wherein,

[0183] The calculation of the damage state and target parameter values ​​of the plurality of particles includes:

[0184] Determine the target particle from among the plurality of particles;

[0185] In response to determining that the target particle's historical damage state is that it is an undamaged particle, the nearest neighbor particle of the target particle is determined based on the previous parameter value and a preset compact support region.

[0186] Based on the Mohr-Coulomb theory, the damage state of the target particle is determined; specifically, refer to the foregoing, which will not be repeated here.

[0187] Based on the historical parameter values, the current damage state of the target particle, and the historical damage state of the nearest neighbor particle, the corresponding preset target parameter control equation is determined and the target parameter change rate of the current target particle is calculated; here, the preset target parameter control equation can be Equation (1), Equation (26), or Equation (27);

[0188] Based on a preset time step, the previous parameter value, and the target parameter change rate, the current target parameter value of the target particle is calculated; wherein, the current target parameter value is used to update the previous parameter value.

[0189] Using the above method, the entire process of crack initiation and propagation (corresponding to the generation of broken particles), contact formation, frictional slip (relative motion between particles) and complete failure of rock slopes can be accurately captured through iterative calculation. Furthermore, the transformation from continuous domain to discontinuous domain is achieved within a single framework, which can significantly improve computational efficiency and is applicable to rock slope engineering problems of different scales.

[0190] In some embodiments, constructing the slope model of the target slope includes:

[0191] A first model is constructed based on the first feature data of the target slope; wherein the first model includes first particles; the first feature data may include geometric feature data; and the first model may be the aforementioned initial model.

[0192] Based on the second feature data of the target slope, a preset search algorithm is used to determine at least a portion of the first particles as second particles; here, the second feature data may be structural feature data.

[0193] By introducing structural particles into the slope model in this way, the actual target slope can be more accurately reproduced.

[0194] In some embodiments, the step of acquiring and determining at least a portion of the first particles as second particles using a preset search algorithm based on the second feature data of the target slope includes:

[0195] Based on the second feature data, a structural feature surface is formed within the region of the first model;

[0196] A search particle is introduced and set on the structural feature surface;

[0197] The search particle is moved along the structural feature surface, and in response to the search domain formed by the search radius of the search particle covering any of the first particle, the first particle is determined as the second particle.

[0198] By forming structural feature surfaces, and then setting search particles on these surfaces, it is helpful to improve search efficiency and reduce computational load.

[0199] In some embodiments, determining the corresponding preset target parameter control equation based on historical parameter values, the current damage state of the target particle, and the historical damage state of the nearest neighbor particle includes:

[0200] In response to the determination that the current target particle and each of the nearest neighbor particles are non-destructive particles, the preset target parameter control equation is a discontinuous smooth particle hydrodynamic control equation, as shown in equation (1); here, the interaction between non-destructive particles is transmitted by constructing virtual bonds;

[0201] In response to determining that the current target particle is a non-destructive particle and at least one of the nearest neighbor particles is a destructive particle, the preset target parameter control equation is a continuous particle control equation, as shown in equation (26).

[0202] In response to determining that the current target particle is a broken particle, the preset target parameter control equation is a discontinuous particle control equation, as shown in equation (27); here, by destroying the above virtual bond through the pseudo-spring method based on the Mohr-Coulomb and maximum tensile stress failure criteria, the technical effect of simulating cracks can be achieved.

[0203] Optionally, the continuous particle control equation and the discontinuous particle control equation are constructed based on the discontinuous smooth particle hydrodynamic equation and the contact force between particles.

[0204] In this way, a new contact force is formed between the broken particles and the contacting particles. The contact force can be used to simulate frictional slippage and separation / contraction between fracture surfaces.

[0205] In some embodiments, calculating the damage state and target parameters of the plurality of particles further includes:

[0206] In response to determining that the historical damage state of the target particle is a damaged particle, at least one contacting particle of the target particle is determined according to a preset critical distance; here, the preset critical distance can be twice the smooth core radius;

[0207] Based on the historical parameter values ​​of the target particle and the contact particle, the rate of change of the target parameter of the current target particle is calculated using the discontinuous particle control equation; here, the discontinuous particle control equation can be as shown in equation (27);

[0208] Based on a preset time step, the previous parameter value, and the target parameter change rate, the current target parameter value of the target particle is calculated; wherein, the current target parameter value is used to update the previous parameter value. It should be noted that historical parameter values ​​include the previous parameter value; the current target parameter value is the value obtainable in this calculation.

[0209] In some embodiments, the discontinuous particle control equation includes contact force; the step of calculating the current rate of change of the target parameter of the target particle using the discontinuous particle control equation based on the historical parameter values ​​of the target particle and the contact particle includes:

[0210] Based on the historical parameter values ​​of the target particle and the contact particle, determine the current change in distance between the target particle and the contact particle;

[0211] In response to determining that the distance change is zero, the static friction between the target particle and the contacting particle is used to determine the contact force of the target particle;

[0212] In response to determining that the change in distance is greater than zero, the sliding friction between the target particle and the contacting particle is used to determine the contact force of the target particle.

[0213] By adopting this approach, frictional slip is considered in the discontinuous particle control equations, making the calculation results more consistent with the actual slope aging process.

[0214] It should be noted that the method of this disclosure embodiment can be executed by a single device, such as a computer or server. The method of this embodiment can also be applied to a distributed scenario, where multiple devices cooperate to complete the task. In such a distributed scenario, one of these devices may execute only one or more steps of the method of this disclosure embodiment, and the multiple devices will interact with each other to complete the method described.

[0215] It should be noted that the above description describes some embodiments of this disclosure. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in a different order than that shown in the above embodiments and still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0216] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this disclosure also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the simulation method described in any of the above embodiments.

[0217] Figure 10 This embodiment illustrates a more specific hardware structure of an electronic device, which may include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040, and a bus 1050. The processor 1010, memory 1020, input / output interface 1030, and communication interface 1040 are interconnected internally via the bus 1050.

[0218] The processor 1010 can be implemented using a general-purpose CPU (Central Processing Unit), microprocessor, application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.

[0219] The memory 1020 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 can store the operating system and other applications. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program code is stored in the memory 1020 and is called and executed by the processor 1010.

[0220] The input / output interface 1030 is used to connect input / output modules to realize information input and output. Input / output modules can be configured as components within the device (not shown in the figure) or externally connected to the device to provide corresponding functions. Input devices may include keyboards, mice, touchscreens, microphones, various sensors, etc., while output devices may include displays, speakers, vibrators, indicator lights, etc.

[0221] The communication interface 1040 is used to connect a communication module (not shown in the figure) to enable communication between this device and other devices. The communication module can communicate via wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.).

[0222] Bus 1050 includes a pathway for transmitting information between various components of the device, such as processor 1010, memory 1020, input / output interface 1030, and communication interface 1040.

[0223] It should be noted that although the above-described device only shows the processor 1010, memory 1020, input / output interface 1030, communication interface 1040, and bus 1050, in specific implementations, the device may also include other components necessary for normal operation. Furthermore, those skilled in the art will understand that the above-described device may only include the components necessary for implementing the embodiments of this specification, and not necessarily all the components shown in the figures.

[0224] The electronic devices described above are used to implement the corresponding simulation methods in any of the foregoing embodiments and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0225] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this disclosure also provides a non-transitory computer-readable storage medium that stores computer instructions for causing the computer to perform the simulation method as described in any of the above embodiments.

[0226] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.

[0227] The computer instructions stored in the storage medium of the above embodiments are used to cause the computer to execute the simulation method as described in any of the above embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0228] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of this disclosure (including the claims) is limited to these examples; within the framework of this disclosure, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of the embodiments of this disclosure as described above, which are not provided in detail for the sake of brevity.

[0229] Although this disclosure has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed.

[0230] This disclosure is intended to cover all such substitutions, modifications, and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this disclosure should be included within the scope of protection of this disclosure.

Claims

1. A method for simulating the failure process of a rock slope, characterized in that, include: Construct a slope model of the target slope; wherein the slope model includes multiple particles that are assigned different parameter values, and the multiple particles belong to at least two types of particles; Based on the preset time step and the parameter values, the damage state and target parameter values ​​of the multiple particles are calculated iteratively until the predetermined conditions are met. Output at least once the damage state of the plurality of particles and the target parameter value; wherein, The calculation of the damage state and target parameter values ​​of the plurality of particles includes: Determine the target particle from among the plurality of particles; In response to determining that the target particle's historical damage state is that it is an undamaged particle, the nearest neighbor particle of the target particle is determined based on the previous parameter value and a preset compact support region. Based on the Mohr-Coulomb theory, the damage state of the target particle is determined. Based on historical parameter values, the current damage state of the target particle, and the historical damage state of the nearest neighboring particle, determine the corresponding preset target parameter control equation and calculate the target parameter change rate of the current target particle; Based on a preset time step, the previous parameter value, and the target parameter change rate, the current target parameter value of the target particle is calculated; wherein, the current target parameter value is used to update the previous parameter value; wherein, The step of determining the corresponding preset target parameter control equation based on historical parameter values, the current damage state of the target particle, and the historical damage state of the nearest neighbor particle includes: In response to determining that the current target particle and each of the nearest neighbor particles are non-destructive particles, the preset target parameter control equation is a discontinuous smooth particle hydrodynamic control equation. In response to determining that the current target particle is a non-destructive particle and at least one of the nearest neighbor particles is a destructive particle, the preset target parameter control equation is a continuous particle control equation; the continuous particle control equation is shown in equation (26); Equation (26); In response to determining that the current target particle is a broken particle, the preset target parameter control equation is a discontinuous particle control equation; the discontinuous particle control equation is shown in equation (27); Equation (27); in, It is particle density; m It is the particle mass; It is the index of the target particle; j It is the index of the nearest neighbor particle; v It is the particle velocity; x These are the particle's position coordinates; t It is a calculation of time; U Represents the set of undamaged particles; α and β Indicates that they can be taken in Cartesian coordinate systems respectively. x , y and z The summation convention, also known as Einstein's summation convention, is applied to duplicate indexes. W ij For smooth kernel functions; Represented as the Cauchy stress tensor; It is the Dirac function; Artificial viscosity; Artificial stress; It is gravity; k n For normal stiffness, U n The normal overlap is... n It is the unit outward normal vector at the point of contact; Calculated based on formula (25); Equation (25); where, It is static friction; It is sliding friction; Normal contact force; μ The coefficient of friction; The continuous particle control equation and the discontinuous particle control equation are constructed based on the discontinuous smooth particle hydrodynamic equation and the contact force between particles.

2. The simulation method according to claim 1, characterized in that, The slope model for constructing the target slope includes: A first model is constructed based on the first feature data of the target slope; wherein the first model includes a first particle; Based on the second feature data of the target slope, a preset search algorithm is used to determine at least a portion of the first particles as second particles.

3. The simulation method according to claim 2, characterized in that, The step of acquiring and determining at least a portion of the first particles as second particles based on the second feature data of the target slope using a preset search algorithm includes: Based on the second feature data, a structural feature surface is formed within the region of the first model; A search particle is introduced and set on the structural feature surface; The search particle is moved along the structural feature surface, and in response to the search domain formed by the search radius of the search particle covering any of the first particle, the first particle is determined as the second particle.

4. The simulation method according to claim 1, characterized in that, The calculation of the damage state and target parameters of the plurality of particles further includes: In response to determining that the historical damage state of the target particle is a damaged particle, at least one contact particle of the target particle is determined according to a preset critical distance. Based on the historical parameter values ​​of the target particle and the contacting particle, the rate of change of the target parameter of the current target particle is calculated using the discontinuous particle control equation. Based on a preset time step, the previous parameter value, and the target parameter change rate, the current target parameter value of the target particle is calculated; wherein, the current target parameter value is used to update the previous parameter value.

5. The simulation method according to claim 4, characterized in that, The discontinuous particle control equation includes contact force; the step of calculating the current rate of change of the target parameter of the target particle using the discontinuous particle control equation based on the historical parameter values ​​of the target particle and the contact particle includes: Based on the historical parameter values ​​of the target particle and the contact particle, determine the current change in distance between the target particle and the contact particle; In response to determining that the distance change is zero, the static friction between the target particle and the contacting particle is used to determine the contact force of the target particle; In response to determining that the change in distance is greater than zero, the sliding friction between the target particle and the contacting particle is used to determine the contact force of the target particle.

6. The simulation method according to claim 1, characterized in that, The parameters of the particles include density, elastic modulus, Poisson's ratio, cohesion, tensile strength, internal friction angle, normal stiffness, coefficient of friction, position, velocity, and smooth core radius; and / or The target parameter includes at least one of velocity, acceleration, and density.

7. The simulation method according to claim 1, characterized in that, The predetermined conditions are selected from a predetermined number of times or the target parameter values ​​of some of the particles satisfying a preset rule.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable by the processor, characterized in that, The processor implements the simulation method according to any one of claims 1 to 7 when executing the computer program.

9. A non-transitory computer-readable storage medium, characterized in that, The non-transitory computer-readable storage medium stores computer instructions for causing the computer to perform the simulation method according to any one of claims 1 to 7.

Images