A method for analyzing the stability of soil slopes based on the exact finite element method

Abstract

The present invention discloses a method for analyzing the stability of soil slopes based on the precise finite element method, which includes establishing a two-dimensional slope model and dividing the grid; using the modified hyperbolic circularized Mohr-Coulomb yield constitutive model to determine whether the soil mass yields. If yielding occurs, an accurate explicit Euler algorithm will be used for stress redistribution, and the stress point will always be ensured to be on or within the yield surface through stress back-pulling; for the solution of the non-linear equations, Newton-Raphson iteration is used for solution; the strength reduction method is introduced, and a failure criterion with sufficient accuracy requirements is adopted. The safety factor of the slope is obtained through iterative search, and data such as the displacement and plastic zone of the slope can also be obtained through post-processing of the results. The present invention is reasonable and accurate for the stability analysis of slopes, can give a high-precision slope safety factor, and at the same time, based on the data such as the slope displacement and plastic zone calculated, it can provide a reasonable basis or suggestion for the implementation of slope engineering prevention and control measures.

Description

Technical Field

[0001] The present invention belongs to the field of computer simulation calculations, involves computer finite element modeling and computational analysis, and specifically relates to a method for analyzing the stability of soil slopes based on the precise finite element method. Background Art

[0002] Landslide disasters often cause damage to rivers, farmland, highways, railways, and other buildings along transportation lines, seriously threatening human life and property safety. How to calculate a more accurate slope safety factor and how to obtain data such as in-situ stress, strain, displacement, or deformation of slopes is crucial for the engineering prevention and control of slopes. Traditional limit equilibrium methods and limit analysis methods usually require prior knowledge or assumption of the slip surface and assume that the sliding mass is a rigid body. Many such assumptions limit the use of these methods and make it difficult to analyze the stability of complex slopes, such as heterogeneous slopes considering spatial variability, and at the same time, the deformation and stress-strain state of slopes cannot be obtained.

[0003] However, with the development of computer technology, the problem of large-scale data calculation has been overcome, and the finite element method has developed rapidly. The finite element method does not have many assumptions of the limit equilibrium method and limit analysis method and is applicable to the stability analysis of various complex slopes, which is well received by scholars.

[0004] In view of this, based on the finite element method and the strength reduction method, this paper simultaneously uses a high-precision explicit Euler algorithm for stress redistribution and performs a pull-back on the redistributed stress each time to ensure that the stress points always remain on the yield surface. In order to avoid the singularity of the traditional Mohr-Coulomb yield constitutive model, it is corrected by hyperbolic circularization. Finally, a more accurate slope safety factor and the slope deformation and stress-strain state are obtained. Summary of the Invention

[0005] Based on a precise finite element method, the present invention invents a method for analyzing the stability of soil slopes, which is applicable to various complex slopes, can also consider the spatial variability of soil, can accurately solve the safety factor of slopes, and at the same time, through post-processing of the program, the stress-strain, displacement, or deformation state at each position of the slope can be obtained, and the potential plastic zone or slip surface position of the slope can be provided. It provides a reliable basis for the engineering prevention and control of landslides.

[0006] The technical solution adopted in this paper is as follows:

[0007] A method for analyzing the stability of soil slopes based on the precise finite element method, comprising the following steps:

[0008] Step 1: Obtain the geometric parameters and boundary conditions of the slope, establish a two-dimensional finite element model of the slope and solve it to obtain the initial displacement, strain, and stress of the elements;

[0009] Step 2: Take the stress at the initial Gauss point as the predicted stress, and consider the modified hyperbolic arc-shaped Mohr-Coulomb yield constitutive model to determine whether the Gauss point yields. If it yields, stress redistribution is required. If not, directly judge whether the next Gauss point yields until all Gauss points of all elements are judged;

[0010] Step 3: Since stress redistribution occurs, the problem of solving the linear equations in the elastic state becomes the problem of solving the non-linear equations in the elastoplastic state. The Newton-Raphson iteration is used for solution. According to the redistributed stress and the plastic flow rule, the elastoplastic stress-strain matrix is obtained, and a new global stiffness equation is constructed;

[0011] Step 4: Calculate the plastic strain at the Gauss point according to the plastic flow rule. At the same time, with the iteration, the element nodal displacements, strains, stresses, and plastic strains are continuously accumulated until the set maximum number of iterations is reached, which is regarded as the slope has failed, or the convergence criterion is met;

[0012] Step 5: Introduce the strength reduction method for calculation. Continuously change the reduction factor by the bisection method and repeat Steps 1 to 4 until the ultimate failure state of the slope is searched and the accuracy requirement is met. The reduction factor at this time is the safety factor of the slope.

[0013] Further, according to the displacement calculated in Step 1, the strain increment Δε i and stress increment Δσ i are obtained, and the predicted stress at the Gauss point is set to judge whether yielding occurs:

[0014] Δε i = BΔu i , Δσ i = D e Δε i

[0015] σ i = σ i-1 + Δσ i

[0016] Δu i is the displacement increment at the i-th iteration, B is the element geometric matrix, D e is the elastic stress matrix of the element, Δε i is the strain increment at the i-th iteration, Δσ i is the stress increment at the i-th iteration,

[0017] Use the modified hyperbolic arc-shaped Mohr-Coulomb yield constitutive model to judge the yielding of the Gauss point:

