Method for determining asymmetric three-dimensional potential slip surface using maximum shear strain increment

Abstract

The application relates to a method for determining an asymmetric three-dimensional potential sliding surface by using maximum shear strain increment. The method comprises the following steps: establishing a geomechanics model and arranging equidistant vertical lines and measuring points, calculating the position of the measuring point with the maximum shear strain increment on each vertical line; determining the circular in-group to which all the measuring points in the measuring point set with the maximum shear strain increment belong through iteration; forming a preliminary slope surface exposed sliding surface boundary line according to the shape presented by the measuring points in the circular in-group, supplementing the missing boundary measuring points, and obtaining the supplemented slope surface exposed sliding surface boundary line; searching for all the measuring points with the maximum shear strain increment located within the supplemented slope surface exposed sliding surface boundary line, grouping the measuring points with the maximum shear strain increment and the measuring points of the supplemented slope surface exposed sliding surface boundary line, and then sequentially connecting the measuring points to obtain a three-dimensional potential sliding surface; and calculating the stability safety factor of the potential sliding surface. The search result accuracy of the slope sliding surface is improved.

Description

Technical Field

[0001] This application relates to the field of geotechnical engineering technology for slopes and tunnels, and in particular to a method for determining an asymmetric three-dimensional potential sliding surface using the maximum shear strain increment. Background Technology

[0002] Slope stability analysis has always been a hot research topic in the field of geotechnical engineering. With the development of the national economy, a large number of complex slope engineering problems are constantly emerging in the construction of related highways, railways, water conservancy, ports, civil engineering and other industries. Conducting slope stability analysis research has profound theoretical significance and broad engineering application prospects.

[0003] Common methods for slope stability analysis include the limit equilibrium method, limit analysis method, and numerical analysis method. The traditional limit equilibrium method is widely used, but it requires assumptions about the shape of the slope sliding surface and is not suitable for analyzing composite sliding surfaces such as rock slopes; moreover, it skips the stress-strain relationship in the slope and does not consider the strength conditions of the slope soil, thus failing to provide key information about displacement.

[0004] With the development of computer technology, numerical analysis methods such as the finite difference method have been widely applied. Many scholars have also developed a set of methods that use only numerical methods to determine the slope slip surface, such as the finite element strength reduction method, the equivalent plastic strain method, the stress influence coefficient method, and the variational method. However, most current three-dimensional slip surface search methods are derived from the extension of two-dimensional slip surface search methods. They work well when applied to symmetrical slip surfaces, but when it comes to true three-dimensional slip surface searches, it is difficult to determine the shear inlet and shear outlet, resulting in significant differences between the search results and the actual slip surface. In addition, when calculating the ultimate limit state of a slope, it is necessary to define the plastic failure slip surface that runs from the toe to the crest. Currently, the general approach is to first determine the approximate range of the slip surface and then smooth it through curve fitting. This further restricts the shape of the slip surface, making it impossible to effectively search for slip surfaces with non-differentiable points, such as polygonal slip surfaces, resulting in low accuracy of the slope slip surface search results. Summary of the Invention

[0005] Therefore, it is necessary to provide a method for determining asymmetric three-dimensional potential sliding surfaces by utilizing the maximum shear strain increment, which can improve the accuracy of search results for slope sliding surfaces and address the aforementioned technical problems.

[0006] A method for determining an asymmetric three-dimensional potential sliding surface using the maximum shear strain increment, the method comprising:

[0007] Step S1: Based on the FLAC3D continuous numerical simulation platform, establish a geomechanical model that considers complex strata distribution and consists of nodes and elements, and initialize it. The initialization includes: assigning deformation and strength parameters, applying self-weight external load and boundary conditions to calculate the initial stress balance, and clearing displacement and velocity to zero.

[0008] Step S2: Perform excavation numerical simulation on the geomechanical model until the stress equilibrium state is reached, retain the strain field results at the stress equilibrium state, arrange equally spaced vertical lines on the geomechanical model, and arrange measuring points at equally spaced intervals from bottom to top along the vertical lines. Based on the strain field results at the stress equilibrium state, calculate the location of the measuring point with the maximum shear strain increment on each vertical line.

[0009] Step S3: Group the measurement points based on whether they are located on the slope surface, and assign the measurement points located on the slope surface to the slope surface group and the measurement points located inside the slope to the slope inside group.

[0010] Step S4: Divide all the measuring points in the slope table group into segments along the X-axis and Y-axis according to continuity, determine the set of measuring points with the maximum value of shear strain increment, and discard the remaining measuring points in the slope table group;

[0011] Step S5: For the measurement points in the set of measurement points with the maximum value of shear strain increment, use a clustering algorithm to group any two measurement points whose distance between them is less than the threshold δ into a group, and determine the inner circle group to which all measurement points belong through iteration to obtain each inner circle group;

[0012] Step S6: Sequentially add the measuring points with the largest and second largest shear strain increment values ​​in each row and column of the bottom surface of the geomechanical model along the X-axis and Y-axis to the retained measuring point set, and delete the remaining measuring points in the inner circle group; if there is only one point in a row or column, only the point with the largest shear strain increment value in that row or column is retained and added to the retained measuring point set.

[0013] Step S7: Then perform step S5 again on the measured points in the retained set of measured points to update the grouping of all points and obtain the updated inner circle groups;

[0014] Step S8: Based on the shape of the measuring points in each updated inner circle group, find the inner circle group combination that is most likely to form a closed loop. Determine this inner circle group combination as the inner circle group combination that constitutes the boundary line of the exposed slip surface of the slope. Discard the other inner circle groups. According to the boundary line connection order, in the measuring points of two adjacent inner circle groups, find the point pair that is closest to a measuring point in one inner circle group and a measuring point in another inner circle group. In the point pair, set the previous measuring point as the starting point and the next measuring point as the ending point according to the boundary line connection order.

[0015] Step S9: Using an improved maze algorithm, the measuring points of each circle group are selectively connected end-to-end according to the starting and ending points. Then, the starting and ending points of adjacent groups are connected according to the connection order, and unconnected measuring points are discarded to form the preliminary slip surface boundary line of the slope surface.

[0016] Step S10: Sort the points connected on the preliminary slip surface boundary line of the slope surface according to the connection order, and determine whether the two adjacent measuring points are continuous along the X-axis and Y-axis. If they are not continuous, supplement the missing boundary measuring points according to the direction of the line connecting the two measuring points to obtain the supplemented slip surface boundary line of the slope surface.

[0017] Step S11: Search for all measurement points with the maximum shear strain increment within the slip surface boundary line exposed by the supplemented slope surface, group them together with the measurement points of the slip surface boundary line exposed by the supplemented slope surface, and then connect them sequentially to obtain a three-dimensional potential slip surface.

[0018] Step S12: Traverse all polygons that constitute the three-dimensional potential sliding surface, calculate the average stress tensor of each polygon, and accumulate the stability safety factor of the potential sliding surface.

[0019] The above method for determining the asymmetric three-dimensional potential sliding surface using the maximum shear strain increment has the following beneficial effects:

[0020] (1) This application addresses the problem of determining the asymmetric complex three-dimensional sliding surface of slopes. It aims to rely on the finite difference numerical simulation platform and utilize the maximum shear strain increment distribution of the model to construct a series of algorithms to select and optimize the sliding surface boundary line exposed on the slope surface, which can solve the problem of characterizing the three-dimensional potential sliding surface of slopes.

[0021] (2) In this application, the maximum value points of the slope surface shear strain increment are grouped and optimized to eliminate the influence of the sliding surface interference section and to evaluate the stability of different regions of the model. Therefore, more realistic and reliable simulation results can be obtained, improving the accuracy of the obtained three-dimensional potential sliding surface, which can be further used to calculate the safety factor of the potential sliding surface. Attached Figure Description

[0022] Figure 1 This is a flowchart illustrating a method for determining an asymmetric three-dimensional potential sliding surface using the maximum shear strain increment in one embodiment.

[0023] Figure 2 This is a schematic diagram of a geomechanical model based on FLAC3D slope excavation numerical values ​​in one embodiment.

[0024] Figure 3 This is a contour plot of the calculated geostress equilibrium shear strain increment in one embodiment;

[0025] Figure 4 This is a schematic diagram of a geomechanical model in one embodiment, in which vertical lines are arranged at equal intervals along the average normal direction of the slope.

[0026] Figure 5This is a schematic diagram of a geomechanical model for calculating and grouping the maximum points of shear strain increments along each vertical line in one embodiment.

