Construction method of quasi-regular double six-interval triangular network for searching potential sliding surface of slope

Abstract

This invention provides a quasi-regular double-six discontinuous triangular mesh construction method for searching potential slope sliding surfaces, comprising: iteratively generating a quasi-regular initial triangular mesh adapted to the slope model boundary; locally dividing into six meshes to generate a six-directional arrayed topological mesh; constructing a data format of six-node triangular elements and high-order discontinuities; and solving the potential slope sliding surface using a second-order cone programming model. This invention adapts to the geometric boundaries of discrete slope models and meets the special requirements of the quasi-regular double-six discontinuous triangular mesh, forming a six-directional continuous discontinuous mesh structure within any small region. Furthermore, the double-six discontinuous mesh can be further pre-defined with local partitions, reducing resource consumption of mesh elements in non-critical areas. Applied to the upper limit finite element method of limit analysis, it achieves high-precision searching of potential slope sliding surfaces through rapid solving of the second-order cone programming problem, thus reducing computational burden while maintaining search accuracy.

Description

Technical Field

[0001] This invention relates to the field of numerical analysis technology in geotechnical engineering, and in particular to a method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes. Background Technology

[0002] Upper limit finite element analysis (UPEM) within the framework of plastic limit analysis is a powerful tool for conducting slope stability analysis, searching for potential sliding surfaces, and identifying failure modes. The UEM for slope stability relies on a triangular mesh formed by discretizing the computational model, and the structural characteristics of this mesh significantly affect the computational accuracy of the stability analysis.

[0003] Unlike the continuous medium elastoplastic finite element method, the upper limit finite element triangulation (TLE) allows for velocity discontinuities on all common edges of elements to accommodate discontinuous deformation and failure phenomena such as shear bands. However, the effectiveness of velocity discontinuities in the TLE is limited by their location and orientation. For a specific slope stability analysis model, whether shear failure will occur at any point (or small region) in the TLE and its dominant slip direction cannot be known in advance, making it impossible to pre-set optimal discontinuity orientations. Therefore, even with global velocity discontinuities, the actual computational efficiency of a conventional TLE remains low, offering limited improvement to the accuracy of slope stability analysis. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of the aforementioned background technology by providing a method for constructing a "double-six" discontinuous triangular network with six-node triangular elements and a six-directional topology. This method applies the quasi-regular "double-six" discontinuous triangular network to the upper limit finite element method of limit analysis, and achieves high-precision search for potential slope sliding surfaces through the rapid solution of the second-order cone programming problem.

[0005] To achieve the above objectives, this invention provides a quasi-regular double-six discontinuous triangulation network construction method for searching potential sliding surfaces of slopes, comprising the following steps:

[0006] S1. Construct the initial point set of the slope model according to the preset unit size, perform Delaunay triangulation to generate the initial mesh, compare the triangulation mesh with the truss structure, apply Hooke's law to calculate the resultant force vector of all truss elastic members, adjust the node coordinates along the direction of the resultant force vector, and iteratively implement node coordinate optimization and mesh reconstruction until an approximate equilateral triangular mesh that adapts to the boundary of the slope model and is arranged in a quasi-regular array is generated.

[0007] S2, a local densification strategy driven by densification edges is applied to the key areas of the approximately equilateral triangular mesh to form a six-directional triangular mesh with a local array arrangement. In the key areas of the slope model, a one-to-six densification is applied. In the transition area, a passive densification strategy is applied based on the number of densification edges. The original units are retained in the non-densification area.

[0008] S3, based on the six-direction triangulation, extracts the coordinate matrix and element information matrix, generates a six-node triangular element with fifteen degrees of freedom by interpolating the midpoint of the edge, constructs a high-order discontinuity line on the common edge, adds non-negative auxiliary variables and constraints to realize the velocity jump between elements, and forms the double six-discontinuity triangulation data format required for the upper limit finite element of limit analysis.

[0009] S4. Based on the upper bound theorem of limit analysis, a second-order cone programming model for slope stability is established to search for the critical instability state of plastic flow in the slope. The model outputs the point cloud of dissipated energy density distribution with high-order triangular elements and velocity discontinuities isomorphic, and obtains the potential sliding surface of the slope limit state.

[0010] Furthermore, S1 includes:

[0011] According to the preset initial grid size h 0. Generate an initial large-area regular array node set, filter the internal nodes of the slope model, merge the coordinates of fixed points on the same boundary, and construct the initial point set;

[0012] Perform Delaunay triangulation on the input point set, update the internal mesh structure of the model, sort the vertex index of the cell edge in ascending order, and remove duplicate edges;