[0018]

[0019]

[0020]

[0021]

[0022]

[0023] s x = σ x -σ m , s y = σ y -σ m , s z = σ z -σ m

[0024] s xy = σ xy , s xz = σ xz , s yz = σ yz

[0025]

[0026]

[0027]

[0028]

[0029] J2 = -(s x s y + s x s z + s y s z ) + s xy 2 + s xz 2 + s yz 2

[0030] J3 = s x s y s z + 2s xy s xz s yz + s x s yz 2 + s y s xz 2 + s z s xy2

[0031] Let \(F\) be the yield function, \(\sigma\) x , \(\sigma\) y , \(\sigma\) z , \(\sigma\) xy , \(\sigma\) xz , \(\sigma\) yz are respectively six stress components, \(c\), are respectively the cohesion and friction angle of the soil mass, \(a\), \(\theta\) T are two parameters for the hyperbolic circular arc approximation, usually \(\theta\) T \( = 25^{\circ}\) can meet the accuracy requirements; \(J_2\) and \(J_3\) are respectively the second and third invariants of the stress deviator tensor.

[0032] Furthermore, in step two, if yielding occurs for the first time, the proportionality factor for the elastic - plastic transition stage needs to be calculated to obtain the plastic stress increment for subsequent stress redistribution calculation and the stress before redistribution:

[0033] Among them, the calculation formula for the proportionality factor is:

[0034]

[0035]

[0036] Let \(z\) be the required proportionality factor, \(\Delta\sigma\) e be the elastic stress increment from the elastic to the plastic stage, \(F_0\), \(F_1\), \(F_2\) are respectively the yield function values at the stresses \(\sigma\) A , \(\sigma\) A +\(\Delta\sigma\) e , \(\sigma\) A + \(z_1\Delta\sigma\) e respectively, is the derivative value of the yield function at \(\sigma\) A + \(z_1\Delta\sigma\) e ;

[0037] If it enters the plastic state for the second time, the proportionality factor \(z\) is taken as zero, and the elastic strain increment is all used as the plastic strain increment for stress redistribution.

[0038] Furthermore, in step two, the explicit Euler algorithm proposed by Abbo is used for stress redistribution. For each Gauss point where yielding occurs, stress redistribution is carried out one by one, which is divided into the following five steps:

[0039] (1) The stress before redistribution is:

[0040] \(\sigma\) i =\(\sigma\) i-1 + \(z\Delta\sigma\) e

[0041] (2) Set the pseudo-time. The initial pseudo-time T = 0, and the pseudo-time increment step ΔT = 1. The stress redistribution ends until T = 1. The stress increment at each pseudo-time step is:

[0042] Δσ i = ΔTΔσ e -Δλ i D e b i

[0043]

[0044] Adopt the ideal elastoplastic model, A = 0, Δλ i The plastic multiplier increment at the i-th pseudo-time step, b i is the partial derivative of the plastic potential function with respect to stress at the i-th pseudo-time step, a i is the partial derivative of the yield function with respect to stress at the i-th pseudo-time step. Considering the associated flow rule, the plastic potential function is equal to the yield function;

[0045] To obtain an accurate stress increment, a higher-order explicit Euler method is adopted:

[0046]

[0047] The error at the current sub-step is:

[0048]

[0049] (3) If R T+ΔT > STOL, the sub-increment step fails, and a smaller pseudo-time increment step needs to be taken. Then, the size of the pseudo-time increment step is modified by a coefficient q:

[0050]

[0051] ΔT = max{qΔT, ΔT min}

[0052] where the allowable error STOL is set to 0.05, and the minimum relative error Eps is set to 10 -16 . If the sub-increment step fails, it is necessary to return to step (2) with a smaller pseudo-time increment to recalculate the stress increment. Otherwise, accept the updated stress in step (2). This pseudo-time increment step is valid and accepted:

[0053] T = T + ΔT

[0054] (4) Calculate the stress at the sub-step that meets the requirements from step (3). Further judge whether this stress point is on the yield surface. If it is outside the yield surface, stress backtracking is required. After the calculation is completed, calculate the next sub-step:

[0055]

[0056] If the stress back-pull fails, let q = min{q, 1}, and set the new pseudo-time increment step for the next iteration:

[0057] ΔT = min{ΔT, ΔT min}, ΔT = min{ΔT, 1 - T}

[0058] ΔT = qΔT

[0059] where the minimum pseudo-time increment step ΔT is set min = 1;

[0060] (5) Repeat steps (2)-(4) until T = 1 is satisfied. At this time, the stress is the redistributed stress.

[0061] Furthermore, it is necessary to judge whether the redistributed stress in step (4) yields. If the redistributed stress is still outside the yield surface, the stress needs to be corrected and back-pulled to the yield surface. The calculation steps are as follows:

[0062] (1) For the stress σ0 that has been redistributed but not corrected at a certain Gauss point, calculate the yield function value F0 at this time;

[0063] (2) Correct the initial uncorrected stress:

[0064] σ = σ0 + δσ = σ0 - δλD e b0

[0065]

[0066] a0 and b0 are the derivatives of the yield function and the plastic potential function under the stress σ0, respectively.