[0027] Figure 6 This is a schematic diagram illustrating the point where the maximum value of the shear strain increment for each segment in each row / column of the slope table is obtained in one embodiment.

[0028] Figure 7 This is a schematic diagram of a geomechanical model that uses a clustering algorithm to group slope surface points in one embodiment;

[0029] Figure 8 This is a schematic diagram of a geomechanical model for optimizing slope surface points based on the incremental values ​​of shear strain in each row / column, as shown in one embodiment.

[0030] Figure 9 This is a schematic diagram of a geomechanical model for determining the start and end positions of each group of slope surface points in one embodiment;

[0031] Figure 10 This is a schematic diagram of a geomechanical model in one embodiment, which uses an improved maze algorithm to connect and form the slip surface boundary line exposed on the slope.

[0032] Figure 11 A schematic diagram of a geomechanical model for supplementing missing points on the slip surface boundary line exposed on a slope in one embodiment;

[0033] Figure 12 This is a schematic diagram of a geomechanical model for locating the bottom slip surface point of a slope in one embodiment;

[0034] Figure 13 This is a schematic diagram of a geomechanical model in one embodiment of sequentially connecting potential three-dimensional sliding surfaces. Detailed Implementation

[0035] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0036] In one embodiment, such as Figure 1 As shown, a method for determining an asymmetric three-dimensional potential sliding surface using the maximum shear strain increment is provided, including the following steps:

[0037] Step S1: Based on the FLAC3D continuous numerical simulation platform, establish a geomechanical model consisting of nodes and elements that considers complex strata distribution and perform initialization. Initialization includes: assigning deformation and strength parameters, applying self-weight external loads and boundary conditions to calculate the initial geostress equilibrium, and clearing displacement and velocity to zero.

[0038] Step S2: Perform excavation numerical simulation on the geomechanical model until the geostress equilibrium state, retain the strain field results at the geostress equilibrium state, arrange equally spaced vertical lines on the geomechanical model, and arrange measuring points at equally spaced intervals from bottom to top along the vertical lines. Based on the strain field results at the geostress equilibrium state, calculate the location of the measuring point with the maximum shear strain increment on each vertical line.

[0039] Step S3: Group the measuring points based on whether they are located on the slope surface. The measuring points located on the slope surface are assigned to the slope surface group, and the measuring points located inside the slope are assigned to the inside slope group.

[0040] Step S4: Divide all the measuring points in the slope table group into segments along the X-axis and Y-axis according to continuity, determine the set of measuring points with the maximum value of shear strain increment, and discard the remaining measuring points in the slope table group.

[0041] Step S5: For the measurement points in the set of measurement points with the maximum shear strain increment, use a clustering algorithm to group any two measurement points whose distance between them is less than the threshold δ into a group. Iterate to determine the inner circle group to which all measurement points belong, and obtain each inner circle group.

[0042] Step S6: Sequentially add the measuring points with the largest and second largest shear strain increment values ​​in each row and column of the bottom surface of the geomechanical model along the X-axis and Y-axis to the retained measuring point set, and delete the remaining measuring points in the inner circle group; if there is only one point in a row or column, only the point with the largest shear strain increment value in that row or column is retained and added to the retained measuring point set.

[0043] Step S7: Then perform step S5 again on the measured points in the retained set of measured points to update the grouping of all points and obtain the updated inner circle groups.

[0044] Step S8: Based on the shape of the measuring points in each updated inner circle group, find the inner circle group combination that is most likely to form a closed loop. Determine this inner circle group combination as the inner circle group combination that constitutes the boundary line of the exposed slip surface of the slope. Discard the other inner circle groups. According to the boundary line connection order, in the measuring points of two adjacent inner circle groups, find the point pair that is closest to a measuring point in one inner circle group and a measuring point in another inner circle group. In the point pair, set the previous measuring point as the starting point and the next measuring point as the ending point according to the boundary line connection order.

[0045] Step S9: Using the improved maze algorithm, selectively connect the measurement points of each circle group end-to-end according to the start and end points. Then connect the start and end points of adjacent groups according to the connection order, discarding unconnected measurement points to form the preliminary slip surface boundary line of the slope surface.

[0046] Step S10: Sort the points connected on the initial slope surface exposed slip surface boundary line according to the connection order, and determine whether the two adjacent measuring points are continuous along the X-axis and Y-axis. If they are not continuous, supplement the missing boundary measuring points according to the direction of the line connecting the two measuring points to obtain the supplemented slope surface exposed slip surface boundary line.

[0047] Step S11: Search for all measurement points with the maximum shear strain increment within the slip surface boundary line exposed on the supplemented slope surface, group them together with the measurement points of the slip surface boundary line exposed on the supplemented slope surface, and then connect them sequentially to obtain the three-dimensional potential slip surface.

[0048] Step S12: Traverse all polygons that constitute the three-dimensional potential sliding surface, calculate the average stress tensor of each polygon, and accumulate the stability safety factor of the potential sliding surface.

[0049] The method described above for determining the asymmetric three-dimensional potential sliding surface using the maximum shear strain increment aims to search for the deformation field generated under different external load variations using a numerical model established on the FLAC3D finite difference numerical simulation platform. It searches for the point with the maximum shear strain increment on each vertical line of the slope, and uses the built-in FISH language to program and determine the boundary line of the exposed sliding surface on the slope. Then, all points within this range are sequentially connected to form a surface to obtain the three-dimensional potential sliding surface of the slope. The landslide stability coefficient under this condition is then calculated to evaluate the slope stability.

[0050] In one embodiment, step S2 includes:

[0051] Step S21: Perform excavation numerical simulation on the geomechanical model until the stress equilibrium state is reached, and retain the strain field results at the stress equilibrium state.

[0052] Step S22: Set the average normal direction of the slope surface of the geomechanical model as the vertical direction, and denote the unit vector in the vertical direction as (v1, v2, v3); arrange vertical lines at equal intervals along the X-axis and Y-axis, with the vertical lines extending from the slope surface to the boundary of the geomechanical model, and the vertical line spacing is smaller than the average size of the unit.

[0053] Step S23: Create a two-dimensional array (xx, yy) to represent each vertical line, where xx is the X-coordinate of the intersection of the vertical line and the bottom surface of the geomechanical model, yy is the Y-coordinate of the intersection of the vertical line and the bottom surface of the geomechanical model, and the spacing between adjacent vertical lines is dd; arrange measuring points at equal intervals from bottom to top along each vertical line inside the geomechanical model, with the coordinates of the measuring points being (q*v1+xx, q*v2+yy, q*v3+zz), where q is the measuring point number on each vertical line, and zz is the Z-coordinate of the intersection of the vertical line and the bottom surface of the geomechanical model;

[0054] Step S24: Based on the strain field results under the stress equilibrium state, use the built-in FISH function to statistically analyze the shear strain increments at each measuring point on the vertical line, and find the measuring point number q corresponding to the maximum shear strain increment on that vertical line. max The location of the measuring point with the maximum shear strain increment on the corresponding vertical line is (q max *v1+xx,q max *v2+yy,q max *v3+zz); Find the maximum value q of the measuring point number on the vertical line. max2 The corresponding measuring point on the perpendicular line is (q) max2 *v1+xx,q max2 *v2+yy,q max2 *v3+zz).

[0055] In one embodiment, step 3 includes:

[0056] Step S31: Traverse each perpendicular line and compare the measurement point number q on each perpendicular line. max and q max2 The value of q, if max =q max2 If the maximum shear strain increment of the vertical line is measured at the slope surface, then the maximum shear strain increment of the vertical line is measured at the slope surface; otherwise, the maximum shear strain increment of the vertical line is measured at the slope interior.

[0057] Step S32: Assign the measuring points located on the slope surface to the slope surface group and the measuring points located inside the slope to the slope inside group.

[0058] In one embodiment, step S4 includes:

[0059] Step S41: Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point in the slope table group. Divide the group into multiple groups by grouping measuring points with the same yy value. Iterate through each group in ascending order of xx value, using the measuring points within that group as the current measuring point. The current measuring point is determined based on its belonging to the current first segment x. i If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the current first segment x. i If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next first segment x. i+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the first segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the first segments x1, x2, x3...x M In the first group that is traversed, the smallest measurement point of xx value belongs to the first segment x1, i is the first segment number, i∈1,2,3……M, and M is the number of the first segments;

[0060] Step S42: Calculate the shear strain increment of each measuring point in the first segment, and select the measuring point corresponding to the maximum shear strain increment to add to the first measuring point set;