[0013] Let the element edges satisfy Hooke's law, and assume the equivalent spring stiffness is... K , computational unit node i The resultant force at the point:

[0014]

[0015] in, For nodes i exist x Component of force in direction, For nodes i exist y Component of force in direction, For nodes i With nodes j The distance between them , They are nodes i , j coordinates For nodes i With nodes j Expected length between:

[0016]

[0017] in, It is the sum of the squares of the lengths of all cell edges in the grid. This is the sum of the squares of the lengths of all cell edges in the desired mesh.

[0018] The node coordinates are iteratively updated based on explicit time integration, causing the nodes to shift along the direction of the resultant force vector;

[0019] For nodes that exceed the boundary after offset, the implicit distance function along the boundary is used. Projecting along the negative gradient direction, we constrain it back to the boundary:

[0020]

[0021] in, g Distance function right x and y The sum of squares of partial derivatives;

[0022] The algorithm iterates through node coordinate offsets and mesh reconstruction, recording the sum of squared offsets of all nodes within the model during each iteration. The process stops when the sum of the squares of the offsets of all nodes in the computational domain is less than a threshold or when the preset number of iterations is met. An approximate equilateral triangular mesh with a quasi-regular array arrangement adapted to the boundary is generated inside the model.

[0023] Furthermore, S2 specifically includes the following sub-steps:

[0024] S21, Traverse all units and select triangular units located within the key area based on the centroid coordinate discrimination criterion;

[0025] S22, extract and mark adjacent pre-encrypted edge data according to the encryption unit index, traverse all units to count the number of pre-encrypted edges contained in them, and mark four types of encryption edge states: no encryption edge unit, one encryption edge unit, two encryption edge units, and three fully encrypted edge units.

[0026] S23. For cells without encrypted edges, the original mesh topology is preserved, and no encryption operation is performed. When a cell contains only one encrypted edge, the midpoint splitting method is used. A new node is inserted at the midpoint of the encrypted edge, and the new node is connected to the vertex of the non-encrypted edge, splitting the original cell into two sub-cells. When a cell contains two encrypted edges, three-way splitting encryption is performed: the longest encrypted edge is identified, its midpoint is taken as the pivot node, and the pivot node is connected to the midpoint and corresponding vertex of the adjacent encrypted edge to generate three sub-cells. When all three edges of a cell need to be encrypted, a one-to-six encryption is performed: the centroid node of the cell is calculated as the core node, and the centroid node is connected to the midpoint and original vertex of each edge in a counterclockwise order, and the original cell is deleted to generate six sub-cells.

[0027] Furthermore, S3 includes the following sub-steps:

[0028] S31, extract the cell coordinate matrix [P] and cell information matrix [Q] as the output of the mesh processing. [Q] stores the information of the three vertices of each cell.

[0029] S32, obtain all element edge numbers, and generate a six-node triangular high-order element with fifteen degrees of freedom by interpolating the edge midpoints;

[0030] S33, all element edges in the model are identified as common edges and boundary edges based on the boundary identifier. High-order velocity discontinuities are constructed on the common edges of the elements. There are jumps in the velocity and plastic multipliers of the nodes on both sides of the discontinuity. The element numbers on both sides of the velocity discontinuity and the positions of the nodes on both sides of the discontinuity are searched, judged and established. Auxiliary variables are introduced on the velocity discontinuity to remove the absolute values ​​in the associated flow rules.

[0031] Furthermore, in the S32 nodal triangle higher-order element, the velocity at any point within the element... Represented as a quadratic function of the nodal coordinates using shape functions:

[0032]

[0033] in, , They are nodes i exist x and y The velocity component in the direction, For unit shape functions;

[0034] Plastic multiplier at any point within the unit Obtained by interpolation of the three vertex values:

[0035]

[0036] in, For nodes i Plastic multipliers at the location.

[0037] Furthermore, the critical instability state of the slope in S4 is obtained through the upper bound theorem of limit analysis. When the slope reaches the limit state, there are infinitely many permissible velocity fields such that the internal dissipated energy is no greater than the work done by the external force.

[0038]

[0039] in, For internal energy dissipation, and These are the plastic stress and strain vectors, respectively. Let be the velocity vector along the discontinuity line. Doing work for external forces and All The function, Let be the yield function. This is the computational domain of the model.

[0040] Furthermore, in S4, the objective function is minimized using the upper bound theorem, forming a second-order cone programming model:

[0041]