[0067] (3) Since the corrected stress may be further away from the yield surface than the initial stress during the stress correction process, it is necessary to make a further judgment on the corrected stress, that is: if |F(σ)| > |F(σ0)|, then discard the previous correction and perform the following correction:

[0068] σ = σ0 - δλa0

[0069]

[0070] (4) Repeat steps (2)-(3) 10 times. If the corrected stress satisfies |F(σ)| ≤ FTOL at this time, exit the loop;

[0071] If the number of cycles is less than 10, it is considered that the stress correction is successful, and the stress σ0 is updated to σ. Otherwise, it is considered that the stress backhaul fails, and a smaller pseudo-time increment step needs to be set for calculation.

[0072] Further, in step three, the stress-strain curve is no longer linear, and the elastoplastic stress-strain matrix D is obtained through the plastic flow rule. ep :

[0073]

[0074]

[0075] Reconstruct a new stiffness matrix, and then solve the large-scale nonlinear equations according to the Newton-Raphson iteration:

[0076] Δu n = -(K T (u n )) -1 ψ(u n ) = (K T n ) -1 (Q - P(u n ))

[0077] (K T ) n is the elastoplastic stiffness matrix obtained by the nth-step iterative calculation, Q is the external load, and ψ(u n ) is the unbalanced force at the nth-step iteration;

[0078] P is the internal force during the nth-step iteration process:

[0079] P = ∫∫∫B T σ n dΩ.

[0080] Further, in step four, the Newton-Raphson iteration method is adopted, the maximum number of iterations is set to 25, and the iteration convergence criteria for displacement and unbalanced force are set:

[0081] ||Δu i || ≤ Tol·||u i ||, ||ψ i || ≤ Tol·||Q||

[0082] where Tol is taken as 0.00001. During the iteration process, displacement, stress, strain, and plastic strain are all continuously accumulated. Among them, the plastic strain is calculated according to the flow rule, and the plastic strain calculation formula is:

[0083]

[0084]

[0085] Q is the plastic potential function.

[0086] Furthermore, in step five, in order to consider the slope safety reserve problem, the strength reduction method is introduced:

[0087]

[0088] Where FS is the reduction coefficient, which is used for binary search in the process of solving the slope safety factor. The calculations of steps one to four are continuously repeated until FS converges to the state where the slope is just at the limit state. This reduction coefficient FS is the slope safety factor.

[0089] Furthermore, according to the stress and plastic strain of each element Gauss point obtained in step four, the stress or strain at each position of the slope can be obtained. First, the stress or plastic strain of each element node is obtained by extrapolation:

[0090]

[0091]

[0092] σ i σ j σ m σ p are the nodal stresses of the element to be solved, and σ1, σ2, σ3, σ4 are the stresses of the 4 Gauss points of the unsmoothed element. At this time, only the nodal stresses of each element are obtained, and the stresses of each node of the slope still need to be obtained by the weighted average method:

[0093]

[0094] is the stress of each element containing the node, and s j is the area of each element containing node i;

[0095] In order to more intuitively reflect the position of the slope plastic zone, in the plane strain problem, the four plastic strain components obtained for each node can be calculated to obtain the equivalent plastic strain through the following formula;

[0096]

[0097] On the other hand, this application also protects an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the computer program, it implements the analysis method described in any one of claims 1 to 9.

[0098] Compared with the prior art, the present invention has the following advantages:

[0099] 1. This application adopts a high - order stress redistribution method: the explicit Euler algorithm, which has higher calculation accuracy than the visco - elastic - plastic method used in traditional finite element methods, and can obtain more accurate slope safety factors and data such as slope displacements, stresses, and strains.

[0100] 2. Compared with other methods, such as slope stability analysis methods like the limit equilibrium method and the limit analysis method, this method does not require prior assumption of the slip surface, can calculate the safety factors of various complex heterogeneous slopes, and can also consider the influence of soil spatial variability and other properties.

[0101] 3. Some numerical software calculations only consider the Mohr - Coulomb yield constitutive model, or it is difficult to directly modify the constitutive model according to requirements because the embedded programs in the software are not visible. This method can modify the yield constitutive model for further solution according to different requirements. For example, this method adopts a modified hyperbolic circularized Mohr - Coulomb yield constitutive model. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the following - described drawings are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0103] Figure 1 It is a calculation flow chart of a soil slope stability analysis method based on the precise finite element method in this application.

[0104] Figure 2 It is a reference diagram of the slope calculation model of a soil slope stability analysis method based on the precise finite element method in this application.

[0105] Figure 3 It is a four - node rectangular element used in the calculation of a soil slope stability analysis method based on the precise finite element method in this application.

[0106] Figure 4 It is the modified hyperbolic circularized Mohr - Coulomb yield surface on the π - plane of a soil slope stability analysis method based on the precise finite element method in this application, which can eliminate the singularity at the tip and edge of the traditional Mohr - Coulomb yield surface.

[0107] Figure 5 It is the transition of the Gaussian point stress from the elastic stage to the plastic stage for a soil slope stability analysis method based on the precise finite element method in this application, used to solve the scale factor.

[0108] Figure 6This application is a method for analyzing the stability of soil slopes based on the precise finite element method, showing the calculation flow chart of stress redistribution under the Euler algorithm.