[0061] Step S43: Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of the measuring points located on the slope. Divide the area into multiple groups by grouping measuring points with the same xx value. Then, iterate through each group, selecting the measuring points within that group in ascending order of their yy values, and use these points as the current measuring point. The current measuring point is determined based on whether it belongs to the second segment y. j If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the second segment y. j If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next second segment y. j+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the second segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the second segments y1, y2, y3...y M In the first group that is traversed, the point with the smallest yy value belongs to the second segment y1, j is the number of the second segment, j∈1,2,3……Q, and Q is the number of the second segments;

[0062] Step S44: Calculate the shear strain increment of each measuring point in the second segment, and select the measuring point corresponding to the maximum shear strain increment to add to the second measuring point set;

[0063] Step S45: Take the union of the first set of measuring points and the second set of measuring points as the set of measuring points with the maximum value of shear strain increment, and discard all unselected measuring points in the slope table group.

[0064] In one embodiment, step S5 includes:

[0065] Step S51: First, randomly select one measurement point from the set of measurement points with the maximum shear strain increment as the current measurement point and assign it to the inner circle group g. n Draw a circle with δ as its center and radius as its threshold. Then, group the ungrouped measurement points within the circle that fall within the set of measurement points with the maximum shear strain increment into the inner circle group g. n And will be classified as the inner group g of the circle. n The measuring points are removed from the set of measuring points with the maximum shear strain increment. Among the measuring points in the set of measuring points with the maximum shear strain increment, the first randomly selected measuring point is assigned to the inner circle group g1.

[0066] Step S52: Next, from the inner circle group g n Select the measurement points that are not currently being measured as the current measurement points. Draw a circle with the current measurement point as the center and a radius of δ as the threshold. Assign measurement points that fall within the circle and are not grouped from the set of measurement points with the maximum shear strain increment to the group g within the circle. n And will be classified as the inner group g of the circle.n The measurement points are deleted from the set of measurement points with the maximum shear strain increment, and so on, until there are no measurement points in group g within the circle that have not been selected as the current measurement points; n in group g;

[0067] Step S53: Arbitrarily select a measurement point from the set of measurement points with the maximum shear strain increment and classify it into group g within the circle n+1 , repeat steps S51 and S52 to determine the measurement points classified into group g within the circle n+1 and delete the measurement points classified into group g within the circle n+1 from the set of measurement points with the maximum shear strain increment. By analogy, all the measurement points in the set of measurement points with the maximum shear strain increment are grouped, and groups g1, g2... gN within the circle are obtained, where gN is the last group within the circle.

[0068] In one embodiment, step S9 includes:

[0069] Step S91: Sequentially select the starting point within a group within the circle as the first measurement point within that group, draw a circle with a radius of threshold δ with it as the center, classify all the measurement points within the circle that fall within this circle into the mapping group m1 of the first measurement point, and sort the measurement points within mapping group m1 in ascending order of the distance from each measurement point to the first measurement point;

[0070] Step S92: Select the measurement point a with the smallest distance within mapping group m1 as the second measurement point, draw a circle with a radius of threshold δ with it as the center, and classify all the measurement points within the circle that fall within this circle and except for the first measurement point into the mapping group m a of measurement point a, and sort the measurement points within mapping group m a in ascending order of the distance from each measurement point to measurement point a;

[0071] Step S93: Select the measurement point b with the smallest distance within mapping group m a , compare the distance s2 from the second measurement point to measurement point b with the distance s1 from the first measurement point to measurement point b. If s2 < s1, then select measurement point b as the third measurement point; if s2 ≥ s1, then select the measurement point with the second smallest distance within mapping group m a as measurement point b, and compare the distance s2 from the second measurement point to measurement point b with the distance s1 from the first measurement point to measurement point b again. If s2 < s1, then select measurement point b as the third measurement point; if s2 ≥ s1, then select the measurement point with the third smallest distance within mapping group m a as measurement point b, and make judgments in ascending order of the distance within mapping group m a by analogy until the measurement point b corresponding to s2 < s1 is found as the third measurement point. Among them, if mapping group m aIf all measuring points within the range have a distance s2 that is not less than a distance s1, then return to the mapping group m1 of the first measuring point, select the measuring point with the second smallest distance within the mapping group m1 as the second measuring point, and so on.

[0072] Step S94: Based on the distance within the mapping group of the previous measuring point, select the next measuring point of the previous measuring point in ascending order, until the end point within the circle group is selected, and then connect the selected measuring points in sequence.

[0073] Step S95: Select the starting and ending points of adjacent inner circles and connect them to form the preliminary slip surface boundary line of the slope surface.

[0074] In one embodiment, step S10 includes:

[0075] Step S101: Based on the preliminary exposed slip surface boundary line of the slope, traverse all measuring points on the preliminary exposed slip surface boundary line in the order of boundary line connection. Calculate the two-dimensional arrays corresponding to the vertical lines of pairwise adjacent measuring points P and Q. The array corresponding to measuring point P is denoted as (Rp, Cp), and the array corresponding to measuring point Q is denoted as (Rq, Cq). If (Rq-Rp) / dd > 1, then the array corresponding to the boundary measuring points to be supplemented is (Rs, Cs), where Rs ∈ [Rp+dd, Rq-dd]. The step size interval between the boundary measuring points to be supplemented is dd. The expression for calculating Cs is:

[0076]

[0077] Determine the parameter q of the corresponding perpendicular line based on the two-dimensional array (Rs, Cs). max2 The location of the supplementary boundary measurement point is (q max2 *v1+Rs,q max2 *v2+Cs,q max2 *v3+zz), after completing all boundary measurement points, update the connection order;

[0078] Step S102: Traverse the supplemented boundary measurement points according to the connection order, and count the two-dimensional arrays corresponding to the perpendicular lines of adjacent measurement points M and N. The array corresponding to measurement point M is denoted as (Rm, Cm), and the array corresponding to measurement point N is denoted as (Rn, Cn). If (Cn-Cm) / dd>1, then the array corresponding to the boundary measurement points to be supplemented is (Rt, Ct), where Ct∈[Cm+dd, Cn-dd]. The step size distance between the boundary measurement points to be supplemented is dd. The expression for calculating Rt is:

[0079]

[0080] Determine the parameter q of the corresponding perpendicular line based on the two-dimensional array (Rt, Ct). max2 The location of the supplementary boundary measurement point is (qmax2 *v1+Rt,q max2 *v2+Ct,q max2 *v3+zz), after supplementing all boundary measurement points, update the connection order to obtain the supplemented slip surface boundary line exposed on the slope surface.

[0081] In one embodiment, step S11 includes:

[0082] Step S111: Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point on the exposed slip surface boundary line of the supplemented slope surface. Divide the area into multiple groups, grouping measuring points with the same yy value. Iterate through each group sequentially, selecting the measuring points within that group in ascending order of xx value, and use these points as the current measuring point. The current measuring point is determined based on its belonging to the current third segment x. k If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the current third segment x. k If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next third segment x. k+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the third segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the third segments x1, x2, x3...x L In the first group that is traversed, the smallest measurement point of xx value belongs to the third segment x1, k is the number of the third segment, k∈1,2,3……L, and L is the number of the third segments;

[0083] Step S112: Calculate the x values ​​of the test points within each third segment, and create a two-dimensional array (x1, x2) and a four-dimensional array (y0, y1, y2, y3) to represent each third segment. x1 is the minimum x value of the test points within each third segment, x2 is the maximum x value of the test points within each third segment, y1 is the yy value of the test point corresponding to the x1 value within each third segment, y2 is the yy value of the test point corresponding to the x2 value within each third segment, y0 is the yy value of the test point preceding the y1 value according to the boundary line connection order, and y3 is the yy value of the test point following the y2 value according to the boundary line connection order.

[0084] Step S113: Create a feature variable V to represent each third segment, and obtain the values ​​of y0, y1, y2 and y3 for each third segment. If (y0-y1)×(y3-y2)<0, set the V value of the third segment to 1; otherwise, set the V value to 0.