[0042] in, For plastic flow constraints on the element and velocity discontinuity line, For boundary velocity constraints, A To satisfy the coefficient matrix corresponding to the associated flow rule, B This is the coefficient matrix corresponding to the boundary constraint equations satisfied by the nodes at the boundary. It is a column vector composed of unknown variables such as the velocity of all unit nodes in the model, auxiliary variables, and plasticity multipliers. Add constraints to the second-order cone programming, parameters , , Let be the angle between a point on the yield surface of the Mohr-Coulomb yield criterion and the Y-axis. For speed constraints related to self-weight, This is the load matrix.

[0043] Furthermore, the internal dissipation energy comprises two parts: the internal dissipation energy of the six-node triangular element and the dissipation energy along the higher-order velocity discontinuity line.

[0044]

[0045] in, For the distribution domain of plastic units, For velocity jumps on discontinuous lines, For velocity discontinuity line sets.

[0046] Furthermore, the dissipation density distribution point cloud is displayed by interpolating the model's dissipation energy density to the element nodes. k No. i Energy dissipation density at each node Represented as:

[0047]

[0048] in, For unit k area, j For unit k In and nodes i Adjacent vertex nodes, c , These represent the cohesion and internal friction angle of the soil, respectively. For the first k Plastic multiplier vectors of each unit , a For constant coefficient terms, For nodes i, j The corresponding auxiliary variables.

[0049] Furthermore, the potential sliding surface of the slope is characterized by a point cloud of normalized dissipative energy density:

[0050]

[0051] in, This represents the maximum energy dissipation density of the model.

[0052] The above-described solution of the present invention has the following beneficial effects:

[0053] The quasi-regular double-six discontinuous triangular mesh construction method for searching potential slope sliding surfaces provided by this invention achieves topology optimization of the mesh by introducing force balance iterative updates, making it adaptable to the geometric boundaries of the discrete slope model and meeting the special requirements of the quasi-regular double-six discontinuous triangular mesh, forming a six-directional continuous discontinuous mesh structure within any small region; at the same time, the double-six discontinuous mesh can further preset local partitions, reducing the resource consumption of non-critical region mesh units. Therefore, applying the quasi-regular double-six discontinuous triangular mesh to the upper limit finite element method of limit analysis and achieving high-precision search of potential slope sliding surfaces through the rapid solution of the second-order cone programming problem can reduce the computational burden while ensuring search accuracy;

[0054] Other beneficial effects of the present invention will be described in detail in the following detailed description section. Attached Figure Description

[0055] Figure 1 This is a flowchart of the steps of the present invention;

[0056] Figure 2 This is a schematic diagram of the triangular unit topology of the present invention, wherein (a) contains one encrypted edge, (b) contains two encrypted edges, and (c) contains three encrypted edges;

[0057] Figure 3 This is a schematic diagram of a local six-directional through-grid of the present invention;

[0058] Figure 4 This is a schematic diagram of a 15-DOF, 6-node triangular element of the present invention;

[0059] Figure 5 This is a schematic diagram of the six-node velocity discontinuity line of the present invention;

[0060] Figure 6This is a schematic diagram of the slope stability analysis model of the present invention;

[0061] Figure 7 The initial mesh of the approximate equilateral triangular mesh is used for the array arrangement of the model boundary in this invention;

[0062] Figure 8 This is a six-directional topological mesh arranged in a local array according to the present invention. Detailed Implementation

[0063] The following specific examples illustrate the implementation of this disclosure. Those skilled in the art can easily understand other advantages and effects of this disclosure from the content disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of this disclosure, and not all of them. This disclosure can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of this disclosure. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. Based on the embodiments in this disclosure, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this disclosure.

[0064] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this disclosure, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0065] It should also be noted that the illustrations provided in the following embodiments are merely schematic representations of the basic concept of this disclosure. The illustrations only show components relevant to this disclosure and are not drawn according to the actual number, shape, and size of components in implementation. In actual implementation, the type, quantity, and proportion of each component can be arbitrarily changed, and the component layout may be more complex. Furthermore, specific details are provided in the following description to facilitate a thorough understanding of the examples. However, those skilled in the art will understand that the described aspects can be practiced without these specific details.

[0066] like Figure 1As shown, embodiments of the present invention provide a quasi-regular double-six-discontinuous triangular mesh construction method for searching potential sliding surfaces of slopes. The ordinary Delaunay triangular mesh is analogized to a truss structure. Iterative calculations of elastic members obeying Hooke's law determine the equilibrium position of the resultant force at each node. The element shape and network structure of the triangular mesh are iteratively optimized to adapt to the boundary of the slope mechanics model, presenting a quasi-regular array-type approximate equilateral triangular mesh. The key areas of the model are divided into six units for densification, forming a quasi-regular triangular mesh with six-directional continuity. The six-directional mesh data is further converted into a format of fifteen-DOF six-node triangular high-order plastic elements. Non-negative auxiliary variables are introduced to fill the common edges of the elements as velocity discontinuities, establishing a finite element second-order cone programming mathematical model for limit analysis. The solution process corresponds to the search for the limit state of slope plastic flow instability. The dissipated energy density cloud map jointly contributed by the obtained triangular elements and velocity discontinuities can easily identify the potential sliding surface of the slope. Specifically, it includes the following steps:

[0067] S1. Construct the initial point set of the slope model according to the preset unit size, perform Delaunay triangulation to generate the initial mesh, analogize the triangulation to a truss structure, apply Hooke's law to calculate the resultant force vector of all truss elastic member nodes, adjust the node coordinates along the direction of the resultant force vector, and iteratively implement node coordinate optimization and mesh reconstruction until an approximately equilateral triangular mesh that adapts to the boundary of the slope model and is arranged in a quasi-regular array is generated.

[0068] In this embodiment, the step specifically includes:

[0069] S11, according to the preset initial grid size h 0. Generate an initial large-area regular array node set, filter the internal nodes of the slope model, merge the coordinates of fixed points on the same boundary, and construct the initial point set;

[0070] S12, perform Delaunay triangulation on the input point set, update the internal mesh structure of the model, sort the vertex index of the cell edge in ascending order, and remove duplicate edges;

[0071] S13, let the element edges satisfy Hooke's law, and assume the equivalent spring stiffness is... K Calculate the element nodes using the following formula. i The resultant force at the point:

[0072]

[0073] in, For nodes i exist x Component of force in direction, For nodes i exist y Component of force in direction, For nodesi With nodes j The distance between them , They are nodes i , j coordinates For nodes i With nodes j The expected length between them is calculated using the following formula:

[0074]

[0075] in, It is the sum of the squares of the lengths of all cell edges in the grid. This is the sum of the squares of the lengths of all cell edges in the desired mesh.

[0076] S14, performs iterative updates of node coordinates based on explicit time integration, causing the nodes to shift along the direction of the resultant force vector. The iterative formula is shown below:

[0077]

[0078] in, , The node after the previous iteration i of x, y coordinate, , For the node after iteration i of x, y coordinate, For the force balance iteration step size, it can generally be taken as... =0.01.

[0079] S15, for nodes that exceed the boundary after offset, the implicit distance function along the boundary is used. Projecting along the negative gradient direction, we constrain it back to the boundary:

[0080]

[0081] in, g Distance function right x and y The sum of squares of partial derivatives.

[0082] implicit distance function of the boundary It relates to the nearest distance from the node to the boundary. For a circular boundary, the specific expression is:

[0083]

[0084] in, Let the coordinates be the center of the circle. The radius is defined in this embodiment. The boundary function value inside the visible model is set to be negative, 0 on the boundary, and positive outside the boundary.

[0085] The specific expression for the rectangular boundary is:

[0086]

[0087] in, , The coordinates of the bottom left and top right vertices of the rectangular boundary are given, and the boundary function value inside the model is negative. Other complex boundaries can be obtained by performing Boolean operations on simple boundaries, which will not be detailed in this embodiment.

[0088] S16, repeat S12~S15 to perform node coordinate offset and mesh reconstruction, and record the sum of squared offsets of all nodes inside the model in each loop. .

[0089] S17, stop when the sum of the squares of the offsets of all nodes in the computational domain is less than the threshold or the preset number of iterations is met. At this time, the system is approximately in force equilibrium, and an approximate equilateral triangular mesh with a quasi-regular array arrangement adapted to the boundary is generated inside the model.

[0090] S2 applies a local densification strategy driven by densification edges to the key areas of the approximately equilateral triangular mesh, forming a six-directional topological mesh with a local array arrangement. That is, the key areas of the slope model are divided into six densifications, the transition area is subjected to a passive densification strategy based on the number of densification edges, and the undensed areas retain the original units.

[0091] It should be noted that, in this embodiment, the key area is defined as an extended area determined through empirical identification based on the slope's geometric characteristics. This area must completely cover the spatial distribution range of the potential slip surface. Furthermore, for complex conditions such as multi-level slopes and other irregular slopes, the area can be progressively extended inwards from the actual terrain contour to determine the key area.

[0092] In this embodiment, this step specifically includes the following sub-steps:

[0093] S21, Traverse all elements and select triangular elements located within the key region based on the centroid coordinate discrimination criterion:

[0094]

[0095]

[0096] in, For unit centroid coordinates, This indicates how to obtain the centroid coordinates. For key areas, Encryption unit to meet encryption conditions Index set.