[0109] Figure 7 This application is a method for analyzing the stability of soil slopes based on the precise finite element method. During the stress redistribution process, it shows the calculation flow chart of stress pulling back to the yield surface.

[0110] Figure 8 This application is a method for analyzing the stability of soil slopes based on the precise finite element method, showing the calculation flow chart of the Newton-Raphson iteration method for calculating non-linear equations.

[0111] Figure 9 This is the slope deformation diagram (magnified) obtained by post-processing the computer program of a method for analyzing the stability of soil slopes based on the precise finite element method in this application.

[0112] Figure 10 This is the cloud diagram of the slope plastic zone obtained by post-processing the computer program of a method for analyzing the stability of soil slopes based on the precise finite element method in this application. Detailed implementation manners

[0113] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0114] The following combines Figure 1-10 , and specifically introduces the embodiments of this application.

[0115] As Figure 1 shown is the calculation flow chart of the present invention, which summarizes and introduces the calculation process. The specific steps are as follows:

[0116] Step 1: Obtain the following parameters of the slope: slope geometric parameters, unit weight parameters, strength parameters, and deformation parameters.

[0117] The slope geometric parameters include the top slope length w1, the horizontal slope length s1, the bottom slope length w2, the slope height h1, the base thickness h2, the number of divided units nx1 at the top of the slope, the number of divided units nx2 at the bottom of the slope, the number of divided units ny1 for the slope height, and the number of divided units ny2 for the base thickness; the specific parameter values are as Figure 2 shown.

[0118] The soil parameters include the unit weight γ of the soil, and the strength parameters include the cohesion c and the friction angle The deformation parameters include the elastic modulus E and the Poisson's ratio μ. If

[0119] Step 2: Establish a two-dimensional slope model, perform mesh division on the model, obtain the element numbers and node coordinate information of each element, assign material parameters to the elements, set the boundary conditions of the model, and construct a complete finite element slope calculation model.

[0120] The establishment of the slope model and mesh division can be generated through self-written programs, MATLAB, or ANSYS software. Obtain data such as the element number sequence, node coordinates, and boundary conditions, and assign material parameters to each element for subsequent finite element calculations. If the spatial variability of the soil mass needs to be considered, it is necessary to assign the soil parameters of each element separately.

[0121] Step 3: According to the element and node information, construct the element displacement field in the local coordinate system ξ-η by interpolation, introduce the shape function, use Gaussian integration to construct the stiffness matrix and equivalent nodal force matrix of the element and integrate them, and consider the boundary conditions of the model to form the initial overall stiffness equation.

[0122] Construct the element displacement field:

[0123]

[0124] For a four-node quadrilateral element, as Figure 3 shown, the displacement field is:

[0125]

[0126] N is the shape function matrix and q is the nodal displacement matrix.

[0127] The shape function expression of a certain element n node in the local coordinate system is:

[0128] N n =(1 + ξ n ξ)(1 + η n η) / 4

[0129] Construct the strain field as:

[0130] ε(x, y)=B·q

[0131] B is the geometric matrix, written in block form:

[0132] B = [B i B j B m B p

[0133] Where:

[0134] ​

[0135] The constructed stress field is as follows:

[0136]

[0137] Considering the plane strain problem, when in the elastic state, the elastic matrix D is:

[0138]

[0139] Construct the element stiffness equation:

[0140] Kq = f

[0141] Among them, both the element stiffness matrix K and the equivalent nodal force matrix f are solved using Gaussian quadrature:

[0142]

[0143] G(ξ, η) = B T DB|J|

[0144]

[0145] Among them, f v is the equivalent nodal force of body force, H is the weight, and J is the Jacobian matrix for the transformation between local coordinates and global coordinates:

[0146]

[0147] Integrate the element stiffness equation into the global stiffness equation according to the node numbers. During the process of solving this linear system of equations, to avoid the singularity of the equations, the large number multiplication method is adopted to handle the boundary conditions and modify the existing stiffness equation. An example of calculation using the large number multiplication method is as follows: Now assume the boundary conditions are:

[0148] u1 = β1 u2 = β3

[0149] β1 and β3 are the known displacement values of the nodes.

[0150] The existing stiffness equation is:

[0151]

[0152] Modify it using the large number multiplication method to:

[0153]

[0154] Step 4: Solve the initial global stiffness equation (large linear system of equations) to obtain the initial displacements, strains, and stresses of the elements.

[0155] The linear equations are solved by using the Gauss elimination method or the Gauss-Seidel iteration method, or the built-in solution method in MATLAB can also be used.

[0156] Step Five: Take the stress at the initial Gauss point as the predicted stress, and consider the modified hyperbolic arc-shaped Mohr-Coulomb yield constitutive model to judge whether the Gauss point yields. If yielding occurs, stress redistribution is required. If no yielding occurs, directly judge whether the next Gauss point yields until all Gauss points of all elements are calculated.

[0157] According to the displacement calculated in the fourth step, the strain increment Δε i and the stress increment Δσ i are further obtained, and the predicted stress at the Gauss point is set to judge whether yielding occurs:

[0158] Δε i = BΔu i Δσ i = D e Δε i

[0159] σ i = σ i-1 +Δσ i

[0160] Adopt the modified hyperbolic arc-shaped Mohr-Coulomb yield constitutive model, as shown in Figure 4, and judge whether each Gauss point yields in turn:

[0161] The expression of the yield function is:

[0162]

[0163] Where:

[0164]

[0165]

[0166]

[0167]

[0168] s x = σ x -σ m ,s y = σ y -σ m ,s z = σ z -σ m

[0169] s xy = σ xy,s xz = σ xz ,s yz = σ yz

[0170]

[0171]

[0172]

[0173]

[0174] Step 6: If yielding occurs for the first time, the proportionality factor for the elastic-to-plastic transition stage needs to be calculated, as shown in Figure 5 .

[0175] The calculation formula for the proportionality factor is

[0176]

[0177]

[0178] If it enters the plastic state for the second time, the proportionality factor z is taken as zero, and the elastic strain increment is all regarded as the plastic strain increment for stress redistribution.

[0179] Step 7: Use the explicit Euler algorithm proposed by Abbo for stress redistribution. The calculation process is as shown in Figure 6 .

[0180] Determine the plastic strain increment according to the proportionality factor obtained in Step 6, and then use the explicit Euler algorithm proposed by Abbo for stress redistribution. Stress redistribution needs to be carried out for each Gauss point where yielding occurs. Now, analyze the stress of a certain Gauss point where yielding occurs, which can be mainly divided into the following five sub-steps:

[0181] (1) The stress before redistribution is:

[0182] σ i = σ i-1 + zΔσ e

[0183] (2) Set the pseudo-time. Initially, T = 0, ΔT = 1, until T = 1 and the stress redistribution ends. The stress increment at each time step is:

[0184] Δσ i = ΔTΔσ e - Δλ i D e b i

[0185]

[0186] Using the ideal elastoplastic model, A = 0, b is the partial derivative of the plastic potential function with respect to stress, a is the partial derivative of the yield function with respect to stress, and considering the associated flow rule, the plastic potential function is equal to the yield function.

[0187] In order to obtain accurate stress increments, a higher-order calculation method is adopted.

[0188]

[0189] The error of the current sub-step is:

[0190]

[0191] (3) If R T+ΔT > STOL, the sub-increment step fails and a smaller time increment step needs to be taken. Then, the size of the time increment step is modified according to the following formula:

[0192]

[0193] ΔT = max{qΔT, ΔT min}

[0194] where the minimum relative error Eps is set to 10 -16 , the allowable error STOL is set to 0.05. If the sub-increment step fails, it is necessary to return to step (2) of step seven with a smaller pseudo-time increment to recalculate the stress increment. Otherwise, accept the updated stress in step (2). This pseudo-time increment step is valid and accepted:

[0195] T = T + ΔT

[0196] (4) Calculate the stress at the sub-step that meets the requirements obtained from step (3) of step seven. Further, judge whether this stress point is on the yield surface. If it is outside the yield surface, perform stress backtracking according to step eight. After the calculation is completed, perform the calculation of the next sub-step:

[0197]

[0198] If the stress backtracking fails, let q = min{q, 1}, and set the new pseudo-time increment step for the next iteration:

[0199] ΔT = min{ΔT, ΔT min},ΔT = min{ΔT, 1 - T}

[0200] ΔT = qΔT

[0201] where the minimum pseudo-time increment step ΔT min = 1.

[0202] (5) Repeat steps (2) - (4) of step seven until T = 1 is satisfied. At this time, the stress is the required redistributed stress.

[0203] Step eight: During the stress redistribution process, if the Gaussian point stress is still outside the yield surface, the stress needs to be corrected and pulled back to the yield surface. The calculation flow chart is as Figure 7 shown.

[0204] According to the redistributed stress obtained in step (4) of step seven, it is also necessary to judge whether it yields. If the redistributed stress is still outside the yield surface, the stress needs to be corrected and pulled back to the yield surface again. The calculation is divided into the following 5 small steps:

[0205] (1) Assume the initial uncorrected stress σ0 for a certain Gaussian point and calculate the yield function value F0 at this time.

[0206] (2) Correct the initial uncorrected stress:

[0207] σ = σ0 + δσ = σ0 - δλD e b0

[0208]

[0209] (3) Since during the stress correction process, it is possible that the corrected stress is further away from the yield surface than the initial stress. At this time, it is necessary to make a further judgment on the corrected stress. If |F(σ)| > |F(σ0)|, discard the previous correction and perform the following correction:

[0210] σ = σ0 - δλa0

[0211]

[0212] (4) Repeat steps (2) - (3) of step eight 10 times. If the corrected stress satisfies |F(σ)| ≤ |FTOL| at this time, exit the loop.

[0213] (5) If the number of loops is less than 10 times, it is considered that the stress correction is successful, and update the stress σ0 = σ. Otherwise, it is considered that the stress pull-back fails, and a smaller pseudo-time increment step will be set according to step (4) of step seven for calculation.

[0214] Step nine: Due to the occurrence of stress redistribution, the problem of solving the linear equation system in the elastic state becomes the problem of solving the non-linear equation system in the elastoplastic state. The Newton-Raphson iteration is used for solving. According to the redistributed stress and the plastic flow rule, the elastoplastic stress-strain matrix is obtained, and a new overall stiffness equation is constructed.

[0215] When the Gauss point yields, the stress-strain curve shows a non-linear relationship, and the elastoplastic stress-strain matrix is obtained through the plastic flow rule:

[0216]

[0217]

[0218] Reconstruct a new stiffness matrix according to the algorithm in step three, and then solve the large non-linear equations set by Newton-Raphson iteration. The specific iteration flowchart is as Figure 8 shown, and the iteration calculation formula is:

[0219] Δu n =-(K T (u n )) -1 ψ(u n )=(K T n ) -1 (F-P(u n ))

[0220] ψ(u n ) is the unbalanced force, F is the external load, and P is the internal force in a certain iteration process:

[0221] P=∫∫∫B T σ n dΩ

[0222] Step ten: Calculate the plastic strain of the Gauss point according to the plastic flow rule. At the same time, with the iteration, the element node displacements, strains, stresses, and plastic strains are continuously accumulated until the set maximum number of iterations (regarded as the slope has failed) or the convergence criterion is met.

[0223] The maximum number of iterations set by the Newton-Raphson iteration method is 25 times, and the iteration convergence criteria for displacement and unbalanced force are set:

[0224] ||Δu i ||≤Tol·||u i ||, ||ψ i ||≤Tol·||Q||

[0225] where Tol is taken as 0.00001.

[0226] During the iteration process, the displacements, stresses, strains, and plastic strains will be continuously accumulated. Among them, the plastic strain is calculated according to the flow rule,

[0227] and the calculation formula is:

[0228]

[0229]

[0230] Step 11: Introduce the strength reduction method to perform the above calculations, and obtain the safety factor of the slope by means of bisection search.

[0231] Considering the problem of slope safety reserve, introduce the strength reduction method to obtain the slope safety factor. The reduction calculation formula is

[0232]

[0233] Continuously change the reduction coefficient by bisection method, set the upper and lower limits of the reduction coefficient. When the slope fails at a certain reduction coefficient, change the upper limit value of the reduction coefficient, and vice versa, change the lower limit value of the reduction coefficient, and continuously bisect until the difference between the upper and lower limits of the reduction coefficient is less than the error limit, usually set to 0.01.

[0234] Specifically: In the process of solving the slope safety factor, set two reduction coefficients, one large and one small, in advance for bisection search (FSmin = 0, FSmax = 5, which can be adjusted according to different slopes). The initial reduction coefficient FS = (FSmin + FSmax) / 2. Then reduce the strength parameters c and φ according to this FS, and calculate Steps 1 to 4. If the calculation does not converge, it is considered that the reduction coefficient FS is too large and the slope has failed, and let FSmax = FS. On the contrary, if the calculation converges, it is considered that the reduction coefficient FS is too small and the slope has not failed, and let FSmin = FS. Then calculate a new reduction coefficient according to FS = (FSmin + FSmax) / 2, and repeat Steps 1 to 4 until the slope is just in the limit state under a certain reduction coefficient. This reduction coefficient is the slope safety factor. The specific implementation of this method is when (FSmax - FSmin) / 2 ≤ err (err = 0.01, the calculation accuracy can be adjusted according to requirements), it is considered that FS = (FSmin + FSmax) / 2 at this time is the slope safety factor.

[0235] Step 12: Perform post-processing of the program, and finally obtain data such as the displacements, strains, stresses, and plastic strains of each node in the model. The deformation diagram and plastic zone nephogram of the slope can be obtained by plotting with Tecplot. The results are shown as Figure 9-10 shown.

[0236] After the calculation of Step 12 is completed, the stresses and plastic strains of each Gauss point of each element can be obtained. In order to more intuitively reflect the stress or strain state at each position of the slope, it is necessary to obtain the state of each node. First, obtain the stress or plastic strain of each element node by extrapolation:

[0237]

[0238]

[0239] σ i ,σ j ,σ m ,σ p is the nodal stress of the unit to be obtained, and σ1, σ2, σ3, σ4 are the stresses at the 4 Gauss points of the unsmoothed unit.

[0240] However, only the nodal stresses of each unit are obtained at this time, and the stresses of each node of the slope still need to be obtained by the method of weighted average:

[0241]

[0242] The method for solving the plastic strain of the unit is similar. Usually, in order to more intuitively reflect the position of the plastic zone of the slope, in plane strain problems, the four plastic strain components obtained at each node can be calculated through the following formula to obtain the equivalent plastic strain.

[0243]

[0244] The method illustrated in the above embodiments can be specifically implemented by a computer chip or an entity, or by a product with certain functions. A typical implementation device is a computer. Specifically, the computer can be, for example, a personal computer, a laptop computer, a cellular phone, a camera phone, a smart phone, a personal digital assistant, a media player, a navigation device, an email device, a game console, a tablet computer, a wearable device, or a combination of any of these devices.

[0245] For the convenience of description, the above devices are described by function as various units respectively. Of course, when implementing this specification, the functions of each unit can be implemented in the same or multiple software and / or hardware.

[0246] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0247] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory produce a manufactured article including an instruction device, and the instruction device implements the process Figure 1 one process or multiple processes and / or blocks Figure 1The functions specified in one or more boxes.