[0085] Step S114: For the yy values ​​of the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point of the exposed slip surface boundary line of the supplemented slope, traverse the yy values ​​in ascending order and assign an initial feature value T(yy) = 0 to each yy value; for each yy value, traverse the xx values ​​of the two-dimensional array (xx, yy) of all vertical lines of the geomechanical model in ascending order and compare them with the x1 values ​​of all third segments within the corresponding group. When there exists x1 = xx, obtain the V value of the corresponding third segment and update T(yy) to the sum of its current value and the V value; if the updated T(yy) is odd, obtain the corresponding third segment. Given the value x2, update the starting position of the xx value ss = x2 + dd. Starting from x2 + dd, continue to traverse the xx values ​​in ascending order and compare them with the x1 values ​​of all third segments in the corresponding group. When there exists x1 = xx, obtain the V value of the corresponding third segment and update T(yy) to the sum of its current value and V. If the updated T(yy) is even, divide the two-dimensional array (xx, yy) in the interval [ss, x1 - dd] into the bottom sliding surface group, then obtain the x2 value of the corresponding third segment, and continue to traverse the xx values ​​in ascending order starting from x2 + dd, and so on, until the xx values ​​have been traversed.

[0086] Step S115: Traverse the two-dimensional array (xx, yy) divided into the bottom sliding surface group, and determine the parameter q of the corresponding perpendicular line according to each two-dimensional array. max Establish a measuring point, the location of which is (q). max *v1+xx,q max *v2+yy,q max *v3+zz), and assign the measuring points to the bottom sliding surface group;

[0087] Using the quadrilateral formed by the pairwise adjacent xx and yy values ​​of the two-dimensional array (xx, yy) of all vertical lines in the geomechanical model, the measuring points in the quadrilaterals of the bottom slip surface group and the measuring points of the slip surface boundary line exposed on the supplemented slope surface are grouped together. All groups are traversed. If a group contains 3 or 4 measuring points, they are connected in sequence to form a triangle or quadrilateral. The surface formed by all triangles and quadrilaterals is the three-dimensional potential slip surface.

[0088] In one embodiment, step S12 includes:

[0089] Step S121: Traverse all polygons constituting the three-dimensional potential sliding surface. Using the stress values ​​of the m vertices on each polygon in the global coordinate system xyz, determine the average stress tensor of each polygon. The average stress tensor of the k-th polygon is:

[0090]

[0091] In the formula, S ki Let S be the stress tensor of the i-th vertex of the k-th polygon. k Composed of 6 stress components {σ xk ,σ yk ,σ zk ,τ xyk ,τ xzk ,τ yzk} constitute, σ xk σ yk σ zk τ represents the normal stress in the x, y, and z directions of the k-th polygon, respectively. xyk τ xzk τ yzk These are the shear stresses on the xy, xz, and yz planes of the k-th polygon, respectively.

[0092] Step S122: For the k-th polygon, define the z′ axis of the local coordinate system as its normal direction, the x′ axis of the local coordinate system as the average displacement direction of the m vertices, and the y′ axis of the local coordinate system as the cross product direction of z′ and x′. Use the formula S′ k =LS k L T The mean stress tensor is transformed from the global coordinate system to the local coordinate system, where... (l1 m1n1) is the direction cosine of the x′ axis of the local coordinate system in the global coordinate system, (l2 m2 n2) is the direction cosine of the y′ axis of the local coordinate system in the global coordinate system, and (l3 m3 n3) is the direction cosine of the z′ axis of the local coordinate system in the global coordinate system.

[0093] Step S123: Take S′ k The stress component σ′ zk For normal stress, τ′ xzk For tangential stress, the anti-slip force and sliding force of all polygons are summed to determine the stability safety factor of the potential sliding surface. The stability safety factor of the potential sliding surface is:

[0094]

[0095] In the formula, F s The stability safety factor is n, where n is the total number of polygons, and c is the total number of polygons. k , Let A be the cohesion and internal friction angle of the k-th polygon. k Let be the area of ​​the k-th polygon.

[0096] The above-mentioned method for determining the asymmetric three-dimensional potential sliding surface using the maximum shear strain increment relies on the FLAC3D finite difference numerical simulation platform. It constructs a geomechanical model of a three-dimensional slope for numerical calculation, discretizes the slope into a series of vertical lines, searches for the coordinates of the position of the maximum shear strain increment on the vertical lines, and uses the built-in FISH language to write a program to determine the asymmetric three-dimensional potential sliding surface using the maximum shear strain increment. This method can quickly and efficiently determine the morphology and safety factor of complex sliding surfaces of three-dimensional slopes.

[0097] In one embodiment, a rock slope has complex strata and a complex deformation field generated by varying external loads. The method of determining the asymmetric three-dimensional potential sliding surface using the maximum shear strain increment of this application is used to search for the sliding surface. The steps are as follows:

[0098] (1) Based on the continuous numerical simulation platform (finite difference software FLAC3D), a geomechanical model was established that can consider complex stratigraphic distribution and consists of nodes and units. It comprises six layers of material, from top to bottom: carbonaceous phyllite, siliceous slate, mudstone, foliated chlorite-bearing sodic schist, limestone, and chlorite-bearing sodic schist. Each layer is further divided into strongly weathered, weakly weathered, and slightly weathered layers according to the degree of weathering. Figure 2 As shown, the geomechanical model consists of np = 124652 nodes and ne = 359216 elements.

[0099] For complex geomechanical models, normal constraints are applied based on their boundary coordinates. By traversing the x, y, and z coordinates of np = 124652 nodes, the minimum x-axis coordinate is obtained as xmin = 0, the maximum x-axis coordinate is xmax = 452, the minimum y-axis coordinate is ymin = 0, the maximum y-axis coordinate is ymax = 341, and the minimum z-axis coordinate is zmin = 110. With a tolerance of error = 0.01, normal constraints are applied to the left boundary (x = xmin), right boundary (x = xmax), front boundary (y = ymax), rear boundary (y = ymin), and bottom boundary (z = zmin) of the geomechanical model.

[0100] The following parameters were set for carbonaceous phyllite (strongly weathered): elastic modulus = 1.5 GPa, Poisson's ratio = 0.31, cohesion = 50 kPa, internal friction angle = 21.2°, tensile strength = 4.5 MPa; carbonaceous phyllite (weakly weathered): elastic modulus = 3.0 GPa, Poisson's ratio = 0.29, cohesion = 150 kPa, internal friction angle = 25.6°, tensile strength = 5.0 MPa; carbonaceous phyllite (slightly weathered): elastic modulus = 4.5 GPa, Poisson's ratio = 0.26, cohesion = 550 kPa, internal friction angle = 31.0°, tensile strength = 6.5 MPa; siliceous slate (strongly weathered) elastic modulus... Elastic modulus = 3 GPa, Poisson's ratio 0.29, cohesion = 250 kPa, internal friction angle = 21.8°, tensile strength = 3.8 MPa; Elastic modulus of siliceous slate (weakly weathered) = 6 GPa, Poisson's ratio 0.26, cohesion = 750 kPa, internal friction angle = 31.0°, tensile strength = 4.5 MPa; Elastic modulus of siliceous slate (slightly weathered) = 9 GPa, Poisson's ratio 0.23, cohesion = 850 kPa, internal friction angle = 35.0°, tensile strength = 5.0 MPa; Elastic modulus of mud = 22 MPa, Poisson's ratio 0.35, cohesion = 15 kPa, internal friction angle = 21.8°, tensile strength = 3.8 MPa. Friction angle = 14.0°, tensile strength = 40 kPa; foliated chlorite-bearing sodium schist has an elastic modulus of 28 MPa, Poisson's ratio of 0.32, cohesion of 18 kPa, internal friction angle of 21.8°, and tensile strength of 120 kPa; limestone has an elastic modulus of 0.2 GPa, Poisson's ratio of 0.27, cohesion of 100 kPa, internal friction angle of 21.8°, and tensile strength of 0.3 MPa; limestone (weakly and slightly weathered) has an elastic modulus of 0.8 GPa, Poisson's ratio of 0.28, cohesion of 200 kPa, internal friction angle of 31.0°, and tensile strength of 0.5 MPa. The elastic modulus of chlorite-albite schist (strongly weathered) is 2.0 GPa, Poisson's ratio is 0.29, cohesion is 0.8 MPa, internal friction angle is 25.6°, and tensile strength is 1.3 MPa; the elastic modulus of chlorite-albite schist (weakly weathered) is 5.0 GPa, Poisson's ratio is 0.28, cohesion is 1.7 MPa, internal friction angle is 31.0°, and tensile strength is 1.8 MPa; the elastic modulus of chlorite-albite schist (slightly weathered) is 8.5 GPa, Poisson's ratio is 0.26, cohesion is 4.2 MPa, internal friction angle is 34.0°, and tensile strength is 3.5 MPa; the shear expansion angle is taken as 0.0.