[0097] S22, extract and mark adjacent pre-encrypted edge data according to the encryption unit index, traverse all units to count the number of pre-encrypted edges contained in them, and mark four types of encryption edge states: unit with no encryption edge, unit with one encryption edge, unit with two encryption edges, and unit with three fully encrypted edges.

[0098] S23. For elements without encrypted edges, the original mesh topology is preserved, and no encryption operation is performed. When an element contains only one encrypted edge, the midpoint splitting method is used. A new node is inserted at the midpoint of the encrypted edge, and the new node is connected to the vertex of the unencrypted edge, splitting the original element into two sub-elements. When an element contains two encrypted edges, a three-way split encryption is performed: the longest encrypted edge is identified, its midpoint is taken as the pivot node, and the pivot node is connected to the midpoint and corresponding vertex of the adjacent encrypted edge, generating three sub-elements. When all three edges of an element need to be encrypted, a one-to-six encryption is performed: the centroid node of the element is calculated as the core node, and the centroid node is connected to the midpoint and original vertex of each edge in a counterclockwise order, and the original element is deleted, generating six sub-elements, with the topology as follows: Figure 2 , Figure 3 As shown, ensure that the order of the sub-unit nodes is consistent.

[0099] S3 extracts the coordinate matrix and element information matrix based on a six-directional triangular mesh, generates a six-node triangular element with fifteen degrees of freedom through edge midpoint interpolation, constructs a high-order discontinuous line on the common edge, adds non-negative auxiliary variables and constraints to realize the velocity jump between elements, and forms the double six-discontinuous triangular mesh data format required for the upper limit finite element method of limit analysis.

[0100] In this embodiment, the step specifically includes:

[0101] S31, extract the cell coordinate matrix [P] and cell information matrix [Q] as the output of the mesh processing, where [Q] stores the three vertex information of each cell.

[0102] S32, obtain all element edge numbers, and generate a six-node triangular high-order element with fifteen degrees of freedom by interpolating the edge midpoints. The node order is as follows: Figure 4 As shown, the velocity at any point within the cell It can be expressed as a quadratic function of the nodal coordinates using shape functions:

[0103]

[0104] in, , They are nodes i exist x and y The velocity component in the direction, It is a unit shape function.

[0105] Plastic multiplier at any point within the unit Obtained by interpolation of the three vertex values:

[0106]

[0107] in, For nodes i Plastic multipliers at the location, such as Figure 4 As shown.

[0108] S33, all element edges in the model are categorized into two types based on boundary identification: common edges and boundary edges. Higher-order velocity discontinuities are constructed on the common edges of the elements. The velocities and plastic multipliers of the nodes on both sides of the discontinuity exhibit abrupt changes, such as... Figure 5 As shown, information such as the element numbers on both sides of the velocity discontinuity line and the positions of the discontinuity line nodes on both sides of the element nodes are searched, and auxiliary variables are attached to the velocity discontinuity line. ,and The function value on the velocity discontinuity line should always be non-negative in order to eliminate the absolute value in the associated flow rules.

[0109] It should be noted that before introducing auxiliary variables, the velocities at the nodes of the discontinuous line must satisfy the associated flow law, and since they contain absolute values, linear programming is not feasible. Therefore, auxiliary variables are introduced, including... And all of them are greater than 0, so as to construct a correlation flow rule without absolute values.

[0110] S4. Based on the upper bound theorem of limit analysis, a second-order cone programming model for slope stability is established to search for the critical instability state of plastic flow in the slope. The model outputs the point cloud of dissipated energy density distribution with high-order triangular elements and velocity discontinuities isomorphic, and obtains the potential sliding surface of the slope limit state.

[0111] It should be noted that the critical instability state of plastic flow on a slope is obtained through the upper bound theorem of limit analysis. That is, when a slope reaches its limit state, there must exist an infinite number of permissible velocity fields such that the internal dissipation energy does not exceed the work done by the external forces.

[0112]

[0113] in, For internal energy dissipation, and These are the plastic stress and strain vectors, respectively. Let be the velocity vector along the discontinuity line. Doing work for external forces and All The function, Let be the yield function. This is the computational domain of the model.

[0114] In this embodiment, the internal dissipation energy includes two parts: the internal dissipation energy of the six-node triangular element and the dissipation energy on the higher-order velocity discontinuity line, namely:

[0115]

[0116] in, For the distribution domain of plastic units, For velocity jumps on discontinuous lines, For velocity discontinuity line sets.

[0117] By applying the upper bound theorem, we can establish a minimum objective function, resulting in the following second-order cone programming (SOCP) model:

[0118]