[0248] In a typical configuration, a computing device includes one or more processors (CPUs), an input / output interface, a network interface, and memory.

[0249] The memory may include non-permanent memory in the form of computer-readable media, random access memory (RAM), and / or non-volatile memory such as read-only memory (ROM) or flash memory (flash RAM). Memory is an example of computer-readable media.

[0250] Accordingly, the present application provides a computer-readable storage medium storing a computer program, which when executed by a processor implements the method described in one of the foregoing.

[0251] The present application further provides an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein the processor implements the analysis method described in one of the foregoing when executing the computer program.

[0252] Computer-readable media includes both permanent and non-permanent, removable and non-removable media and can be implemented by any method or technology for storing information. The information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, compact disc read-only memory (CD-ROM), digital versatile discs (DVD) or other optical storage, magnetic cassettes, magnetic tape magnetic disk storage or other magnetic storage devices, or any other non-transmission media that can be used to store information accessible by a computing device. As defined herein, computer-readable media does not include transitory media such as modulated data signals and carrier waves.

[0253] Those skilled in the art should understand that the embodiments of this specification can be provided as a method, system, or computer program product. Therefore, this specification can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Moreover, this specification can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0254] This specification may be described in the general context of computer-executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. The specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices connected through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including storage devices.

[0255] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some of the technical features. However, such modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for analyzing the stability of soil slopes based on the precise finite element method, characterized in that: It includes the following steps: Step 1: Obtain the slope geometric parameters and boundary conditions, establish a two-dimensional finite element model of the slope and solve it to obtain the initial displacements, strains, and stresses of the elements; Step 2: Take the stress at the initial Gauss point as the predicted stress, and consider the modified hyperbolic circularized Mohr-Coulomb yield constitutive model to judge whether this Gauss point yields. If it yields, stress redistribution is required. If it does not yield, directly judge whether the next Gauss point yields until all Gauss points of all elements are judged; Step 3: Due to stress redistribution, the problem of solving the linear equation system in the elastic state becomes the problem of solving the non-linear equation system in the elastoplastic state. The Newton-Raphson iteration is used for solution. According to the redistributed stress and the plastic flow rule, the elastoplastic stress-strain matrix is obtained, and a new overall stiffness equation is constructed; Step 4: Calculate the plastic strain of the Gauss point according to the plastic flow rule. At the same time, with the iteration, the element node displacements, strains, stresses, and plastic strains are continuously accumulated until the set maximum number of iterations is reached, which is regarded as the slope has failed, or the convergence criterion is met; Step 5: Introduce the strength reduction method for calculation. Continuously change the reduction factor by the bisection method and repeat Steps 1 to 4 until the ultimate failure state of the slope is searched and the accuracy requirement is met. At this time, the reduction factor is the safety factor of the slope.

2. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 1, characterized in that: Based on the displacement calculated in Step 1, the strain increment Δε is further obtained. i and the stress increment Δσ i , and the Gaussian point predicted stress is set to judge whether yielding occurs: Not i = BΔu i , Δσ i = D e Δε i σ i = σ i-1 + Δσ i Δu i is the displacement increment at the i-th iteration step, B is the element geometric matrix, D e is the elastic stress matrix of the element, Δε i is the strain increment at the i-th iteration step, Δσ i is the stress increment at the i-th iteration step, The modified hyperbolic circularized Mohr-Coulomb yield constitutive model is used to judge the yield of the Gauss point: s x = σ x - σ m ,s y = σ y - σ m ,s z = σ z - σ m s xy = σ xy ,s xz = σ xz ,s yz = σ yz J2 = -(s x s y + s x s z + s y s z ) + s xy 2 + s xz 2 + s yz 2 J3 = s x s y s z + 2s xy s xz s yz + s x s yz 2 + s y s xz 2 + s z s xy 2 F is the yield function, and σ x , σ y , σ z , σ xy , σ xz , σ yz are the six stress components respectively, and c, are the cohesion and friction angle of the soil mass respectively, and a, θ T are two parameters for the hyperbolic circular arc approximation. Usually θ T = 25° can meet the accuracy requirements; J2 and J3 are the second and third invariants of the stress deviator tensor respectively.

3. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 2, characterized in that: In Step 2, if it is the first time to yield, it is necessary to calculate the proportionality factor in the elastic-to-plastic transition stage to obtain the plastic stress increment for subsequent stress redistribution calculation and the stress before redistribution: The calculation formula for the proportionality factor is: F2 = F(σ A + z1Δσ e ) z is the required scale factor, Δσ e is the elastic stress increment from the elastic to the plastic stage, and F0, F1, and F2 are the values of the yield function under the stresses σ A , σ A +Δσ e , σ A +z1Δσ e respectively, and is the derivative value of the yield function under σ A +z1Δσ e ; If it enters the plastic state for the second time, the proportionality factor z is taken as zero, and the elastic strain increment is all used as the plastic strain increment for stress redistribution.

4. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 1, characterized in that: In Step 2, the explicit Euler algorithm proposed by Abbo is used for stress redistribution. Stress redistribution is carried out one by one for each yielding Gauss point, which is divided into the following five steps: (1) The stress before redistribution is: σ i = σ i-1 + zΔσ e (2) Set the pseudo-time. The initial pseudo-time T = 0, and the pseudo-time increment step ΔT = 1 until T = 1 and the stress redistribution ends. The stress increment at each pseudo-time step is: Δσ i = ΔTΔσ e -Δλ i D e b i Using the ideal elastoplastic model, A = 0, Δλ i Increment of the plastic multiplier at the i-th pseudo-time step, b i Partial derivative of the plastic potential function with respect to stress at the i-th pseudo-time step, a i Partial derivative of the yield function with respect to stress at the i-th pseudo-time step. Considering the associated flow rule, the plastic potential function is equal to the yield function; In order to obtain an accurate stress increment, a higher-order explicit Euler method is used: The error at the current sub-step is: (3) If R T+ΔT > STOL, the sub-increment step fails and a smaller pseudo-time increment step needs to be taken. Then, the size of the pseudo-time increment step is modified by a coefficient q: ΔT = max{qΔT, ΔT min} where the allowable error STOL is set to 0.05 and the minimum relative error Eps is set to 10 -16 , if the sub-increment step fails, it is necessary to return to step (2) with a smaller pseudo-time increment to recalculate the stress increment. Otherwise, accept the updated stress in step (2). This pseudo-time increment step is valid and accepted: T = T + ΔT (4) Calculate the stress at the sub-step that meets the requirements from Step (3), and further judge whether this stress point is on the yield surface. If it is outside the yield surface, stress backtracking is required. After the calculation is completed, calculate the next sub-step: If the stress backtracking fails, let q = min{q, 1}, and set the new pseudo-time increment step for the next iteration: ΔT = min{ΔT, ΔT min}, ΔT = min{ΔT, 1 - T} ΔT = qΔT wherein a minimum pseudo-time increment step ΔT is set min = 1; (5) Repeat steps (2)-(4) until T = 1 is satisfied. At this time, the stress is the redistributed stress.

5. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 4, characterized in that: In step (4), it is also necessary to judge whether the redistributed stress yields. If the redistributed stress is still outside the yield surface, the stress needs to be corrected and pulled back to the yield surface. The calculation steps are as follows: (1) For a certain Gauss point, calculate the yield function value F0 of the redistributed but uncorrected stress σ0. (2) Correct the initial uncorrected stress: σ = σ0 + δσ = σ0 - δλD e b0 a0 and b0 are respectively the derivatives of the yield function and the plastic potential function under the stress σ0. (3) Since the corrected stress may be further away from the yield surface than the initial stress during the stress correction process, it is necessary to further judge the corrected stress, that is: if |F(σ)| > |F(σ0)|, then discard the previous correction and perform the following correction: σ = σ0 - δλa0 (4) Repeat steps (2)-(3) 10 times. If the corrected stress at this time satisfies |F(σ)| ≤ |FTOL, then exit the loop; (5) If the number of loops is less than 10 times, it is considered that the stress correction is successful, and the stress σ0 = σ is updated. Otherwise, it is considered that the stress pull-back fails, and a smaller pseudo-time increment step needs to be set for calculation.

6. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 1, characterized in that: In step three, the stress-strain curve no longer shows a linear relationship, and the elastoplastic stress-strain matrix D is obtained through the plastic flow rule ep : Reconstruct a new stiffness matrix, and then solve the large-scale nonlinear equations according to the Newton-Raphson iteration: Δu n = -(K T (u n )) -1 ψ(u n ) = (K T n ) -1 (Q - P(u n )) (K T ) n is the elastoplastic stiffness matrix obtained from the n-th step iterative calculation, Q is the external load, and ψ(u n ) is the unbalanced force at the n-th step iteration; P is the internal force during the nth step iteration. P = ∫∫∫ B T σ n dΩ。 7. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 1, characterized in that: In step four, the Newton-Raphson iteration method is adopted, and the maximum number of iterations is set to 25 times. At the same time, the iteration convergence criteria for displacement and unbalanced force are set: ||Δu i ||≤Tol·||u i ||,||ψ i ||≤Tol·||Q|| Where Tol is taken as 0.00001. During the iteration process, displacement, stress, strain, and plastic strain are continuously accumulated. Among them, the plastic strain is calculated according to the flow rule, and the plastic strain calculation formula is: Q is the plastic potential function.

8. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 7, characterized in that: In step five, in order to consider the slope safety reserve problem, the strength reduction method is introduced: Where FS is the reduction factor. During the solution process of the slope safety factor, it is used for dichotomy search, and the calculations of steps one to four are continuously repeated until FS converges to the slope just being in the limit state. This reduction factor FS is the slope safety factor.

9. A method for analyzing the stability of soil slopes based on the precise finite element method according to claim 1, characterized in that: According to the stress and plastic strain of each unit Gauss point obtained in step four, the stress or strain at each position of the slope can be obtained. First, the stress or plastic strain of each unit node is obtained by extrapolation: σ i ,σ j ,σ m ,σ p are the nodal stresses of the required element, and σ1, σ2, σ3, σ4 are the stresses at the 4 Gauss points of the unsmoothed element. At this time, only the stresses at the nodes of each element are obtained, and the stresses at each node of the slope still need to be obtained by the method of weighted average: is for containing the stress of each element of the node, s j is the area of each element containing node i; In order to more intuitively reflect the position of the slope plastic zone, in the plane strain problem, the four plastic strain components obtained for each node can be calculated to obtain the equivalent plastic strain through the following formula; 10. An electronic device, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, When the processor executes the computer program, the analysis method according to any one of claims 1 to 9 is implemented.

Images