Slope stability safety factor calculation method fusing finite element method and limit equilibrium method
By integrating the finite element method and the limit equilibrium method, and combining the coupled calculation of unsaturated seepage stress in slopes, an index mode is established and integrals are performed along the slip zone path. This solves the problem of multi-physics field influence in slope stability calculation in existing technologies, and realizes accurate evaluation and reliable calculation of slope stability variation law.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- CHINA THREE GORGES UNIV
- Filing Date
- 2023-02-09
- Publication Date
- 2026-07-03
AI Technical Summary
Existing technologies, when calculating the stability safety factor of reservoir landslides, cannot effectively consider the slope stability variation law under the coupling of multiple factors and multiple physical fields, especially the influence of seepage field and stress deformation. Furthermore, the finite element method and limit equilibrium method have shortcomings in the reliability of calculation results.
By integrating the finite element method and the limit equilibrium method, and extracting the results of coupled finite element calculations of unsaturated seepage stress on the slope, an indexing mode of line segment set and finite element mesh is established. The normal force and shear force are solved by integral along the slip zone path, and the total anti-sliding force and total sliding force are calculated to finally obtain the stability safety factor.
It enables a reasonable evaluation of the slope stability variation law under the coupled action of multiple factors and multiple physical fields, and provides accurate calculation of the total sliding force and total anti-sliding force on the sliding surface. The calculation results are reliable and have clear physical meaning.
Smart Images