[0119] in, For plastic flow constraints on the element and velocity discontinuity line, For boundary velocity constraints, A To satisfy the coefficient matrix corresponding to the associated flow rule, B This is the coefficient matrix corresponding to the boundary constraint equations satisfied by the nodes at the boundary. It is a column vector composed of unknown variables such as the velocity of all unit nodes in the model, auxiliary variables, and plasticity multipliers. Add constraints to the second-order cone programming, parameters , , Let Y be the angle between a point on the yield surface of the Mohr-Coulomb yield criterion and the Y-axis. For speed constraints related to self-weight, This is the load matrix.

[0120] In this embodiment, the dissipation density distribution point cloud is displayed by interpolating the model's dissipation energy density to the element nodes. Wherein, the element... k No. i Energy dissipation density at each vertex node It can be represented as:

[0121]

[0122] in, For unit k area, j For unit k In and nodes i Adjacent vertex nodes, c , These represent the cohesion and internal friction angle of the soil, respectively. For the firstk The plastic multiplier vector of each unit, i.e. , a For constant coefficient terms, a =[1,1,4,4,1,1] T , For nodes i, j The corresponding auxiliary variables.

[0123] It should be noted that in this embodiment, the potential sliding surface of the slope is characterized by a normalized dissipated energy density point cloud, i.e. (Values ​​range from 0 to 1) This represents the maximum energy dissipation density of the model.

[0124] Therefore, the potential sliding surface of the slope limit state can be obtained by using the point cloud of dissipated energy density distribution.

[0125] The effects of the present invention are further illustrated below through specific examples. It is assumed that the rock and soil mass is a homogeneous isotropic material, and a model is established as follows: Figure 6 The slope stability analysis model is shown below. The slope gradient is... =70°, cohesion is =10 kPa, internal friction angle is =35°, critical unit weight is The slope height is =10m, stability coefficient .

[0126] First, the computational domain is discretized using a local quasi-regular double-six discontinuous triangular mesh: This is established via S1. Figure 7 The initial mesh shown is an approximate equilateral triangle arranged in a quasi-regular array to adapt to the model boundary, where the number of elements is... =2040; Set the vertical boundary parameters via S2. =1.5 Horizontal boundary parameters =2 In the key region, the initial mesh is divided into six sub-mesh units, and the slope model is discretized into an array of six-directional topological meshes. At this point, the number of elements is... =8570, such as Figure 8 As shown.

[0127] Then, a double six-discontinuous triangular mesh is constructed using S3, consisting of six-node triangular elements and higher-order discontinuities. Full constraints are applied to the left and right boundaries and the bottom, while the slopes and top surface are free. Plastic-related flow constraints are considered within the higher-order elements and at the velocity discontinuities. The six-node velocity... v And the three-vertex plastic multiplier, and the corresponding auxiliary variable , To optimize the variables, the total number of velocity discontinuities, nodes, and optimization variables in the model are 15452, 5671, and 604312, respectively. An upper limit finite element second-order cone programming model is established.

[0128] Finally, using the developed second-order cone programming optimization solution program, based on S4, the solution is obtained for a given number of elements. When the slope model reaches a critical instability state at a value of 8570, the upper limit solution of the stability coefficient is obtained. =14.13, compared to the upper bound solution obtained by existing technologies at the same mesh density. =14.45, with an error of less than 3%. For the upper limit solution, the smaller the result, the closer it is to the true solution, which proves the computational advantage of this method in searching for potential sliding surfaces of slopes.

[0129] To further illustrate the computational advantages of this method in upper bound finite element analysis of slopes, a global double-six discontinuous mesh and a Delaunay mesh are introduced for comparison. The normalized dissipative energy density point cloud distribution of the model under different mesh partitioning modes is plotted. It can be seen that the dissipated energy corresponding to the three mesh partitioning modes is concentrated at the slope toe. The difference lies in the fact that the local and global double-six discontinuous meshes exhibit smooth transitions within the potential sliding surface and the extended range on both sides, and the high dissipative energy density bands are narrower than those of the Delaunay mesh, indicating that the double-six discontinuous mesh can more accurately search for the potential sliding surface of the slope. Simultaneously, the comparison reveals that the local and global double-six discontinuous meshes, under the same initial mesh density conditions, yield the same upper bound solution and dissipated energy density cloud map, indicating that this invention only performs a one-to-six mesh refinement operation on the key areas of the model, which can reduce the computational burden while ensuring the computational accuracy of the model.

[0130] To further illustrate the discretization effect of the six-directional discontinuous lines of this invention on discontinuous velocity fields, point cloud maps of dissipated energy density under three grid discretization modes are output. It was found that in both double-six-discontinuous grids, a distinct crack propagating from the slope toe to the slope surface can be identified in the point cloud of dissipated energy distribution along the dissipated lines, with the dissipated energy density symmetrically distributed on both sides of the crack. However, the crack generated at the slope toe by the Delaunay grid is not obvious. This indicates that the quasi-regular double-six-discontinuous triangular mesh in this invention is essentially aligned with the direction of the potential sliding surface of the slope, which helps to improve the search accuracy of the slope sliding surface.