[0101] Gravity is applied in the negative z-axis direction, and the initial stress field is calculated to achieve equilibrium. The nodes np = 124652 are traversed, and their velocity and displacement values ​​are set to 0 to represent the initial geostress state.

[0102] (2) Perform excavation numerical simulation using the geomechanical model, calculate the geostress balance, and retain the strain field results. Figure 3 The image shown is a contour map of the shear strain increment of the geomechanical model.

[0103] like Figure 4 As shown, the average normal direction of the slope surface of the geomechanical model is set as the vertical direction, and the unit vector in the vertical direction is denoted as (v1, v2, v3); vertical lines are arranged at equal intervals along the X-axis and Y-axis, extending from the slope surface to the boundary of the geomechanical model, and the spacing between the vertical lines is smaller than the average size of the element.

[0104] Create a two-dimensional array (xx, yy) to represent each vertical line, where xx is the X-coordinate of the intersection of the vertical line and the bottom surface of the model, yy is the Y-coordinate of the intersection of the vertical line and the bottom surface of the model, and the spacing between adjacent vertical lines is dd. Arrange measuring points at equal intervals from bottom to top inside the geomechanical model along each vertical line, with the coordinates of the measuring points being (q*v1+xx, q*v2+yy, q*v3+zz), where q is the measuring point number on each vertical line, and zz is the Z-coordinate of the intersection of the vertical line and the bottom surface of the model.

[0105] Using the built-in FISH function to count the shear strain increments at each measuring point on the vertical line, ① find the measuring point number q corresponding to the maximum value. max The location of the measuring point with the maximum shear strain increment on the corresponding vertical line is (q) max *v1+xx,q max *v2+yy,q max *v3+zz); ② Find the maximum value q of the measuring point number. max2 The corresponding measuring point on the perpendicular line is (q) max2 *v1+xx,q max2 *v2+yy,q max2 *v3+zz).

[0106] (3) The locations of the measuring points for the maximum shear strain increment are grouped according to whether they are located on the slope surface, such as... Figure 5 As shown, the measuring points located on the slope surface are divided into the slope surface group and the measuring points located inside the slope are divided into the slope inside group. The slope surface group is named "Slope Surface" and the slope inside group is named "Slope Inside".

[0107] (4) Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of the measuring points located on the slope. Group the measuring points according to the yy values, and the yy values ​​of the measuring points in each group are the same. Iterate through all the measuring points in each group in ascending order of xx values. First, assign the first measuring point to segment x1. If the vertical line of the next measuring point is adjacent to the vertical line of the previous measuring point, then the next measuring point is also assigned to segment x1, and so on. If the vertical line of the next measuring point is not adjacent to the vertical line of the previous measuring point, then the next measuring point is assigned to segment x2, and so on.

[0108] Calculate the shear strain increments at each measuring point within each row and segment, and select the measuring point corresponding to the maximum shear strain increment, such as... Figure 6 As shown by the red dot in the middle.

[0109] Then, the two-dimensional array (xx, yy) corresponding to the perpendicular lines of the measuring points located on the slope is calculated. The measuring points are grouped according to the xx value, and the xx values ​​of the measuring points in each group are the same. In turn, all measuring points in each group are traversed in ascending order of yy value. The first measuring point is assigned to segment y1. If the perpendicular line of the next measuring point is adjacent to the perpendicular line of the previous measuring point, the next measuring point is also assigned to segment y1, and so on. If the perpendicular line of the next measuring point is not adjacent to the perpendicular line of the previous measuring point, the next measuring point is assigned to segment y2, and so on.

[0110] Calculate the shear strain increments at each measuring point within each column and segment, and select the measuring point corresponding to the maximum shear strain increment, such as... Figure 6 As shown by the red dot in the middle.

[0111] (5) In each selected segment, take any one of the measurement points with the maximum value of shear strain increment and classify it into the inner circle group g1. With g1 as the center, draw a circle with a radius threshold δ = 3 and classify the ungrouped measurement points in the circle into the inner circle group g1.

[0112] Traverse the measurement points that are classified into the inner group g1. Draw a circle with a radius threshold δ = 3 centered on the circle. Classify the ungrouped measurement points in the circle into the inner group g1. Continue in this manner until all measurement points in the circle are classified into the inner group g1.

[0113] Then, randomly select one point from the ungrouped measurement points and assign it to the inner circle group g2. Repeat the first two steps to assign all ungrouped measurement points within the circle to the inner circle group g2, and so on, until the inner circle group gN is reached, at which point all measurement points have been grouped. Figure 7 As shown.

[0114] (6) For each inner circle group, select the points with the largest and second largest shear strain increment values ​​in each row and column along the X and Y axes, respectively, and discard the remaining points in that inner circle group; if there is only one point in a row or column, select only the point with the largest shear strain increment value in that row or column. Then repeat step (5) to update the grouping of all points and obtain the updated inner circle groups, such as... Figure 8 As shown.

[0115] (7) Figure 8 As shown, based on the morphology of the measuring points within each updated inner circle group, the inner circle groups most likely to form closed loops are groups numbered 1, 5, 7, 9, and 11. Therefore, groups numbered 1, 5, 7, 9, and 11 are selected to construct the slip surface boundary line exposed on the slope, while the remaining groups are discarded. The group numbers are then updated to 1, 2, 3, 4, and 5. Following the boundary line connection order, the closest pair of points in adjacent groups is found and set as the start and end points, as shown below. Figure 9 As shown, the names are set to "start" and "end" respectively.

[0116] (8) Select the measurement point with the name "start" in the group as the first measurement point, draw a circle with a radius threshold of δ = 3 centered on it, classify all the points inside the circle into the mapping group m1 of the first measurement point, and sort them in ascending order of the distance from the first measurement point.

[0117] Then select the second measurement point with the smallest distance in the mapping group m1, draw a circle with a radius threshold of δ = 3 centered on it, classify all the measurement points inside the circle except the first measurement point into the mapping group m2 of the second measurement point, and sort them in ascending order of the distance from the second measurement point. And so on, until the point with the name "end" in the group is selected, and then connect all the selected points in sequence, as Figure 10 shown.

[0118] When selecting the third measurement point, that is, the measurement point with the smallest distance in the mapping group m2, compare the distance s2 from the second measurement point to this measurement point with the distance s1 from the first measurement point to this measurement point: If s2 < s1, then select this measurement point as the third measurement point; If s2 ≥ s1, then select the measurement point with the second smallest distance from this measurement point in the mapping group m2, and compare the distance s2 from the second measurement point to this measurement point with the distance s1 from the first measurement point to this measurement point again, until the measurement point corresponding to s2 < s1 is found as the third measurement point, and so on, to determine the fourth, fifth... until the measurement point with the name "end".

[0119] If, after traversing all the measurement points in the mapping group mB, the distance s2 from the second measurement point to this measurement point is not less than the distance s1 from the second measurement point to this measurement point, then return to the mapping group m1 of the first measurement point, select the point with the second smallest distance in the mapping group m1 as the second measurement point, and so on.

[0120] Finally, select the starting point "start" and the ending point "end" of adjacent groups, connect them, and form a preliminary slip surface boundary line of the slope outcrop, as Figure 10 shown.

[0121] (9) According to the preliminary slip surface boundary line of the slope outcrop, traverse all the measurement points in the connection order, and count the two-dimensional arrays corresponding to the perpendicular lines where two adjacent measurement points P and Q are located. The array corresponding to the measurement point P is denoted as (Rp, Cp), and the array corresponding to the measurement point Q is denoted as (Rq, Cq). If (Rq - Rp) / dd > 1, then boundary measurement points need to be supplemented: Rs ∈ [Rp + dd, Rq - dd], with a step spacing of dd, and calculate Cs according to the following formula:

[0122]

[0123] Determine the parameter q corresponding to the perpendicular line according to the two-dimensional array (Rs, Cs) max2 , and the position of the supplemented boundary measurement point is (q max2 * v1 + Rs, q max2*v2+Cs,q max2 *v3+zz). Update the connection order after adding all boundary measurement points.

[0124] Then, traverse the supplemented boundary measurement points in the connection order, and count the two-dimensional arrays corresponding to the perpendicular lines of adjacent measurement points M and N. The array corresponding to measurement point M is denoted as (Rm, Cm), and the array corresponding to measurement point N is denoted as (Rn, Cn). If (Cn-Cm) / dd>1, then boundary measurement points need to be supplemented: Ct∈[Cm+dd, Cn-dd], and the step size interval dd is calculated as follows: Rt is calculated as follows:

[0125]

[0126] Determine the parameter q of the corresponding perpendicular line based on the two-dimensional array (Rt, Ct). max2 The location of the supplementary boundary measurement point is (q max2 *v1+Rt,q max2 *v2+Ct,q max2 *v3+zz). After completing all boundary measurement points, update the connection order to obtain the supplemented slip surface boundary line exposed on the slope surface. Figure 11 The image shows the exposed slip surface boundary line of the slope surface after the missing measuring points have been supplemented.

[0127] (10) Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point on the exposed slip surface boundary line of the slope surface after statistical analysis. Divide the data into multiple groups by grouping measuring points with the same yy value. Then, iterate through the measuring points in each group in ascending order of xx value as the current measuring point. Based on the fact that the current measuring point belongs to the current third segment x... k If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the current third segment x. k If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next third segment x. k+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the third segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the third segments x1, x2, x3...x L In the first group that is traversed, the smallest measurement point of xx value belongs to the third segment x1, k is the number of the third segment, k∈1,2,3……L, and L is the number of the third segments;

[0128] Calculate the x values ​​of the test points within each third segment, and create a two-dimensional array (x1, x2) and a four-dimensional array (y0, y1, y2, y3) to represent each third segment. x1 is the minimum x value of the test points within each third segment, x2 is the maximum x value of the test points within each third segment, y1 is the yy value of the test point corresponding to the x1 value within each third segment, y2 is the yy value of the test point corresponding to the x2 value within each third segment, y0 is the yy value of the test point preceding the y1 value according to the boundary line connection order, and y3 is the yy value of the test point following the y2 value according to the boundary line connection order.

[0129] Create a feature variable V to represent each third segment, and obtain the values ​​of y0, y1, y2 and y3 for each third segment. If (y0-y1)×(y3-y2)<0, set the V value of the third segment to 1; otherwise, set the V value to 0.

[0130] For each measuring point on the exposed slip surface boundary line of the supplemented slope, the yy values ​​of the two-dimensional array (xx, yy) corresponding to the vertical line are traversed in ascending order, and an initial feature value T(yy) = 0 is assigned to each yy value. For each yy value, the xx values ​​of the two-dimensional array (xx, yy) of all vertical lines in the geomechanical model are traversed in ascending order, and compared with the x1 values ​​of all third segments within the corresponding group. When x1 = xx exists, the V value of the corresponding third segment is obtained, and T(yy) is updated to the sum of its current value and the V value. If the updated T(yy) is odd, the x2 value of the corresponding third segment is obtained. The value is updated by setting the starting position of the xx value ss = x2 + dd. Starting from x2 + dd, the xx values ​​are traversed in ascending order and compared with the x1 values ​​of all third segments in the corresponding group. When x1 = xx exists, the V value of the corresponding third segment is obtained, and T(yy) is updated to the sum of its current value and V. If the updated T(yy) is even, the two-dimensional array (xx, yy) in the interval [ss, x1 - dd] is divided into the bottom sliding surface group, and then the x2 value of the corresponding third segment is obtained. Starting from x2 + dd, the xx values ​​are traversed in ascending order, and so on, until the xx values ​​are traversed.

[0131] Iterate through the two-dimensional array (xx, yy) divided into the bottom sliding surface group, and determine the parameter q of the corresponding perpendicular line based on each two-dimensional array. max Establish a measuring point, the location of which is (q). max *v1+xx,q max *v2+yy,q max *v3+zz), and assign the measuring points to the bottom sliding surface group, such as Figure 12 As shown.

[0132] Using the quadrilaterals formed by the pairwise adjacent xx and yy values ​​of the two-dimensional array (xx, yy) of all vertical lines in the geomechanical model, the measuring points of the bottom slip surface group and the measuring points of the slip surface boundary line exposed on the supplemented slope surface are grouped together. All groups are traversed; if a group contains 3 or 4 measuring points, they are connected sequentially to form triangles or quadrilaterals. The surface formed by all triangles and quadrilaterals is the three-dimensional potential slip surface. Figure 13 As shown.

[0133] (11) Traverse all polygons that constitute the three-dimensional potential sliding surface, and use the built-in FISH function to calculate the area A of the k-th polygon. k And using the stress values ​​of the m vertices of the polygon in the global coordinate system xyz, calculate the average stress tensor of the k-th polygon:

[0134]

[0135] In the formula, S ki Let S be the stress tensor of the i-th node of the k-th polygon. k Composed of 6 stress components {σ xk ,σ yk ,σ zk ,τ xyk ,τ xzk ,τ yzk} constitute, σ xk σ yk σ zk τ represents the normal stress in the x, y, and z directions of the k-th polygon, respectively. xyk τ xzk τ yzk These are the shear stresses on the xy, xz, and yz planes of the k-th polygon, respectively.

[0136] For the k-th polygon, the z′ axis of the local coordinate system is defined as its normal direction, the x′ axis is defined as the average displacement direction of the m vertices, and the y′ axis is defined as the cross product direction of z′ and x′. The mean stress tensor is transformed from the global coordinate system to the local coordinate system using the following formula:

[0137] S′ k =LS k L T

[0138] In the formula, (l1m1 n1) is the direction cosine of the x′ axis of the local coordinate system in the global coordinate system, (l2m2 n2) is the direction cosine of the y′ axis of the local coordinate system in the global coordinate system, and (l3 m3 n3) is the direction cosine of the z′ axis of the local coordinate system in the global coordinate system.

[0139] Take S′ k The stress component σ′ zk For normal stress, τ′ xzk For tangential stress, the anti-slip force and sliding force of the polygon are accumulated, and the stability safety factor is calculated using the following formula:

[0140]

[0141] In the formula, F s The stability safety factor is n, where n is the total number of polygons, and c is the total number of polygons. k , Let A be the cohesion and internal friction angle of the k-th polygon. k Let be the area of ​​the k-th polygon.

[0142] The calculated stability safety factor F of the three-dimensional potential sliding surface s =1.32, such as Figure 13 As shown.

[0143] This application utilizes an existing continuous numerical simulation platform to establish a geomechanical model that considers complex geological formations. The strain field generated under varying external loads is calculated through numerical simulation. The model is then discretized into a series of vertical lines, and the location of the maximum shear strain increment on these lines is searched. A self-developed algorithm is used to search for the complex asymmetric three-dimensional potential sliding surface of the soil and rock mass and calculate the safety factor. The entire calculation method is fast, reasonable, and logically sound, accurately reflecting the three-dimensional potential sliding surface and safety factor of the slope. This provides a more realistic assessment of slope stability and offers significant insights for searching for sliding surfaces, calculating safety factors, and determining remediation schemes for soil-rock boundary slopes.

[0144] 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.

[0145] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. 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 patent application should be determined by the appended claims.

Claims