Figure CN116151072B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of hydraulic and geotechnical engineering, specifically to a method for calculating the slope stability safety factor that integrates the finite element method and the limit equilibrium method. Background Technology
[0002] Reservoir landslides pose a significant safety hazard and even property loss to people's lives and property. The stability of reservoir landslides is influenced by multiple factors, including rainfall and reservoir water level, as well as the coupling effects of multiple physical fields such as water vapor seepage and stress deformation. This often results in complex and varied patterns, making the calculation of the stability safety factor difficult.
[0003] The slope stability safety factor is one of the most important indicators for evaluating slope stability. The limit equilibrium method is widely used in calculating this factor. This method divides the slope into several blocks, performs static equilibrium analysis on the soil blocks, and then calculates the slope safety factor based on force equilibrium or moment equilibrium theory. However, the limit equilibrium method does not consider the influence of soil deformation and unsaturated water vapor pressure on slope stability. The finite element method (FEM) is typically used to analyze the multi-physics field variations within the slope. Stability calculation methods based on FEM are often used in evaluating landslide stability involving seepage fields or seepage-stress coupling. The FEM can consider the influence of slope stress and deformation under the action of water vapor migration. However, existing FEM-based methods such as the strength reduction method and the surcharge method require subjective judgment of the landslide instability state when calculating the stability safety factor, affecting the reliability of the calculation results and limiting the widespread application of FEM-based stability safety factor calculation methods. This method combines the advantages of finite element method and limit equilibrium method to calculate the stability safety factor. It can take into account the effects of landslide seepage and stress deformation, and has clear physical meaning and reliable calculation results. It is a slope stability safety factor calculation method with broad application prospects. Summary of the Invention
[0004] The purpose of this invention is to overcome the shortcomings of existing calculation methods. Based on the results of finite element calculations and the limit equilibrium method, this invention provides a method for calculating the slope stability safety factor that integrates the finite element method and the limit equilibrium method. This method helps to reasonably evaluate the slope stability variation law under multi-factor and multi-physical field coupling conditions.
[0005] To achieve the above objectives, this invention proposes a method for calculating the slope stability safety factor that integrates the finite element method and the limit equilibrium method, comprising the following steps:
[0006] 1) Provide the results of the coupled finite element analysis of unsaturated seepage stress on the slope, and extract the effective stress {σ} at each element node. x ′、σ y ′、τ x ′ y}, including the normal stress σ in the x-directionx Normal stress σ in the y direction y and shear stress τ x ′ y ;
[0007] 2) Obtain the set of line segments that intersect the sliding strip and the finite element mesh, and establish an indexing pattern between the line segment set and the finite element mesh;
[0008] 3) Solve for the normal and shear forces on each line segment set by integrating along the sliding track path;
[0009] 4) Calculate the total anti-slip force and total sliding force on the sliding belt;
[0010] 5) Calculate the stability safety factor.
[0011] Preferably, in step 2), the following steps are used to obtain the set of line segments intersecting the slide and the finite element mesh, and to establish an indexing pattern between the line segment set and the finite element mesh:
[0012] 2.1) Polylines are used to describe the sliding surface. Based on the position of the sliding surface, a coordinate system identical to the finite element mesh is established to obtain the coordinates of the control points of the sliding surface. The control points are connected in order from top to bottom to form polylines describing the sliding surface.
[0013] 2.2) Find the intersection points of the polyline and the finite element mesh segment by segment along the polyline numbering sequence;
[0014] 2.3) Sort the intersection points from top to bottom to obtain the line segment set. The total number of line segments in the line segment set is N. Number the line segments and intersection points to establish an index mode between the line segment set and the finite element mesh.
[0015] Preferably, in step 2.2), the intersection points of the polyline and the finite element mesh are found segment by segment using the coordinate rotation method. The specific steps are as follows:
[0016] a. Number each segment of the polyline as 1, 2, ..., i, ..., n;
[0017] b. For the i-th line segment, the coordinates of its endpoints are A(x) and A(x). i-1 y i-1 ) and B(x i y i Establish a new coordinate system with point A as the origin and A→B as the positive X-axis.
[0018] c. Calculate the rotation matrix from the original coordinate system to the new coordinate system. The rotation matrix R can be expressed as:
[0019]
[0020] In the formula, θ is the rotation angle, which can be expressed as:
[0021]
[0022]
[0023] d. Transform the finite element node coordinates to the new coordinate system (x’, y’) through translation and rotation. The transformation formula is
[0024]
[0025] e. Determine whether the four sides of the finite element mesh cell intersect with the line segment. In the new coordinate system, loop through the four sides of the finite element mesh cell. If the coordinates of the two endpoints of the j-th side are (x’ j-1 , y’ j-1 ) and (x’ j , y’ j ), and if y’ j-1 ×y’ j < 0, and the x’0 in the intersection point (x′0, 0) of the line segment and the x-axis satisfies x′0 < x′ i (x′ i is the abscissa of point B after rotation transformation), then the j-th side of the cell intersects with the i-th line segment of the polyline, and the intersection point coordinates are (x′0, 0);
[0026]
[0027] f. Calculate the coordinates of the intersection point in the original coordinate system The calculation formula is:
[0028]
[0029] g. If i = n, the calculation ends; if i < n, i = i + 1, and go to step b.
[0030] Preferably, in step 3), the integral method is used to solve the normal force and tangential force on the sliding zone line segment k. The calculation steps are as follows:
[0031] 3.1) According to the effective stress components of the finite element calculation results, find the normal stress and shear stress on the inclined section corresponding to the line segment k;
[0032] Normal stress:
[0033] Shear stress:
[0034]
[0035] In the formula, and are the coordinates of the two endpoints of the line segment k;
[0036] 3.2) Find the coordinates of the Gaussian integration points on the line segment k. The calculation formula is:
[0037]
[0038] In the formula, N 11 =N 22 =0.5773503, N 12 =N 21 =0.4226497;
[0039] 3.3) Calculate the weights of the Gaussian integration points of the line segment in the corresponding finite element mesh using the inverse distance weighted average method. The calculation formula is as follows:
[0040]
[0041] In the formula, d ij Let be the distance between the i-th Gaussian point (i = k-1, k) of the line segment and the j-th node of the finite element mesh element. The calculation formula is:
[0042]
[0043] In the formula, (x j ,y j Let be the coordinates of the j-th node. Let i be the coordinates of the i-th Gaussian point;
[0044] 3.4) Find the normal stress at the i-th Gaussian point of the line segment. and tangential stress
[0045] Normal stress:
[0046] Tangential stress:
[0047] In the formula, σ′ n,j Let τ′ be the normal stress at the j-th node of the element. n,j Let W′ be the shear stress at the j-th node of the element. ij Let be the weight of the i-th Gaussian point on the j-th node of the line segment;
[0048] 3.5) Use Gaussian integral to find the normal force F of the k-th segment of the sliding surface segment set. N,k and tangential force F S,k The calculation formula is:
[0049] Normal force:
[0050] Tangential force:
[0051] In the formula, W iW represents the weight of the Gaussian integral of the line segment unit. For the case of two Gaussian points, W... i =1.0.
[0052] Preferably, in step 4), the total anti-slip force of the sliding surface is the sum of the normal force of each segment and the tangent of the internal friction angle multiplied by the cohesive force. The total sliding force is the sum of the shear forces of each segment, and the calculation formula is as follows:
[0053] Total anti-skid force:
[0054] Total downward force:
[0055] In the formula, c k For line segment cohesion, Let l be the internal friction angle of the line segment. k The length of the line segment.
[0056] Preferably, in step 5), the stability safety factor is equal to the total anti-skid force divided by the total sliding force.
[0057]
[0058] The beneficial effects of this invention patent: Compared with the prior art, the advantages of this invention are as follows:
[0059] 1. It can obtain the variation law of slope stability safety factor under the coupling effect of multiple physical fields such as seepage field and stress field, considering the influence of multiple factors such as rainfall and reservoir water level. It can combine the advantages of seepage stress coupling finite element method and limit equilibrium method, and effectively overcome the problems and shortcomings of existing finite element method and limit equilibrium method in calculating safety factor.
[0060] 2. The slope stability calculation method of the present invention provides a practical approach. It adopts the method of precise integration along the path to achieve accurate calculation of the total sliding force and the total anti-sliding force on the sliding surface. The physical meaning is clear, the derivation is rigorous, and the calculation results are reliable.
[0061] 3. The method of this invention combines the advantages of the seepage stress coupling finite element method and the limit equilibrium method. By using precise integration along the path, it achieves accurate calculation of the sliding force and anti-sliding force of the sliding surface. The physical meaning is clear, the derivation is rigorous, and the calculation results are reliable. It helps to reasonably evaluate the slope stability variation law under the coupling of multiple factors and multiple physical fields. Attached Figure Description
[0062] Figure 1 This is a flowchart of a slope stability safety factor calculation method that integrates finite element calculation and limit equilibrium method according to the present invention.
[0063] Figure 2 The diagram shows the finite element mesh and slide position in the embodiment.
[0064] Figure 3 A schematic diagram of the indexing pattern of line segment sets and finite element meshes.
[0065] Figure 4 A schematic diagram of the Gaussian integration points of the line segment set and their relationship with the finite element mesh.
[0066] Figure 5 The slope stability safety factor calculation results are shown in the example. Detailed Implementation
[0067] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.
[0068] A method for calculating the slope stability safety factor that integrates the finite element method and the limit equilibrium method includes the following steps:
[0069] 1) Provide the results of the coupled finite element analysis of unsaturated seepage stress on the slope, and extract the effective stress {σ} at each element node. x ′、σ y ′、τ x ′ y}, including the normal stress σ in the x-direction x Normal stress σ in the y direction y and shear stress τ x ′ y ;
[0070] 2) Obtain the set of line segments that intersect the sliding strip and the finite element mesh, and establish an indexing pattern between the line segment set and the finite element mesh;
[0071] 3) Solve for the normal and shear forces on each line segment set by integrating along the sliding track path;
[0072] 4) Calculate the total anti-slip force and total sliding force on the sliding belt;
[0073] 5) Calculate the stability safety factor.
[0074] Preferably, in step 2), the following steps are used to obtain the set of line segments intersecting the slide and the finite element mesh, and to establish an indexing pattern between the line segment set and the finite element mesh:
[0075] 2.1) Polylines are used to describe the sliding surface. Based on the position of the sliding surface, a coordinate system identical to the finite element mesh is established to obtain the coordinates of the control points of the sliding surface. The control points are connected in order from top to bottom to form polylines describing the sliding surface.
[0076] 2.2) Find the intersection points of the polyline and the finite element mesh segment by segment along the polyline numbering sequence;
[0077] 2.3) Sort the intersection points from top to bottom to obtain the line segment set. The total number of line segments in the line segment set is N. Number the line segments and intersection points to establish an index mode between the line segment set and the finite element mesh.
[0078] Preferably, in step 2.2), the intersection points of the polyline and the finite element mesh are found segment by segment using the coordinate rotation method. The specific steps are as follows:
[0079] a. Number each segment of the polyline as 1, 2, ..., i, ..., n;
[0080] b. For the i-th line segment, the coordinates of its endpoints are A(x) and A(x). i-1 y i-1 ) and B(x i y i Establish a new coordinate system with point A as the origin and A→B as the positive X-axis.
[0081] c. Calculate the rotation matrix from the original coordinate system to the new coordinate system. The rotation matrix R can be expressed as:
[0082]
[0083] In the formula, θ is the rotation angle, which can be expressed as:
[0084]
[0085]
[0086] d. Transform the finite element node coordinates to a new coordinate system (x', y') through translation and rotation. The transformation formula is as follows:
[0087]
[0088] e. Determine whether the four edges of the finite element mesh element intersect with the line segment. In the new coordinate system, iterate through the four edges of the finite element mesh element. If the coordinates of the two endpoints of the j-th edge are (x'...) j -1, y' j-1 ) and (x' j y' j If y' j-1 ×y' j <0, and x′0 is the intersection point of the line segment and the x-axis (x′0, 0). <x i ′(x i (where x′ is the x-coordinate of point B after rotation transformation), then the j-th edge of the unit intersects the i-th segment of the polyline, and the coordinates of the intersection point are (x′0, 0).
[0089]
[0090] f. Find the coordinates of the intersection point in the original coordinate system. The calculation formula is as follows:
[0091]
[0092] g. If i = n, the calculation ends; if i < n, i = i + 1, and go to step b.
[0093] Preferably, in step 3), the integral method is used to solve the normal force and tangential force on the sliding zone segment k, and the calculation steps are as follows:
[0094] 3.1) According to the effective stress components of the finite element calculation results, find the normal stress and shear stress of the inclined section corresponding to segment k;
[0095] Normal stress:
[0096] Shear stress:
[0097]
[0098] In the formula, and are the coordinates of the two end points of segment k;
[0099] 3.2) Find the coordinates of the Gauss integration points of segment k. The calculation formula is as follows:
[0100]
[0101] In the formula, N 11 = N 22 = 0.5773503, N 12 = N 21 = 0.4226497;
[0102] 3.3) Find the weights of the Gauss integration points of the segment in the corresponding finite element mesh, and calculate them by the inverse distance weighted average method. The calculation formula is as follows:
[0103]
[0104] In the formula, d ij is the distance between the i-th Gauss point (i = k - 1, k) of the segment and the j-th node of the finite element mesh element. The calculation formula is as follows:
[0105]
[0106] The calculation formula is as follows: In the formula, (x j , y j ) are the coordinates of the j-th node, are the coordinates of the i-th Gauss point;
[0107] 3.4) Find the normal stress of the i-th Gauss point of the segment and tangential stress
[0108] Normal stress:
[0109] Tangential stress: In the formula, σ n ′ ,j Let τ be the normal stress at the j-th node of the element. n ′ ,j Let W be the shear stress at the j-th node of the element. ij ′ represents the weight of the i-th Gaussian point on the j-th node of the line segment;
[0110] 3.5) Use Gaussian integral to find the normal force F of the k-th segment of the sliding surface segment set. N,k and tangential force F S,k The calculation formula is:
[0111] Normal force:
[0112] Tangential force: In the formula, W i W represents the weight of the Gaussian integral of the line segment unit. For the case of two Gaussian points, W... i =1.0.
[0113] Preferably, in step 4), the total anti-slip force of the sliding surface is the sum of the normal force of each segment and the tangent of the internal friction angle multiplied by the cohesive force. The total sliding force is the sum of the shear forces of each segment, and the calculation formula is as follows:
[0114] Total anti-skid force:
[0115] Total downward force: In the formula, c k For line segment cohesion, Let l be the internal friction angle of the line segment. k The length of the line segment.
[0116] Preferably, in step 5), the stability safety factor is equal to the total anti-skid force divided by the total sliding force.
[0117]
[0118] The working process of this embodiment is as follows:
[0119] 1. The seepage field and stress field of a slope were simulated and calculated using the unsaturated seepage stress coupled finite element method to obtain the effective stress of each element node.
[0120] 2. Find the set of line segments that intersect the sliding zone with the finite element mesh.
[0121] ① Use polylines to describe the sliding band. In the same coordinate system as the finite element mesh, use polylines to draw the position of the sliding band in CAD, such as... Figure 2 As shown, a program written in Fortran reads the control coordinates of the slider position.
[0122] ② Find the intersection points of the polyline with the finite element mesh segment by segment along the polyline numbering sequence. Obtain the coordinates of the intersection points of the sliding polyline and the finite element mesh by using coordinate rotation.
[0123] ③ Sort the intersection points from top to bottom to obtain the line segment set, number the line segments and intersection points, and establish an index mode between the line segment set and the finite element mesh.
[0124] like Figure 3 The polyline is divided into multiple segments AB, BC, ..., and the intersection points 1, 2, 3, 4, 5, 6 of segment AB with the finite element mesh and the intersection points 6, 7, 8 of segment BC with the finite element mesh are obtained respectively. The segments and intersection points are numbered along the polyline direction, and the segment number, the two endpoint numbers of the segment, and the element number of the intersection with the finite element mesh are recorded. In this way, an index mode for the coordinates of the segment and the intersection point is established.
[0125] 3. Use the integration method along the path to find the normal and tangential forces of each sliding strip segment.
[0126] ① Use the formula for calculating the stress of an inclined section to determine the normal stress and shear stress of the inclined section corresponding to the line segment;
[0127] ② Find the coordinates of the Gaussian integral point of the line segment;
[0128] ③ The weights of the Gaussian integral points of the line segment in the corresponding finite element mesh are calculated using the near-point inverse distance weighting method;
[0129] ④ Find the normal and tangential stresses at the Gaussian integration point of the line segment;
[0130] ⑤ Use Gaussian integrals to calculate the normal and shear forces of the sliding strip segment set.
[0131] like Figure 4 The k-th segment of the sliding band intersects the i-th unit at points F1 and F2, with coordinates respectively. and Calculate the coordinates of the Gaussian point of the line segment as follows: and The weights of finite element nodes N1, N2, N3, and N4 at each Gaussian point are calculated using the near-point inverse distance weighting method. The normal and tangential stresses at the Gaussian integration points of the line segments are the sum of the normal and tangential stresses of the finite element nodes multiplied by the weights of each node. The normal and shear forces of the sliding strip line segment set are calculated using Gaussian integration.
[0132] 4. Calculate the total anti-slip force and the total sliding force. The total anti-slip force of the sliding surface is the sum of the product of the normal force of each segment and the tangent of the internal friction angle, plus the cohesion. The total sliding force is the sum of the tangential forces of each segment of the sliding strip.
[0133] 5. Calculate the stability safety factor, which is equal to the total anti-skid force divided by the total sliding force. For example... Figure 5 The slope stability safety factor and the change process of rainfall over time were calculated. It can be seen that the slope stability decreased significantly in the early stage of rainfall, and gradually stabilized as rainfall continued. The calculation results are reliable.
Claims
1. A method for calculating the slope stability safety factor that integrates the finite element method and the limit equilibrium method, characterized in that, Includes the following steps: 1) Provide the results of the coupled finite element analysis of unsaturated seepage stress on the slope, and extract the effective stress at each element node. , , }, including normal stress in the x-direction Normal stress in the y direction and shear stress ; 2) Obtain the set of line segments that intersect the sliding strip and the finite element mesh, and establish an indexing pattern between the line segment set and the finite element mesh; 3) Solve for the normal and shear forces on each line segment set by integrating along the sliding track path; 4) Calculate the total anti-slip force and total sliding force on the sliding belt; 5) Calculate the stability safety factor; In step 2), the following steps are used to obtain the set of line segments intersecting the slide and the finite element mesh, and to establish an indexing pattern between the line segment set and the finite element mesh: 2.1) Polylines are used to describe the sliding surface. Based on the position of the sliding surface, a coordinate system identical to the finite element mesh is established to obtain the coordinates of the control points of the sliding surface. The control points are connected in order from top to bottom to form polylines describing the sliding surface. 2.2) Find the intersection points of the polyline and the finite element mesh segment by segment along the polyline numbering sequence; 2.3) Sort the intersection points from top to bottom to obtain the line segment set. The total number of line segment sets is N. Number the line segments and intersection points and establish an index pattern between the line segment set and the finite element mesh. In step 2.2), the intersection points of the polyline and the finite element mesh are found segment by segment using the coordinate rotation method. The specific steps are as follows: a. Number each segment of the polyline as 1, 2, ... i , ..., n ; b. Regarding the first i Line segment, with endpoint coordinates A ( x i-1 , y i-1 ) and B ( x i , y i A new coordinate system is established with point A as the origin and A→B as the positive X-axis. c. Calculate the rotation matrix from the original coordinate system to the new coordinate system. The rotation matrix R can be expressed as: ; In the formula, Let be the rotation angle, which can be expressed as: ; d. Transform the finite element node coordinates to a new coordinate system through translation and rotation. x’ , y’ The transformation formula is: ; e. Determine whether the four edges of the finite element mesh element intersect the line segment. In the new coordinate system, iterate through the four edges of the finite element mesh element. If the first edge intersects the line segment... j The coordinates of the two endpoints of the strip are respectively ( x’ j-1 , y’ j-1 )and( x’ j , y’ j ),like y’ j-1 × y’ j < 0 And the line segment intersects the x-axis at ( ,0) , Let B be the x-coordinate after rotation transformation. Then the element's x-coordinate is... j Edge and polyline i The line segments intersect, and the coordinates of the intersection point are ( ,0); ; f. Find the coordinates of the intersection point in the original coordinate system. , The calculation formula is: ; g. If i = n The calculation is complete; like i < n , i = i +1, proceed to step b.
2. The method for calculating the slope stability safety factor by integrating the finite element method and the limit equilibrium method according to claim 1, characterized in that, In step 3), the integral method is used to solve for the sliding strip segment. k The calculation steps for the normal and tangential forces on the surface are as follows: 3.1) Based on the effective stress components from the finite element analysis results, determine the line segment. k The normal stress and shear stress corresponding to the inclined section; Normal stress: ; Shear stress: ; ; In the formula, ( , )and( , ) is a line segment k The coordinates of the two endpoints; 3.2) Find the line segment k The coordinates of the Gaussian integral point are calculated using the following formula: ; In the formula, , ; 3.3) Calculate the weights of the Gaussian integration points of the line segment in the corresponding finite element mesh using the inverse distance weighted average method. The calculation formula is as follows: ; In the formula, For line segment number 1 i The Gaussian point and the finite element mesh element j The distance between nodes, i = k-1 , k ; The calculation formula is: ; In the formula, ( , ) is the first j The coordinates of each node, ( , ) is the first i The coordinates of a Gaussian point; 3.4) Find the first segment of line segment 1 i Normal stress at a Gaussian point and tangential stress ; Normal stress: ; Tangential stress: ; In the formula, For unit number j Normal stress at each node, For unit number j Shear stress at each node, For line segment number 1 i The Gaussian point at the th... j The weights on each node; 3.5) Using Gaussian integrals to find the first set of sliding surface segments k Normal force of line segment and tangential force The calculation formula is: Normal force: ; Tangential force: ; In the formula, The weights are the Gaussian integrals of the line segment unit; for the case of two Gaussian points, .
3. The slope stability safety factor calculation method integrating the finite element method and the limit equilibrium method according to claim 1, characterized in that, in step 4), the total anti-slip force of the sliding surface is the product of the normal force of each segment and the tangent of the internal friction angle, and the total sliding force is the sum of the shear forces of each segment, calculated by the following formula: Total anti-skid force: ; Total downward force: ; In the formula, For line segment cohesion, The internal friction angle of the line segment. The length of the line segment.
4. The method for calculating the slope stability safety factor by integrating the finite element method and the limit equilibrium method according to claim 1, characterized in that, In step 5), the stability safety factor equals the total anti-skid force divided by the total sliding force. 。