[0131] Based on the same inventive concept, this embodiment also provides an apparatus, which includes at least one processor and a memory communicatively connected to the at least one processor. The memory stores instructions executable by the at least one processor, which, when executed, enables the at least one processor to perform the aforementioned method for constructing a quasi-regular double-six discontinuous triangulation network for searching potential slope sliding surfaces.

[0132] Based on the same inventive concept, this embodiment also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the aforementioned method for constructing a quasi-regular double-six discontinuous triangular network for searching potential slope sliding surfaces.

[0133] The computer-readable medium includes, but is not limited to, any type of disk (including floppy disks, hard disks, optical disks, CDs). ROM (and magneto-optical disk), ROM, RAM, EPROM (Erasable Programmable Read-Only ROM) Computer-readable media includes erasable programmable read-only memory (EPROM), flash memory, magnetic cards, or optical cards. In other words, this computer-readable medium includes any medium by which a device (such as a printer) stores or transmits information in a readable form.

[0134] The apparatus, computer-readable storage medium, etc. provided in this embodiment have the same inventive concept and the same beneficial effects as the aforementioned method, and will not be described again here.

[0135] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0136] The above embodiments are merely illustrative of several implementation methods of this application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.

Claims

1. A method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes, characterized in that, Includes the following steps: S1. Construct the initial point set of the slope model according to the preset unit size, perform Delaunay triangulation to generate the initial mesh, compare the triangulation mesh with the truss structure, apply Hooke's law to calculate the resultant force vector of all truss elastic members, adjust the node coordinates along the direction of the resultant force vector, and iteratively implement node coordinate optimization and mesh reconstruction until an approximate equilateral triangular mesh that adapts to the boundary of the slope model and is arranged in a quasi-regular array is generated. In S1, node coordinates are iteratively updated based on explicit time integration, causing the nodes to shift along the direction of the resultant force vector. For nodes that have shifted beyond the boundary, they are projected along the negative gradient direction of the implicit distance function of the boundary to constrain them back onto the boundary. S2, a local densification strategy driven by densification edges is applied to the key areas of the approximately equilateral triangular mesh to form a six-directional triangular mesh with a local array arrangement. In the key areas of the slope model, a one-to-six densification is applied. In the transition area, a passive densification strategy is applied based on the number of densification edges. The original units are retained in the non-densification area. S2 includes the following sub-steps: S21, Traverse all units and select triangular units located within the key area based on the centroid coordinate discrimination criterion; S22, extract and mark adjacent pre-encrypted edge data according to the encryption unit index, traverse all units to count the number of pre-encrypted edges contained in them, and mark four types of encryption edge states: no encryption edge unit, one encryption edge unit, two encryption edge units, and three fully encrypted edge units. S23. For cells without encrypted edges, the original mesh topology is preserved, and no encryption operation is performed. When a cell contains only one encrypted edge, the midpoint splitting method is used. A new node is inserted at the midpoint of the encrypted edge, and the new node is connected to the vertex of the non-encrypted edge, splitting the original cell into two sub-cells. When a cell contains two encrypted edges, three-way splitting encryption is performed: the longest encrypted edge is identified, its midpoint is taken as the pivot node, and the pivot node is connected to the midpoint and corresponding vertex of the adjacent encrypted edge to generate three sub-cells. When all three edges of a cell need to be encrypted, a one-to-six encryption is performed: the centroid node of the cell is calculated as the core node, and the centroid node is connected to the midpoint and original vertex of each edge in a counterclockwise order, and the original cell is deleted to generate six sub-cells. S3, based on the six-direction triangulation, extracts the coordinate matrix and element information matrix, generates a six-node triangular element with fifteen degrees of freedom by interpolating the midpoint of the edge, constructs a high-order discontinuity line on the common edge, adds non-negative auxiliary variables and constraints to realize the velocity jump between elements, and forms the double six-discontinuity triangulation data format required for the upper limit finite element of limit analysis. S3 includes the following sub-steps: S31, extract the cell coordinate matrix [P] and cell information matrix [Q] as the output of the mesh processing. [Q] stores the information of the three vertices of each cell. S32, obtain all element edge numbers, and generate a six-node triangular high-order element with fifteen degrees of freedom by interpolating the edge midpoints; S33, all element edges in the model are identified as common edges and boundary edges according to the boundary identifier. High-order velocity discontinuity lines are constructed on the common edges of the elements. There are jumps in the velocity and plastic multipliers of the nodes on both sides of the discontinuity line. The element numbers on both sides of the velocity discontinuity line and the positions of the nodes of the discontinuity line in the nodes of the elements on both sides are searched, judged and established. Auxiliary variables are introduced on the velocity discontinuity line to remove the absolute values ​​in the associated flow rules. S4. Based on the upper bound theorem of limit analysis, a second-order cone programming model for slope stability is established to search for the critical instability state of plastic flow in the slope. The model outputs the point cloud of dissipated energy density distribution with high-order triangular elements and velocity discontinuities isomorphic, and obtains the potential sliding surface of the slope limit state.

2. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 1, characterized in that, S1 includes: According to the preset initial grid size h 0. Generate an initial large-area regular array node set, filter the internal nodes of the slope model, merge the coordinates of fixed points on the same boundary, and construct the initial point set; Perform Delaunay triangulation on the input point set, update the internal mesh structure of the model, sort the vertex index of the cell edge in ascending order, and remove duplicate edges; Let the element edges satisfy Hooke's law, and assume the equivalent spring stiffness is... K , computational unit node i The resultant force at the point: in, For nodes i exist x Component of force in direction, For nodes i exist y Component of force in direction, For nodes i With nodes j The distance between them Adjacent nodes, , They are nodes i , j coordinates For nodes i With nodes j Expected length between: in, It is the sum of the squares of the lengths of all cell edges in the grid. This is the sum of the squares of the lengths of all cell edges in the desired mesh. For nodes that offset beyond the boundary along the resultant force vector direction, the implicit distance function along the boundary is... Projecting along the negative gradient direction, we constrain it back to the boundary: in, g Distance function right x and y The sum of squares of partial derivatives; The algorithm iterates through node coordinate offsets and mesh reconstruction, recording the sum of squared offsets of all nodes within the model during each iteration. The process stops when the sum of the squares of the offsets of all nodes in the computational domain is less than a threshold or when the preset number of iterations is met. An approximate equilateral triangular mesh with a quasi-regular array arrangement adapted to the boundary is generated inside the model.

3. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 2, characterized in that, The velocity at any point within the nodal triangle high-order element in S32. Represented as a quadratic function of the nodal coordinates using shape functions: in, , They are nodes i exist x and y The velocity component in the direction, For unit shape functions; Plastic multiplier at any point within the unit Obtained by interpolation of the three vertex values: in, For nodes i Plastic multipliers at the location.

4. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 3, characterized in that, The critical instability state of the slope in S4 is obtained through the upper bound theorem of limit analysis. When the slope reaches the limit state, there are infinitely many permissible velocity fields such that the internal dissipated energy is no greater than the work done by the external forces. in, For internal energy dissipation, and These are the plastic stress and strain vectors, respectively. Let be the velocity vector along the discontinuity line. Doing work for external forces and All The function, Let be the yield function. This is the computational domain of the model.

5. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 4, characterized in that, In S4, the objective function is minimized using the upper bound theorem, forming a second-order cone programming model: in, For plastic flow constraints on the element and velocity discontinuity line, For boundary velocity constraints, A To satisfy the coefficient matrix corresponding to the associated flow rule, B This is the coefficient matrix corresponding to the boundary constraint equations satisfied by the nodes at the boundary. It is a column vector consisting of the velocities of all unit nodes in the model, auxiliary variables, and unknown plastic multipliers. Add constraints to the second-order cone programming, parameters , , Let be the angle between a point on the yield surface of the Mohr-Coulomb yield criterion and the Y-axis. For speed constraints related to self-weight, This is the load matrix.

6. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 4, characterized in that, The internal dissipation energy consists of two parts: the internal dissipation energy of the six-node triangular element and the dissipation energy on the higher-order velocity discontinuity line. in, For the distribution domain of plastic units, For velocity jumps on discontinuous lines, For velocity discontinuity line sets.

7. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 4, characterized in that, The dissipation density distribution point cloud is displayed by interpolating the model's dissipation energy density to the element nodes. k No. i Energy dissipation density at each node Represented as: in, For unit k area, j For unit k In and nodes i Adjacent vertex nodes, c , These represent the cohesion and internal friction angle of the soil, respectively. For the first k Plastic multiplier vectors of each unit , a For constant coefficient terms, For nodes i, j The corresponding auxiliary variables.

8. The method for constructing a quasi-regular double-six discontinuous triangular network for searching potential sliding surfaces of slopes according to claim 7, characterized in that, The potential sliding surface of the slope is characterized by a point cloud of normalized dissipative energy density: in, This represents the maximum energy dissipation density of the model.

Images