1. A method for determining an asymmetric three-dimensional potential sliding surface using the maximum shear strain increment, characterized in that, The method includes: Step S1: Based on the FLAC3D continuous numerical simulation platform, establish a geomechanical model that considers complex strata distribution and consists of nodes and elements, and initialize it. The initialization includes: assigning deformation and strength parameters, applying self-weight external load and boundary conditions to calculate the initial stress balance, and clearing displacement and velocity to zero. Step S2: Perform excavation numerical simulation on the geomechanical model until the stress equilibrium state is reached, retain the strain field results at the stress equilibrium state, arrange equally spaced vertical lines on the geomechanical model, and arrange measuring points at equally spaced intervals from bottom to top along the vertical lines. Based on the strain field results at the stress equilibrium state, calculate the location of the measuring point with the maximum shear strain increment on each vertical line. Step S3: Group the measurement points based on whether they are located on the slope surface, and assign the measurement points located on the slope surface to the slope surface group and the measurement points located inside the slope to the slope inside group. Step S4: Divide all the measuring points in the slope table group into segments along the X-axis and Y-axis according to continuity, determine the set of measuring points with the maximum value of shear strain increment, and discard the remaining measuring points in the slope table group; Step S5: For the measurement points in the set of measurement points with the maximum value of shear strain increment, use a clustering algorithm to group any two measurement points whose distance between them is less than the threshold δ into a group, and determine the inner circle group to which all measurement points belong through iteration to obtain each inner circle group; Step S6: Sequentially add the measuring points with the largest and second largest shear strain increment values ​​in each row and column of the bottom surface of the geomechanical model along the X-axis and Y-axis to the retained measuring point set, and delete the remaining measuring points in the inner circle group; if there is only one point in a row or column, only the point with the largest shear strain increment value in that row or column is retained and added to the retained measuring point set. Step S7: Then perform step S5 again on the measured points in the retained set of measured points to update the grouping of all points and obtain the updated inner circle groups; Step S8: Based on the shape of the measuring points in each updated inner circle group, find the inner circle group combination that is most likely to form a closed loop. Determine this inner circle group combination as the inner circle group combination that constitutes the boundary line of the exposed slip surface of the slope. Discard the other inner circle groups. According to the boundary line connection order, in the measuring points of two adjacent inner circle groups, find the point pair that is closest to a measuring point in one inner circle group and a measuring point in another inner circle group. In the point pair, set the previous measuring point as the starting point and the next measuring point as the ending point according to the boundary line connection order. Step S9: Using an improved maze algorithm, the measuring points of each circle group are selectively connected end-to-end according to the starting and ending points. Then, the starting and ending points of adjacent groups are connected according to the connection order, and unconnected measuring points are discarded to form the preliminary slip surface boundary line of the slope surface. Step S10: Sort the points connected on the preliminary slip surface boundary line of the slope surface according to the connection order, and determine whether the two adjacent measuring points are continuous along the X-axis and Y-axis. If they are not continuous, supplement the missing boundary measuring points according to the direction of the line connecting the two measuring points to obtain the supplemented slip surface boundary line of the slope surface. Step S11: Search for all measurement points with the maximum shear strain increment within the slip surface boundary line exposed by the supplemented slope surface, group them together with the measurement points of the slip surface boundary line exposed by the supplemented slope surface, and then connect them sequentially to obtain a three-dimensional potential slip surface. Step S12: Traverse all polygons that constitute the three-dimensional potential sliding surface, calculate the average stress tensor of each polygon, and accumulate the stability safety factor of the potential sliding surface.

2. The method according to claim 1, characterized in that, Step S2 includes: Step S21: Perform excavation numerical simulation on the geomechanical model until the stress equilibrium state is reached, and retain the strain field results at the stress equilibrium state. Step S22: Set the average normal direction of the slope surface of the geomechanical model as the vertical direction, and denote the unit vector in the vertical direction as (v1, v2, v3); arrange vertical lines at equal intervals along the X-axis and Y-axis, with the vertical lines extending from the slope surface to the boundary of the geomechanical model, and the spacing between the vertical lines is smaller than the average size of the unit. Step S23: Create a two-dimensional array (xx, yy) to represent each vertical line, where xx is the X-coordinate of the intersection point of the vertical line and the bottom surface of the geomechanical model, yy is the Y-coordinate of the intersection point of the vertical line and the bottom surface of the geomechanical model, and the spacing between adjacent vertical lines is dd; arrange measuring points at equal intervals from bottom to top inside the geomechanical model along each vertical line, with the coordinates of the measuring points being (q*v1+xx, q*v2+yy, q*v3+zz), where q is the measuring point number on each vertical line, and zz is the Z-coordinate of the intersection point of the vertical line and the bottom surface of the geomechanical model; Step S24: Based on the strain field results under the stress equilibrium state, use the built-in FISH function to statistically analyze the shear strain increments at each measuring point on the vertical line, and find the measuring point number q corresponding to the maximum shear strain increment on the vertical line. max The location of the measuring point with the maximum shear strain increment on the corresponding vertical line is (q max *v1+xx,q max *v2+yy,q max *v3+zz); Find the maximum value q of the measuring point number on the vertical line. max2 The corresponding measuring point on the perpendicular line is (q) max2 *v1+xx,q max2 *v2+yy,q max2 *v3+zz).

3. The method according to claim 2, characterized in that, Step 3 includes: Step S31: Traverse each perpendicular line and compare the measurement point number q on each perpendicular line. max and q max2 The value of q, if max =q max2 If the maximum shear strain increment of the vertical line is measured at the slope surface, then the maximum shear strain increment of the vertical line is measured at the slope surface; otherwise, the maximum shear strain increment of the vertical line is measured at the slope interior. Step S32: Assign the measuring points located on the slope surface to the slope surface group and the measuring points located inside the slope to the slope inside group.

4. The method according to claim 1, characterized in that, Step S4 includes: Step S41: Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point in the slope table group. Divide the group into multiple groups by grouping measuring points with the same yy value. Iterate through the measuring points within each group in ascending order of xx value as the current measuring point, based on whether the current measuring point belongs to the current first segment x. i If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the current first segment x. i If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next first segment x. i+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the first segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the first segments x1, x2, x3...x M In the first group that is traversed, the smallest measurement point of xx value belongs to the first segment x1, i is the first segment number, i∈1,2,3……M, and M is the number of the first segments; Step S42: Calculate the shear strain increment of each measuring point in the first segment, and select the measuring point corresponding to the maximum shear strain increment to add to the first measuring point set; Step S43: Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of the measuring points located on the slope surface. Divide the data into multiple groups, grouping measuring points with the same xx value as a group. Then, iterate through each group, selecting the measuring points within that group in ascending order of yy value as the current measuring point. The current measuring point is determined based on its belonging to the second segment y. j If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the second segment y. j If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next second segment y. j+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the second segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the second segments y1, y2, y3...y M In the first group that is traversed, the point with the smallest yy value belongs to the second segment y1, j is the number of the second segment, j∈1,2,3……Q, and Q is the number of the second segments; Step S44: Calculate the shear strain increment of each measuring point in the second segment, and select the measuring point corresponding to the maximum shear strain increment to add to the second measuring point set; Step S45: Take the union of the first set of measuring points and the second set of measuring points as the set of measuring points with the maximum value of shear strain increment, and discard all unselected measuring points in the slope table group.

5. The method according to claim 1, characterized in that, Step S5 includes: Step S51: First, randomly select one measurement point from the set of measurement points with the maximum value of shear strain increment as the current measurement point and assign it to the inner circle group g. n Draw a circle with δ as its center and radius as its threshold. Then, group the ungrouped measurement points within the circle that fall within the set of measurement points with the maximum shear strain increment into the inner circle group g. n And will be classified as the inner group g of the circle. n The measuring point is removed from the set of measuring points with the maximum value of shear strain increment, wherein the first randomly selected measuring point in the set of measuring points with the maximum value of shear strain increment is assigned to the inner circle group g1. Step S52: Next, from the inner circle group g n Select the measurement point that is not currently used as the current measurement point, and draw a circle with the circle as the center and the radius as the threshold δ. Then, group the measurement points that fall within the circle and are not grouped from the set of measurement points with the maximum shear strain increment into the group g within the circle. n And will be classified as the inner group g of the circle. n The measuring points are removed from the set of measuring points with the maximum shear strain increment, and so on, until the inner circle group g. n There are no measurement points that have not been selected as the current measurement point; Step S53: Then, randomly select one measurement point from the set of measurement points with the maximum value of the shear strain increment and assign it to the inner circle group g. n+1 Repeat steps S51 and S52 to determine which group belongs to the inner circle group g. n+1 The measuring points will be assigned to the inner circle group g. n+1 The measuring point is removed from the set of measuring points with the maximum shear strain increment, and so on. All measuring points in the set of measuring points with the maximum shear strain increment are grouped to obtain inner circle groups g1, g2...gN, where gN is the last inner circle group.

6. The method according to claim 1, characterized in that, Step S9 includes: Step S91: Select the starting point in a circle group as the first measurement point in the circle group, draw a circle with the starting point as the center and the radius as the threshold δ, and classify all measurement points in the circle group that fall into the circle into the mapping group m1 of the first measurement point, and sort them in ascending order of distance between each measurement point in the mapping group m1 and the first measurement point. Step S92: Select the measuring point a with the smallest distance in the mapping group m1 as the second measuring point. Draw a circle with radius δ centered on a point a. Assign all measuring points in the group that fall within this circle and are not the first measuring point to the mapping group m of measuring point a. a and according to mapping group m a The distances between each measuring point and measuring point a are sorted from smallest to largest; Step S93: Select mapping group m a Select the measurement point b with the smallest internal distance. Compare the distance s2 from the second measurement point to the measurement point b with the distance s1 from the first measurement point to the measurement point b. If s2 < s1, then select the measurement point b as the third measurement point; if s2 ≥ s1, then select the mapping group m a Select the measurement point b with the second smallest internal distance as the measurement point b. Compare the distance s2 from the second measurement point to the measurement point b with the distance s1 from the first measurement point to the measurement point b again. If s2 < s1, then select the measurement point b as the third measurement point; if s2 ≥ s1, then select the mapping group m a Select the measurement point b with the third smallest internal distance as the measurement point b, and follow the mapping group m a Judge in ascending order of the internal distance until the measurement point b corresponding to s2 < s1 is found as the third measurement point. Among them, if all the measurement points in the mapping group m a are traversed and the distance s2 is not less than the distance s1, then return to the mapping group m1 of the first measurement point, select the measurement point with the second smallest internal distance in the mapping group m1 as the second measurement point, and so on; Step S94: Based on the distance within the mapping group of the previous measuring point, select the next measuring point of the previous measuring point in ascending order, until the end point within the circle group is selected, and then connect the selected measuring points in sequence. Step S95: Select the starting and ending points of adjacent inner circles and connect them to form the preliminary slip surface boundary line of the slope surface.

7. The method according to claim 1, characterized in that, Step S10 includes: Step S101: Based on the preliminary exposed slip surface boundary line of the slope surface, traverse all measuring points on the preliminary exposed slip surface boundary line in the order of boundary line connection. Count the two-dimensional arrays corresponding to the perpendicular lines of pairwise adjacent measuring points P and Q. The array corresponding to measuring point P is denoted as (Rp, Cp), and the array corresponding to measuring point Q is denoted as (Rq, Cq). If (Rq-Rp) / dd>1, then the array corresponding to the boundary measuring points to be supplemented is (Rs, Cs), where Rs∈[Rp+dd, Rq-dd]. The step size interval between the boundary measuring points to be supplemented is dd. The expression for calculating Cs is: Determine the parameter q of the corresponding perpendicular line based on the two-dimensional array (Rs, Cs). max2 The location of the supplementary boundary measurement point is (q max2 *v1+Rs,q max2 *v2+Cs,q max2 *v3+zz), update the connection order after adding all boundary measurement points; Step S102: Traverse the supplemented boundary measurement points according to the connection order, and count the two-dimensional arrays corresponding to the perpendicular lines of adjacent measurement points M and N. The array corresponding to measurement point M is denoted as (Rm, Cm), and the array corresponding to measurement point N is denoted as (Rn, Cn). If (Cn-Cm) / dd>1, then the array corresponding to the boundary measurement points to be supplemented is (Rt, Ct), where Ct∈[Cm+dd, Cn-dd]. The step size distance between the boundary measurement points to be supplemented is dd. The expression for calculating Rt is: Determine the parameter q of the corresponding perpendicular line based on the two-dimensional array (Rt, Ct). max2 The location of the supplementary boundary measurement point is (q max2 *v1+Rt,q max2 *v2+Ct,q max2 *v3+zz), after supplementing all boundary measurement points, update the connection order to obtain the supplemented slip surface boundary line exposed on the slope surface.

8. The method according to claim 1, characterized in that, Step S11 includes: Step S111: Calculate the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point on the exposed slip surface boundary line of the supplemented slope surface. Divide the area into multiple groups, grouping measuring points with the same yy value. Iterate through each group in ascending order of xx value, using the measuring points within that group as the current measuring point. The current measuring point is determined based on its belonging to the current third segment x. k If the perpendicular line to the next measuring point is adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the current third segment x. k If the perpendicular line to the next measuring point is not adjacent to the perpendicular line to the current measuring point, then the next measuring point is assigned to the next third segment x. k+1 This process continues, selecting the next measurement point as the current measurement point, and then determining the third segment to which each measurement point belongs, until all measurement points in all groups have been selected, thus obtaining the third segments x1, x2, x3...x L In the first group that is traversed, the smallest measurement point of xx value belongs to the third segment x1, k is the number of the third segment, k∈1,2,3……L, and L is the number of the third segments; Step S112: Calculate the xx values ​​of the test points within each third segment, and create a two-dimensional array (x1, x2) and a four-dimensional array (y0, y1, y2, y3) to represent each third segment. x1 is the minimum xx value of the test points within each third segment, x2 is the maximum xx value of the test points within each third segment, y1 is the yy value of the test point corresponding to the x1 value within each third segment, y2 is the yy value of the test point corresponding to the x2 value within each third segment, y0 is the yy value of the test point preceding the test point corresponding to the y1 value according to the boundary line connection order, and y3 is the yy value of the test point following the test point corresponding to the y2 value according to the boundary line connection order. Step S113: Create a feature variable V to represent each of the third segments, and obtain the values ​​of y0, y1, y2 and y3 of each third segment. If (y0-y1)×(y3-y2)<0, set the V value of the third segment to 1; otherwise, set the V value to 0. Step S114: For the yy values ​​of the two-dimensional array (xx, yy) corresponding to the vertical lines of each measuring point of the supplemented slope surface exposed slip surface boundary line, traverse the yy values ​​in ascending order and assign an initial feature value T(yy) = 0 to each yy value; for each yy value, traverse the xx values ​​of the two-dimensional array (xx, yy) of all vertical lines of the geomechanical model in ascending order and compare them with the x1 values ​​of all third segments within the corresponding group. When x1 = xx exists, obtain the V value of the corresponding third segment and update T(yy) to the sum of its current value and V value; if the updated T(yy) is odd, obtain the corresponding third segment. For each segment x2 value, update the starting position xx value ss = x2 + dd. Starting from x2 + dd, continue traversing the xx values ​​in ascending order and compare them with the x1 values ​​of all third segments in the corresponding group. When x1 = xx exists, obtain the V value of the corresponding third segment and update T(yy) to the sum of its current value and V. If the updated T(yy) is even, divide the two-dimensional array (xx, yy) in the interval [ss, x1 - dd] into the bottom sliding surface group, then obtain the corresponding third segment x2 value, and continue traversing the xx values ​​in ascending order starting from x2 + dd, and so on, until all xx values ​​have been traversed. Step S115: Traverse the two-dimensional array (xx, yy) divided into the bottom sliding surface group, and determine the parameter q of the corresponding perpendicular line according to each two-dimensional array. max Establish a measuring point, the location of which is (q). max *v1+xx,q max *v2+yy,q max *v3+zz), and the measuring points are assigned to the bottom sliding surface group; Using the quadrilateral formed by the pairwise adjacent xx and yy values ​​of the two-dimensional array (xx, yy) of all the vertical lines of the geomechanical model, the measuring points of the bottom slip surface group and the measuring points of the slip surface boundary line exposed by the supplemented slope surface are grouped into groups. All groups are traversed. If a group contains 3 or 4 measuring points, they are connected in sequence to form a triangle or quadrilateral. The surface formed by all the triangles and quadrilaterals is the three-dimensional potential slip surface.

9. The method according to claim 1, characterized in that, Step S12 includes: Step S121: Traverse all polygons constituting the three-dimensional potential sliding surface, and determine the average stress tensor of each polygon using the stress values ​​of the m vertices on each polygon in the global coordinate system xyz. The average stress tensor of the k-th polygon is: In the formula, S ki Let S be the stress tensor of the i-th vertex of the k-th polygon. k Composed of 6 stress components {σ xk ,σ yk ,σ zk ,τ xyk ,τ xzk ,τ yzk } constitute, σ xk σ yk σ zk τ represents the normal stress in the x, y, and z directions of the k-th polygon, respectively. xyk τ xzk τ yzk These are the shear stresses on the xy, xz, and yz planes of the k-th polygon, respectively. Step S122: For the k-th polygon, define the z′ axis of the local coordinate system as its normal direction, the x′ axis of the local coordinate system as the average displacement direction of the m vertices, and the y′ axis of the local coordinate system as the cross product direction of z′ and x′. Use the formula S′ k =LS k L T The mean stress tensor is transformed from the global coordinate system to the local coordinate system, where... (l1m1 n1) is the direction cosine of the x′ axis of the local coordinate system in the global coordinate system, (l2 m2 n2) is the direction cosine of the y′ axis of the local coordinate system in the global coordinate system, and (l3 m3 n3) is the direction cosine of the z′ axis of the local coordinate system in the global coordinate system. Step S123: Take S′ k The stress component {σ′ zk For normal stress, τ′ xzk For tangential stress, the anti-slip force and sliding force of all polygons are accumulated to determine the stability safety factor of the potential sliding surface. The stability safety factor of the potential sliding surface is: In the formula, F s The stability safety factor is n, where n is the total number of polygons, and c is the total number of polygons. k , Let A be the cohesion and internal friction angle of the k-th polygon. k Let be the area of ​​the k-th polygon.

Images