Grid-based calculation method and system for yield acceleration of multi-layer soil slopes

Abstract

The present invention discloses a method and system for calculating the yield acceleration of multi-layer soil slopes based on grids, belonging to the technical field of geotechnical structure stability simulation. Existing geotechnical structure stability analysis solutions do not consider multi-layer soil structures, are not applicable to multi-layer soil slopes, and cannot accurately calculate the seismic yield acceleration. The method for calculating the yield acceleration of multi-layer soil slopes based on grids in the present invention obtains the most dangerous slope sliding surface and the corresponding yield acceleration data of a certain multi-layer soil slope by constructing a slope digital simulation unit, a hydro-mechanical coupling unit, a circular slip surface calculation unit, the most dangerous slope sliding surface, and an acceleration calculation unit. Therefore, the geological characteristics of the multi-layer soil slope structure can be fully considered to ensure that the obtained most dangerous sliding surface is consistent with the actual situation; at the same time, the hydro-mechanical coupling effect is considered, so that the seismic slope stability of water retaining structures can be accurately analyzed. Therefore, the present invention can achieve accurate simulation calculation of the yield acceleration.

Description

Technical Field

[0001] The present invention relates to a method and system for calculating the yield acceleration of multi-layer soil slopes based on grids, belonging to the technical field of geotechnical structure stability simulation. Background Art

[0002] Geotechnical structures such as earth-rock dams, levees, and natural slopes are generally composed of granular structures such as soil materials or sand and gravel materials, and are prone to sliding failures. External factors that cause sliding instability of soil slopes include increased pore water pressure in the soil caused by seepage, top cracking, top loading, sudden drawdown of the water level, and seismic action, etc. However, due to the periodicity, short-term nature, and randomness of seismic action, it is difficult to quantitatively evaluate the risk of slope sliding instability that may occur during an earthquake through numerical methods.

[0003] According to the requirements of the "Seismic Design Code for Hydraulic Structures in Hydropower Projects" NB / T35047-2015, it is necessary to calculate the seismic stability along the potential slip surface and the permanent deformation of the dam body according to the seismic action effect. In this analysis process, before calculating the permanent deformation of the slope sliding body, it is necessary to accurately calculate the potential slip surface (the most dangerous slope sliding surface) and the yield acceleration of the slope under seismic conditions.

[0004] However, in the actual calculation process, since soil slopes are often composed of multiple soil layers and the physical and mechanical properties between the soil layers are discontinuous, the failure conditions of different soil layers are different, forming a complex sliding surface composed of multiple segments of linear and circular arcs, making it difficult to accurately obtain the most dangerous sliding surface, and thus it is impossible to accurately calculate the seismic yield acceleration.

[0005] Furthermore, Chinese Patent Application (Publication No.: CN116663274A) discloses a method for evaluating the seismic stability of a slope reinforced with prestressed anchor rods, including the following steps: Step 1, establish a dynamic failure mechanism of a slope reinforced with prestressed anchor rods under seismic action; Step 2, establish an instantaneous displacement mode of a slope reinforced with prestressed anchor rods under seismic action, and obtain the time-recursive formulas of the anchor rod axial force, anchoring angle, slope height, and slope angle; Step 3, based on the upper bound method of limit analysis, deduce the variation of the yield acceleration with time; Step 4, decouple the coupling relationship between the yield acceleration and the seismic displacement based on the time history analysis method, so as to obtain the variation of the yield acceleration of the slope reinforced with prestressed anchor rods and the seismic permanent displacement; Step 5, use the obtained seismic permanent displacement to evaluate the seismic stability of the slope reinforced with prestressed anchor rods.

[0006] The disadvantages of the above solution are as follows: assuming the shape of the sliding surface as a logarithmic spiral and not considering the multi-layer soil structure, it is not applicable to the scenario of multi-layer soil slopes. The most dangerous sliding surface obtained may deviate greatly from the actual situation, so the seismic yield acceleration still cannot be accurately calculated. In addition, the above solution does not consider the hydro-mechanical coupling effect, and there are certain limitations in the seismic slope stability analysis of water retaining structures.

[0007] The information disclosed in this background technology is only used to understand the background of the inventive concept, so it may include information that does not constitute the prior art. Summary of the Invention

[0008] Aiming at the above problems or one of the above problems, the first object of the present invention is to provide a method and system for calculating the yield acceleration of multi-layer soil slopes based on grids. By constructing a slope digital simulation unit, a hydro-mechanical coupling unit, a circular arc sliding surface calculation unit, the most dangerous slope sliding surface, and an acceleration calculation unit, the most dangerous slope sliding surface and the corresponding yield acceleration data of a certain multi-layer soil slope can be obtained. Therefore, the geological characteristics of the multi-layer soil slope structure can be fully considered to ensure that the obtained most dangerous sliding surface is consistent with the actual situation; at the same time, the hydro-mechanical coupling effect is considered, so that the seismic slope stability of water retaining structures can be accurately analyzed, and thus the accurate simulation calculation of the yield acceleration can be realized.

[0009] Aiming at the above problems or one of the above problems, the second object of the present invention is to provide a method and system for calculating the yield acceleration of multi-layer soil slopes based on grids, which can realize the calculation of the yield acceleration of multi-layer soil structures under complex working conditions, accurately obtain the most dangerous slope sliding surface and the corresponding yield acceleration data, have a wide range of applications, high efficiency, be applicable to quickly evaluating the seismic slope stability of multi-layer soil slopes under hydro-mechanical coupling conditions, and be convenient for the integrated simulation design of multi-layer soil slopes and the rapid iteration of design schemes.

[0010] To achieve one of the above objects, the first technical solution of the present invention is:

[0011] A method for calculating the yield acceleration of multi-layer soil slopes based on grids, comprising the following steps:

[0012] Step 1, through a pre-constructed slope digital simulation unit, according to the geological data of a certain multi-layer soil slope, and coupling the hydraulic parameters and mechanical shear resistance information, a geometric grid model is generated;

[0013] Step 2, using a pre-constructed hydro-mechanical coupling unit, adopting the unsaturated and unsteady seepage algorithm, processing the geometric grid model to obtain a steady seepage field, and extracting the pore pressure results;

[0014] Step 3, using a pre-constructed circular arc sliding surface calculation unit, according to the pore pressure results, determining the most dangerous circular arc sliding surface;

[0015] Step 4: Use the pre - constructed slope sliding surface calculation unit. Based on the most dangerous circular arc sliding surface, solve for the yield acceleration coefficients of several slope sliding surfaces, and determine the most dangerous slope sliding surface according to the yield acceleration coefficients.

[0016] Step 5: Based on the pre - constructed acceleration calculation unit, use the local grid encryption algorithm to process the most dangerous slope sliding surface to obtain the yield acceleration data of a multi - soil - layer slope.

[0017] The present invention constructs a slope digital simulation unit, a hydro - mechanical coupling unit, a circular arc sliding surface calculation unit, the most dangerous slope sliding surface, and an acceleration calculation unit, obtaining the most dangerous slope sliding surface of a multi - soil - layer slope and the corresponding yield acceleration data. Therefore, the geological characteristics of the multi - soil - layer slope structure can be fully considered to ensure that the obtained most dangerous sliding surface is consistent with the actual situation. At the same time, considering the hydro - mechanical coupling effect, the seismic slope stability of the water - retaining structure can be accurately analyzed. Thus, accurate simulation calculation of the yield acceleration can be realized, and the solution is scientific, reasonable, and feasible.

[0018] Furthermore, the present invention can realize the calculation of the yield acceleration of the multi - soil - layer structure under complex working conditions, accurately obtain the most dangerous slope sliding surface and the corresponding yield acceleration data, with a wide application range and high efficiency. It is applicable to quickly evaluating the seismic slope stability of the multi - soil - layer slope structure under hydro - mechanical coupling conditions, facilitating the integrated simulation design of the multi - soil - layer slope and the rapid iteration of the design scheme.

[0019] As a preferred technical measure:

[0020] Step 1: Through the pre - constructed slope digital simulation unit, according to the geological data of a multi - soil - layer slope, and coupling the hydraulic parameters and mechanical shear resistance information, the method for generating a geometric grid model is as follows:

[0021] Obtain the geological data of a multi - soil - layer slope, which includes the structural material zoning characteristics, the foundation layering situation, and the slope rate change of the slope surface.

[0022] According to the structural material zoning characteristics, the foundation layering situation, and the slope rate change of the slope surface, determine the slope sliding position information.

[0023] The slope sliding position information at least includes the position information of each level of platform on the slope surface, the position information of slope rate mutation, and the position information of the intersection of the layer and the slope surface.

[0024] Based on the slope sliding position information, adopt the finite - element strength reduction method to select the search area for the sliding - in point and the sliding - out point of the sliding surface.

[0025] According to the search area, establish line elements and surface elements.

[0026] Encrypt the line elements and surface elements to obtain the first-order discrete grid elements for the multi-layer soil slope.

[0027] Obtain the hydraulic parameters and mechanical shear resistance information, including the permeability coefficient, water retention curve, density, internal friction angle, and cohesion; among them, the water retention curve is used to describe the relationship between the effective saturation and pore pressure.

[0028] Couple the hydraulic parameters, mechanical shear resistance information with the first-order discrete grid elements to form a geometric grid model.

[0029] As a preferred technical measure:

[0030] In step two, use the pre-constructed hydraulic coupling element, adopt the unsaturated and unsteady seepage algorithm to process the geometric grid model to obtain the steady-state seepage field, and the method for extracting the pore pressure results is as follows:

[0031] Adopt the unsaturated and unsteady seepage algorithm, based on the permeability coefficient, relative permeability coefficients of liquid and gas phases, gravitational acceleration, and gradient function, to calculate the mass flow rate of the liquid phase and the mass flow rate of the gas phase.

[0032] According to the mass flow rate of the liquid phase and the mass flow rate of the gas phase, and based on the upstream and downstream water levels of the multi-layer soil slope during an earthquake, establish the flow conservation equation for porous media.

[0033] The flow conservation equation for porous media is used to calculate the liquid phase density and gas phase density according to the porosity, time, saturation, mass flow rate of the liquid phase, and mass flow rate of the gas phase.

[0034] Based on the water level and the liquid phase density and gas phase density, construct the upstream and downstream water pressure boundary conditions.

[0035] According to the upstream and downstream water pressure boundary conditions, obtain the steady-state seepage field and extract the pore pressure results.

[0036] As a preferred technical measure:

[0037] In step three, use the pre-constructed circular arc sliding surface calculation unit, and the method for determining the most dangerous circular arc sliding surface according to the pore pressure results is as follows:

[0038] In the first step, according to the geometric grid model, use the grid adaptive method to calculate the upper and lower limits of the radius search for the potential circular arc sliding surface, so as to obtain the upper and lower limits of the circular arc radius.

[0039] In the second step, according to the upper and lower limits of the circular arc radius, calculate the radius of each circular arc sliding surface, and generate the potential circular arc sliding surfaces of each combination of sliding-in points and sliding-out points.

[0040] In the third step, for the potential circular arc slip surface, according to the pore pressure results, the Bishop method is used to calculate the corresponding yield acceleration coefficient; and multiple potential slip surfaces and the corresponding yield acceleration coefficients are summarized to form an array of circular arc slip surfaces.

[0041] In the fourth step, based on the array of circular arc slip surfaces, the minimum value of the yield acceleration coefficient is obtained, and through the local bisection method, a search is carried out near the minimum value of the yield acceleration coefficient to find the most dangerous circular arc slip surface.

[0042] As a preferred technical measure:

[0043] In the first step, according to the geometric grid model, the upper and lower limits of the radius search for the potential circular arc slip surface are calculated. The method for obtaining the upper and lower limits of the circular arc radius is as follows:

[0044] Set the geometric control parameters of the potential circular arc slip surface, which include the coordinates of the entry point, the coordinates of the exit point, and the circular arc radius.

[0045] For any combination of the entry point and the exit point, select the circular arc radius at which the circular arc tangent at the entry point is perpendicular to the slope at the entry point as the lower limit of the circular arc radius.

[0046] According to the lower limit of the circular arc radius and the distance between the entry point and the exit point, calculate the maximum point-arc distance.

[0047] Based on the maximum point-arc distance and the model parameters, calculate the minimum point-arc distance.

[0048] According to the minimum point-arc distance, calculate the upper limit of the sliding arc radius.

[0049] Based on the grid characteristics, the upper limit of the sliding arc radius is iteratively corrected so that the slip surface generated by it maintains a distance from the slope to obtain the final upper limit of the circular arc radius.

[0050] Summarize the upper limit and the lower limit of the circular arc radius to obtain the upper and lower limits of the circular arc radius.

[0051] As a preferred technical measure:

[0052] In the third step, for the potential circular arc slip surface, the method for calculating the corresponding yield acceleration coefficient according to the pore pressure results using the Bishop method is as follows:

[0053] Project the pore pressure results onto the midpoints of the bottom edges of each soil strip to obtain the relevant parameters of the seepage pressure.

[0054] According to the relevant parameters of the seepage pressure, the anti-sliding force generated by the soil cohesion, the anti-sliding force generated by the normal stress of the slip surface, and the matrix suction generated by the unsaturated soil mass above the phreatic line, establish a calculation formula for the anti-sliding force of the soil strip.

[0055] Considering a multi-layer soil structure, calculate the equivalent cohesion at the bottom of a soil strip, the equivalent internal friction angle of a soil strip, together with the resultant force of hydrostatic pressure on a soil strip, the equivalent slope angle of a soil strip, the length of the bottom edge of a soil strip, and the effective saturation of the soil at the midpoint of the bottom edge of a soil strip, and input them into the soil strip anti-sliding force calculation formula to calculate the anti-sliding force of the soil strip considering seepage and mechanical loads on the soil strip;

[0056] Based on the Bishop method, according to the anti-sliding force of the soil strip, the inclination angle of the bottom edge of a soil strip, the distance from the centroid of a soil strip to the center of the slip arc, and the radius of the slip arc, obtain the peak horizontal seismic acceleration coefficient, that is, the yield acceleration coefficient.

[0057] As a preferred technical measure:

[0058] Step four, use the pre-constructed slope sliding surface calculation unit, based on the most dangerous circular arc sliding surface, solve to obtain the yield acceleration coefficients of several slope sliding surfaces, and according to the yield acceleration coefficients, obtain the method for determining the most dangerous slope sliding surface as follows:

[0059] Based on the most dangerous circular arc sliding surface, use the grid adaptive recursive algorithm to generate several sliding surfaces to ensure the integrity of the soil strips in the geometric grid model;

[0060] Using several sliding surfaces, bias the control parameters through a random algorithm to generate multiple potential slope sliding surfaces for solution;

[0061] For all potential slope sliding surfaces, solve to obtain the yield acceleration coefficients of several slope sliding surfaces by alternately iterating the explicit soil strip force and moment equilibrium equations;

[0062] Screen the yield acceleration coefficients of several slope sliding surfaces to obtain the minimum yield acceleration coefficient;

[0063] Take the slope sliding surface corresponding to the minimum yield acceleration coefficient as the most dangerous slope sliding surface.

[0064] As a preferred technical measure:

[0065] Step five, based on the pre-constructed acceleration calculation unit, use the local grid encryption algorithm to process the most dangerous slope sliding surface to obtain the yield acceleration data of a multi-layer soil slope as follows:

[0066] Step 51, for the most dangerous slope sliding surface, perform edge grid iterative encryption to obtain an encrypted grid;

[0067] Step 52, according to the encrypted grid, recalculate the yield acceleration coefficient and judge whether the yield acceleration coefficient converges. When the yield acceleration coefficient does not converge, execute Step 51; when the yield acceleration coefficient converges, execute Step 53;

[0068] Step 53: Convert the yield acceleration coefficient to the yield acceleration to obtain the yield acceleration data of a multi-layer soil slope.

[0069] To achieve one of the above purposes, the second technical solution of the present invention is as follows:

[0070] A method for calculating the yield acceleration of a multi-layer soil slope based on a grid, including the following:

[0071] Generate a geometric grid model according to the geological data of a multi-layer soil slope and couple the hydraulic parameters with the mechanical shear resistance information.

[0072] Process the geometric grid model to obtain a steady-state seepage field.

[0073] Determine the most dangerous circular arc slip surface according to the steady-state seepage field.

[0074] Based on the most dangerous circular arc slip surface, obtain the most dangerous slope slip surface.

[0075] Process the most dangerous slope slip surface to obtain the yield acceleration data of a multi-layer soil slope.

[0076] To achieve one of the above purposes, the third technical solution of the present invention is as follows:

[0077] A system for calculating the yield acceleration of a multi-layer soil slope based on hydraulic coupling, including:

[0078] One or more processors;

[0079] A storage device for storing one or more programs;

[0080] When the one or more programs are executed by the one or more processors, the one or more processors implement the above-mentioned method for calculating the yield acceleration of a multi-layer soil slope based on a grid.

[0081] Compared with the prior art solutions, the present invention has the following beneficial effects:

[0082] The present invention constructs a slope digital simulation unit, a hydraulic coupling unit, a circular arc slip surface calculation unit, the most dangerous slope slip surface, and an acceleration calculation unit to obtain the most dangerous slope slip surface and the corresponding yield acceleration data of a multi-layer soil slope. Therefore, it can fully consider the geological characteristics of the multi-layer soil slope structure to ensure that the obtained most dangerous slip surface is consistent with the actual situation; at the same time, considering the hydraulic coupling effect, it can accurately analyze the seismic slope stability of water retaining structures. Therefore, it can achieve accurate simulation calculation of the yield acceleration, and the solution is scientific, reasonable, and practical.

[0083] Furthermore, the present invention can calculate the yield acceleration of a multi-layer soil structure under complex working conditions, accurately obtain the most dangerous slope sliding surface and the corresponding yield acceleration data, with a wide range of applications and high efficiency. It is applicable to quickly evaluate the seismic slope stability of a multi-layer soil slope under the condition of hydro-mechanical coupling, and is convenient for the integrated simulation design of a multi-layer soil slope and the rapid iteration of the design scheme. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 FIG. is a schematic flow chart of a method for simulating and calculating the yield acceleration of a multi-layer soil slope of the present invention;

[0085] Figure 2 FIG. is a schematic diagram of a relationship curve between the yield acceleration coefficient and the slip arc radius of the present invention;

[0086] Figure 3 FIG. is a schematic diagram of a corresponding relationship between the isopotential line of the steady seepage pore pressure field and the most dangerous sliding surface of the slope when encountering an earthquake. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0087] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention, and are not used to limit the present invention.

[0088] On the contrary, the present invention covers any alternatives, modifications, equivalent methods and solutions made within the spirit and scope of the present invention defined by the claims. Further, in order to enable the public to better understand the present invention, some specific details are described in detail in the following detailed description of the present invention. Those skilled in the art can fully understand the present invention without the description of these details.

[0089] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the technical field to which the present invention belongs. The terms used herein are only for the purpose of describing specific embodiments, and are not intended to limit the present invention. The term "or / and" used herein includes any and all combinations of one or more of the related listed items.

[0090] As Figure 1 shown, the first specific embodiment of the method for calculating the yield acceleration of a multi-layer soil slope based on a grid of the present invention:

[0091] The method for calculating the yield acceleration of a multi-layer soil slope based on a grid includes the following steps:

[0092] Step 1, through a pre-constructed slope digital simulation unit, according to the geological data of a certain multi-layer soil slope, and coupling the hydraulic parameters and the mechanical shear resistance information, a geometric grid model is generated;

[0093] Step 2: Using the pre-constructed hydro-coupling unit and adopting the unsaturated and unsteady seepage algorithm, process the geometric grid model to obtain the steady seepage field and extract the pore pressure results;

[0094] Step 3: Using the pre-constructed circular arc slip surface calculation unit, determine the most dangerous circular arc slip surface according to the pore pressure results;

[0095] Step 4: Using the pre-constructed slope slip surface calculation unit, based on the most dangerous circular arc slip surface, solve to obtain the yield acceleration coefficients of several slope slip surfaces, and determine the most dangerous slope slip surface according to the yield acceleration coefficients;

[0096] Step 5: Based on the pre-constructed acceleration calculation unit, adopt the local grid refinement algorithm to process the most dangerous slope slip surface to obtain the yield acceleration data of a multi-layer soil slope.

[0097] The second specific embodiment of the grid-based yield acceleration calculation method for multi-layer soil slopes of the present invention:

[0098] The grid-based yield acceleration calculation method for multi-layer soil slopes includes the following:

[0099] According to the geological data of a multi-layer soil slope, and coupling the hydraulic parameters and mechanical shear resistance information, generate a geometric grid model;

[0100] Adopt the unsaturated and unsteady seepage algorithm to process the geometric grid model to obtain the steady seepage field and extract the pore pressure results;

[0101] According to the pore pressure results, determine the most dangerous circular arc slip surface and obtain the corresponding yield acceleration coefficient;

[0102] Based on the most dangerous circular arc slip surface and the yield acceleration coefficient, obtain the most dangerous slope slip surface and solve to obtain the corresponding yield acceleration;

[0103] Adopt the local grid refinement algorithm to process the most dangerous slope slip surface and the corresponding yield acceleration to obtain the yield acceleration data of a multi-layer soil slope.

[0104] A specific embodiment of applying the present invention to simulate and calculate a natural water-retaining slope in the western region:

[0105] A natural water-retaining slope in the west is composed of three main soil layers. The bottom and middle layers are clay with strong cohesion, while the upper layer has weak soil cohesion and a small friction angle in the middle layer. The multi-layer soil slope consists of two slope ratios, which are 2:1 and 2.5:1 from bottom to top. There is a wide platform with a width of 8 meters in the middle. The heights of the two slopes are 20m and 15m respectively, and the total height is 35m. The interfaces between the soil layers are approximately parallel and slightly inclined with a slope of 22:1.

[0106] The dynamic method needs to be used to analyze the seismic slope stability of the soil slope, calculate the permanent deformation of the slope sliding body. According to the method for solving the seismic yield acceleration of multi-layer soil hydraulic coupling based on discrete geometry proposed by the present invention, calculate the yield acceleration and its corresponding most dangerous sliding body for dynamic analysis, which specifically includes the following steps:

[0107] Step 1: According to the geotechnical structure design drawings and geological exploration data, establish a grid model including the zoning of each main material, define the hydraulic parameters and mechanical shear resistance performance of each area, and select the search area for the endpoints of the sliding surface.

[0108] Step 2: Adopt the unsaturated and unsteady seepage algorithm, and through the finite element method, calculate the steady seepage field of the geotechnical structure before encountering an earthquake, and extract the pore pressure results as the input conditions for the yield acceleration calculation.

[0109] Step 3: Adopt the improved simplified Bishop method, use the circular arc sliding surface generation and search algorithm based on discrete geometry to determine the most dangerous circular arc sliding surface, and solve the seismic yield acceleration coefficient, and use this circular arc sliding surface as the initial condition for Step 4.

[0110] Step 4: Adopt the improved Morgenstern-Price method, use the non-circular arc sliding surface generation and search algorithm based on discrete geometry to determine the most dangerous non-circular arc sliding surface, that is, the most dangerous slope sliding surface, and solve its corresponding seismic yield acceleration.

[0111] Step 5: For the most dangerous slope sliding surface, adopt the local grid encryption algorithm to iteratively correct the seismic yield acceleration.

[0112] In this embodiment, the present invention proposes a sliding surface generation algorithm considering the size of discrete boundary elements, which satisfies the kinematic rationality of the sliding surface and the integrity condition of soil strips while improving the efficiency and robustness of the sliding surface search algorithm. For the complex soil layer structure of the slope, the present invention proposes a layer intersection identification algorithm applicable to discrete geometry, which is used to calculate the equivalent shear strength parameters at the bottom of soil strips across soil layers, thereby improving the calculation accuracy of the yield acceleration. In terms of solving the yield acceleration, the present invention respectively derives and proposes a calculation format for directly solving the yield acceleration coefficient based on the simplified Bishop method and the Spencer method calculation frameworks, comprehensively utilizes the high solution efficiency of the circular arc method and the calculation accuracy of the non-circular arc method, and fully considers complex working conditions such as seepage pore pressure, unsaturated matrix suction, and external pressure load, so as to improve the solution efficiency. The method of the present invention has a wide application range and high efficiency, and is applicable to quickly evaluating the seismic slope stability of multi-layer soil slope structures under hydro-mechanical coupling conditions.

[0113] In this embodiment, step one specifically includes the following steps:

[0114] First step, obtain the partition characteristics of the main structural materials and the foundation stratification situation from the geotechnical structure design drawings and geological exploration reports. After removing the detailed structures, divide the first-order discrete grid units for each material area, requiring that the maximum grid size is less than 5% of the characteristic size of the area where it is located, and densify the line elements and surface elements around the dam slope to obtain a discrete geometric grid model suitable for seepage and other necessary finite element simulation analyses. The discrete geometric grid model contains a total of 4,605 nodes, 552 line elements and 8,890 first-order triangular plane elements, and the maximum element characteristic size of the slope line element and its adjacent area is 0.5 m.

[0115] Second step, obtain the regional hydraulic parameters, soil mechanics and shear strength parameters of the main structure and each material of the foundation from the geotechnical structure design drawings or safety appraisal reports, specifically including permeability coefficient, water retention curve, density, internal friction angle, cohesion, etc. Among them, the water retention curve is used to describe the relationship between the effective saturation and the pore pressure, and the Van Genuchten model is adopted, and its expression is as follows:

[0116]

[0117]

[0118]

[0119] Where represents the effective saturation, , are model parameters, is the pore pressure, is the air pressure, is the water pressure.

[0120] Through soil sampling for laboratory tests, the approximate Van Genuchten model parameters of each soil layer are calibrated, and the specific values are as follows:

[0121]

[0122]

[0123] The specific values of the permeability coefficient, density, internal friction angle, and cohesion of each soil layer are shown in Table 1.

[0124] Table 1: Hydraulic and mechanical material property parameters of each soil layer of a multi-layer soil slope

[0125]

[0126] In the third step, considering the changes in the slope ratio of the slope surface and the structural stratification, near the positions such as each level platform on the slope surface, the slope ratio mutation point, and the intersection of the layer and the slope surface, according to expert experience or the trial calculation results of the finite element strength reduction method under the action of seismic inertial force, the search areas for the sliding-in point and the sliding-out point of the slip surface are selected as the input for the calculation of the yield acceleration.

[0127] Since the internal friction angle of the middle layer is small and there is a large difference from the other two layers, under the action of the self-weight of the soil mass, the sliding-out point of the most dangerous slip surface may be located at the intersection of the interface between the middle layer and the lower layer and the slope surface, and the sliding-in point may be located at the top of the slope. Therefore, taking 20 m outside the slope angle as the origin, the abscissa of the sliding-out point of the slip surface is selected as 35 m, and the coordinate search range of the sliding-in point is between 110 m and 130 m.

[0128] In this embodiment, in the second step, according to the upstream and downstream water levels of the soil retaining geotechnical structure during an earthquake, based on the finite element method, the Newton-Raphson iteration is used to solve the porous medium flow conservation equation for calculation and solution, and the specific formula is as follows:

[0129]

[0130] Where and represent the liquid phase and gas phase densities respectively; represents the porosity, ignoring the influence of pore pressure on the porosity; represents time; represents the saturation; and represent the mass flow rates of the liquid phase and gas phase respectively, which are calculated according to Darcy's law.

[0131] The calculation formulas for the mass flow rates of the liquid phase and gas phase are as follows:

[0132]

[0133]

[0134] wherein represents the permeability coefficient, respectively represent the relative permeability coefficients of the liquid phase and the gas phase, represents the acceleration due to gravity, represents the gradient function, .

[0135] The upstream and downstream water pressure boundary conditions are expressed as:

[0136]

[0137] wherein represents the water level, represents the ordinate, represents taking the minimum value.

[0138] The upstream water level on the slope is set to 40 m, the water depth is 20 m, the downstream water depth is 0 m, and the groundwater level is 20 m. To correctly calculate the pore pressure of unsaturated soil, first calculate from the initial condition where the upstream and downstream water levels are both at the top of the slope. Taking 1 day as the time step, calculate the upstream and downstream water levels to drop to the stable water level within 3 days, and automatically extend the calculation time step until the relative change amount of the maximum pore pressure in the model between two adjacent time steps is less than 1%. Extract the steady-state pore pressure field and the gas phase pressure field from the seepage results of the last time step as the input conditions for calculating the yield acceleration.

[0139] In this embodiment, in the third step, the following steps are specifically included:

[0140] The first step is to calculate the upper and lower limits of the search for the radius of the potential circular arc slip surface for each combination of the sliding-in point and the sliding-out point according to the discrete geometry input, which includes the following content:

[0141] The circular arc slip surface can be determined by three geometric parameters: the coordinates of the sliding-in point, the coordinates of the sliding-out point, and the radius of the circular arc. Therefore, these three geometric parameters are defined as the control parameters of the circular arc slip surface.

[0142] For any combination of the sliding-in point and the sliding-out point, the method for selecting the lower limit of the circular arc radius is: select the circular arc radius where the tangent of the circular arc at the sliding-in point is perpendicular to the slope at that point as the lower limit of the radius, denoted as , if this circular arc has other intersections with the calculation domain boundary except the sliding-out point, then select the circular arc radius tangent to the relevant boundary as the lower limit.

[0143] For any combination of the sliding-in point and the sliding-out point, the method for selecting the upper limit of the circular arc radius is as follows:

[0144] Define The distance between the straight line passing through the sliding-in point and the sliding-out point and the tangent line of the arc parallel to it can be used as a control parameter for the arc sliding surface instead of the arc radius. Calculate the distance The minimum value of The formula is as follows:

[0145]

[0146] Where Represents The maximum value of, which is calculated from and is a model parameter, representing the number of potential sliding surfaces generated by each combination of sliding-in points and sliding-out points.

[0147] Calculate the distance The maximum value of The formula is:

[0148]

[0149] Where Represents half of the distance between the sliding-in point and the sliding-out point.

[0150] Upper limit of the sliding arc radius Is obtained by transforming The formula is as follows:

[0151]

[0152] Finally, iteratively correct the upper limit of the arc radius so that the sliding surface generated by it keeps a certain distance from the slope, and each soil strip of the sliding body contains at least one discrete geometric element. Specifically, in each iteration, first check whether the minimum distance between the sliding surface corresponding to the upper limit of the radius and the slope is greater than the minimum unit size of the discrete geometry. If it is not greater than the minimum size, then correct the upper limit of the radius As:

[0153]

[0154] Where And Respectively represent the upper limit of the radius in the (n + 1)-th and n-th iterations, Represents the minimum unit size of the discrete geometry. Secondly, check whether the vertical distance between the bottom endpoints of each soil strip and the slope is greater than the discrete geometric unit size at the corresponding slope. If it is not greater than the discrete geometric unit, then correct the ordinate As:

[0155]

[0156] Where respectively represent the vertical coordinates of the right endpoints of the bottom sides of the i-th soil strip in the (n + 1)-th and n-th iterations, represents the maximum size of the discrete geometric unit at the slope surface.

[0157] Meanwhile, regarding the sliding-in point and the sliding-out point, as well as the point the arc radius is used as the upper limit of the arc radius after the (n + 1)-th iteration correction, represents the abscissa of the sliding-in point, represents the width of the soil strip, is an arbitrary integer, represents the vertical coordinate of the intermediate node.

[0158] Step 2: Generate potential circular arc sliding surfaces for each combination of sliding-in points and sliding-out points, which includes the following:

[0159] Firstly, according to the upper limit of the corrected arc radius, calculate the minimum value of the distance The formula for is as follows:

[0160]

[0161] where represents the upper limit of the arc radius.

[0162] Calculate the radius of each circular arc sliding surface , and its calculation formula is:

[0163]

[0164]

[0165] Step 3: For the potential sliding surfaces of each combination of sliding-in points and sliding-out points, use the improved simplified Bishop method to calculate the yield acceleration coefficient, which includes the following:

[0166] Firstly, calculate the self-weight, hydrostatic pressure of each soil strip in the sliding body, and their moments about the center of the sliding arc. According to the material properties of each soil layer area input in Step 1, perform numerical integration within the soil strip to calculate the self-weight of the soil strip, and its calculation formula is as follows:

[0167]

[0168] where represents the self-weight of soil strip i, and represent the density and area of unit k, represents the gravitational acceleration.

[0169] Calculate the hydrostatic pressure of each soil strip and its moment about the center of the sliding arc by integrating on the one-dimensional element of the soil strip slope , and its calculation formula is as follows:

[0170]

[0171]

[0172] Among them represents the resultant vector of the hydrostatic pressure acting on soil strip i, represents the moment of the hydrostatic pressure acting on soil strip i about the center of the slip arc, represents the normal vector of the slope surface, represents the hydrostatic pressure distribution function, respectively represent the coordinate vectors of the acting point of the hydrostatic pressure and the center of the slip arc, respectively represent the distances from the starting points of the slope surfaces of soil strip i and soil strip i - 1, is the integration variable, is the two-dimensional vector product operator.

[0173] Secondly, since the input discrete geometric grid model only contains the spatial discrete information of the soil layer area and does not contain the geometric information of the soil layer interface, a recursive algorithm is required to identify the intersection points of the bottom edges of the soil strips across soil layers and the soil layer interface for calculating the equivalent shear strength parameters.

[0174] For the two ends of the bottom of soil strip i, judge whether the shear strength parameters of the discrete geometric elements where they are located are the same, that is, judge whether they belong to the same soil layer. If they do not belong to the same soil layer, take the midpoint of the bottom edge and judge whether it belongs to the same soil layer as the left and right endpoints respectively. The bottom edge of the soil strip is gradually bisected according to the above method until the distance between the two points participating in the discrimination is less than the predefined threshold, that is, the following inequality is satisfied in the nth layer of recursion:

[0175]

[0176] Among them represent the two coordinate points on the bottom edge of the soil strip participating in the soil layer parameter judgment in the nth layer of recursion, represents the distance threshold, and its value is 1 / 20 of the length of the bottom edge of the soil strip.

[0177] After the recursion ends, intersection points of the bottom edges of the soil strips and the soil layer interface are obtained, represents the number of soil layer areas spanned by the bottom edge of the soil strip. Then calculate the equivalent shear strength cohesion and the equivalent internal friction angle , and their calculation formulas are as follows:

[0178]

[0179]

[0180] Among them respectively represent the material cohesion and the internal friction angle of soil layer k spanned by soil slice i. represents the central angle of soil layer k within the bottom edge of soil slice i. Under the approximate assumption of uniform vertical normal stress distribution at the bottom edge of the soil slice, the above equivalent shear strength parameter algorithm makes the simplified numerical calculation results of the anti-sliding force and moment on the final slip surface close to the exact integration results, thus being more able to reflect the true limit equilibrium state of the sliding body.

[0181] The seepage pressure related parameters are obtained by projecting the pore pressure field and the gas phase pressure field in the finite element unsaturated seepage simulation calculation results onto the midpoint of the bottom edge of each soil slice. Specifically, for the midpoint of the bottom edge of each soil slice, search for the finite element where it is located in the input discrete geometric grid model, and calculate the pore pressure at this point using the element shape function. and the gas phase pressure , and their calculation formulas are as follows:

[0182]

[0183]

[0184] where represents the shape function of Gauss point k, which satisfies the value of 1 at Gauss point k and the value of 0 at other Gauss points. respectively represent the gas phase pressure and the pore pressure at Gauss point k of the element in the seepage finite element simulation results.

[0185] Finally, based on the simplified Bishop algorithm framework, the peak ground acceleration coefficient when the slope stability coefficient is 1 is derived as the yield acceleration coefficient , and its calculation formula is:

[0186]

[0187] where represents the yield acceleration coefficient, represents the inclination angle of the bottom edge of soil slice i, represents the distance from the centroid of soil slice i to the center of the slip arc, represents the radius of the slip arc, represents the moment of the hydrostatic pressure on soil slice i about the center of the slip arc, represents the anti-sliding force calculated according to the seepage effect and the vertical component of the mechanical load on soil slice i , and its specific expression is:

[0188]

[0189] where respectively represent the equivalent cohesion and the internal friction angle at the bottom of soil slice i, represents the resultant force of the hydrostatic pressure on soil slice i, Denote the equivalent slope angle of soil strip \(i\). Denote the length of the bottom of the soil strip. is the effective saturation of the soil mass at the midpoint of the bottom of soil strip \(i\). Denote the self-weight of soil strip \(i\).

[0190] In the expression, the first term is the anti-sliding force generated by soil cohesion, the second term is the anti-sliding force generated by the normal stress on the slip surface, and the third part is the matrix suction generated by the unsaturated soil mass above the phreatic line.

[0191] Effective saturation of soil mass The calculation formula is as follows:

[0192]

[0193] Fourthly, through the local binary search method, accurately search for the most dangerous circular arc slip surface and its yield acceleration coefficient. According to the combination of the slip-in point and the slip-out point where the minimum yield acceleration coefficient is obtained in the third step, the connection line between the slip-in and slip-out points, and the distance of the circular arc tangent parallel to the connection line between the slip-in and slip-out points , calculate the yield acceleration coefficient of the corresponding slip arc respectively according to the yield acceleration coefficient solution method described in the third step where:

[0194]

[0195] If a new minimum yield acceleration coefficient is generated, update and continue the iterative calculation the yield acceleration coefficient of the corresponding slip arc until no new minimum value is generated, where Denote the number of iterations, is the control parameter of the most dangerous slip arc updated in the \((k - 1)\)-th iteration.

[0196] Take the coordinates of the slip-in point, the coordinates of the slip-out point and the radius of the slip arc of the most dangerous slip surface as the initial conditions of the non-circular arc slip surface search algorithm in the fourth step.

[0197] In this specific embodiment, the number of soil strips divided is 10, and calculate the most dangerous slip arc radius and yield acceleration coefficient under the following four working conditions:

[0198] Working condition 1: Only seismic action;

[0199] Working condition 2: On the basis of working condition 1, consider the effect of seepage pore pressure;

[0200] Working condition 3: On the basis of working condition 2, consider the effect of hydrostatic pressure;

[0201] Working condition 4: On the basis of working condition 3, consider the unsaturated matrix suction of the soil mass.

[0202] The calculation results are shown in Table 2. It can be seen that in Case 2, if only the seepage pore pressure effect is considered, due to the high water level and the small self-weight of the sliding body, the yield acceleration coefficient is significantly lower than that in other cases, and the seismic instability risk of the slope is relatively high. After considering the hydrostatic pressure load in Case 3, the yield acceleration coefficient is slightly higher than that in Case 1, indicating that the stabilizing effect of water pressure on the slope is stronger than the destabilizing effect of the pore pressure at the current water level, which is beneficial to the slope stability. After further considering the matrix suction of unsaturated soil in Case 4, the yield acceleration coefficient is further increased, reflecting the strengthening effect of soil matrix suction on the shear strength, and it is consistent with the actual situation.

[0203] Table 2: Calculation results of the yield acceleration coefficient and the most dangerous slip arc radius for each case

[0204]

[0205] In Case 4, the equivalent shear strength algorithm of soil strips across soil layers proposed in the present invention is compared with the existing non-equivalent shear strength calculation method. The comparison results can be seen in Figure 2 . According to Figure 2 The relationship curve shown, the seismic yield acceleration coefficient and the slip arc radius calculated directly using the midpoint of the bottom edge of the soil strip by the non-equivalent shear strength calculation method result in a non-smooth relationship curve between the slip arc radius and the yield acceleration coefficient; while the method of the present invention eliminates the local oscillation of the yield acceleration coefficient caused by the discontinuity of the shear strength parameters of multi-layer soil materials, making the curve smoother and facilitating the accurate identification of the most dangerous slip surface in complex soil layer structures; the finally obtained most dangerous slip arc radius is 160.4 m, which is approximately tangent to the boundary between the middle and lower soil layers, indicating that during the earthquake, shear failure is likely to occur at the interface of the middle and lower soil layers, leading to slope sliding instability.

[0206] Based on the above calculation results, the yield acceleration solving method based on discrete geometry proposed in the present invention is applicable to the search for the most dangerous circular slip surface and the calculation of seismic yield acceleration. The calculation results are reasonable and applicable to complex geotechnical structures with multi-layers. The most dangerous circular slip surfaces obtained in Cases 1, 3, and 4 are the same, so this slip surface can be used as the input and initial condition for Step 4.

[0207] In this embodiment, Step 4 specifically includes the following steps:

[0208] First step, judge whether the minimum value of the yield acceleration coefficient in the current iteration step converges. The criterion is as follows:

[0209]

[0210] Where respectively represent the minimum values of the yield acceleration coefficient in the k-th and k-1-th iterations, It represents the convergence threshold, with a value of 0.002.

[0211] If it converges, input the control parameters of the most dangerous slip surface corresponding to the minimum yield acceleration coefficient into Step 5; otherwise, continue with the current step.

[0212] Step 2: Generate or screen potential non - circular slip surfaces, which include the following:

[0213] Select the number of soil strips , and define a total of control parameters for the non - circular slip surface. Their expressions are as follows:

[0214]

[0215] where and are the abscissas of the left and right endpoints of the slip surface, and is the ordinate of the middle node of the slip surface.

[0216] The soil strip dividing lines are vertical, and the horizontal widths of each soil strip are equal. Their calculation formula is as follows:

[0217] ,

[0218] According to the horizontal widths of each soil strip, calculate the abscissas of each middle node. Their calculation formula is as follows:

[0219]

[0220] If it is the first iteration, generate an initial non - circular slip surface set. One of the potential slip surfaces takes the circular slip surface obtained in Step 3, and the other slip surfaces are generated using a recursive algorithm to ensure the kinematic rationality of the slip surface and the integrity of the soil strips in the discrete geometric grid model.

[0221] Specifically, first generate a pair of coordinates of the slip - in point and the slip - out point in the search range of the slip - in point and the slip - out point according to a uniform distribution. Starting from the slip - out point, calculate the upper and lower limits of the ordinate of the middle node in turn. Select the lower limit such that is above the extension line of the bottom edge of the previous soil strip, and the angle between the left end of the slip surface and the slope is less than 45 degrees. Select the upper limit such that is on the line connecting and the slip - in point, and at the same time satisfy , where and respectively represent the vertical coordinate of the slope surface corresponding to the intermediate node \(i\) and the longitudinal unit size of the discrete geometric grid model at that location.

[0222] Check whether the upper and lower limits of the generated intermediate nodes are reasonable. If then return to the previous intermediate node and randomly generate again , otherwise randomly generate according to the uniform distribution, and its expression is as follows:

[0223]

[0224] If it is currently in the second and subsequent iterations, then screen sliding surfaces from all potential sliding surfaces in the previous iteration through the roulette algorithm to enter the next iteration, and the probability of each sliding surface being selected is:

[0225]

[0226] where respectively represent the maximum and minimum values of the yield acceleration coefficient in the previous iteration, represents the yield acceleration coefficient of the \(k\)-th sliding surface, represents the probability that the sliding surface \(k\) is selected in the roulette algorithm.

[0227] The third step is to generate multiple potential sliding surfaces for solution by randomly biasing the control parameters based on each sliding surface generated in the second step. For the sliding surface , the amplitude of the random bias is:

[0228]

[0229] where are model parameters, and the values are 120m, 100%, and 10% of the soil strip width respectively, is the maximum number of iterations, with a value of 50, is the current iteration number, represents the yield acceleration coefficient of the sliding surface \(i\), represents the yield acceleration coefficient of the sliding surface \(k\), represents the minimum yield acceleration coefficient, represents the maximum value function.

[0230] Furthermore, the number of sliding surfaces generated for solution is:

[0231]

[0232] where represents the maximum yield acceleration coefficient, Denotes the yield acceleration coefficient of the sliding surface k. Denotes the yield acceleration coefficient of the sliding surface i.

[0233] Where Is a model parameter with a value of 40.

[0234] The calculation formula for the randomly selected sliding surface control parameter is:

[0235]

[0236] Where Denotes the control parameter after offset, Denotes the control parameter before offset. Similarly, after randomly offsetting the i-th control parameter, starting from the (i + 1)-th control parameter, the subsequent control parameters are checked for legality and adjusted according to the recursive algorithm described in the second step.

[0237] Step 4: For all potential sliding surfaces generated in Step 3, based on the Morgenstern-Price calculation framework, solve for the yield acceleration coefficient.

[0238] First, calculate the self-weight and hydrostatic pressure of each soil strip of the sliding body according to the method described in the third step of Step 3. The calculation formula for the moment of the hydrostatic pressure on the bottom edge of each soil strip is:

[0239]

[0240] Where Denotes the moment of the hydrostatic pressure on the center of the bottom edge of soil strip i about the bottom edge center, Denotes the normal vector of the slope surface, Denotes the hydrostatic pressure distribution function, Denote the coordinates of the acting point of the hydrostatic pressure and the center of the bottom edge of soil strip i respectively, Denote the distances from the starting point of the slope surface of soil strip i and soil strip i - 1 respectively, Is the integration variable, Is the two-dimensional vector product operator.

[0241] Similarly, calculate the equivalent shear strength parameters of each soil strip according to the method described in the third step of Step 3, and obtain the pore pressure and gas phase pressure parameters through discrete geometric element projection from the finite element seepage simulation results.

[0242] Finally, based on the Morgenstern-Price calculation framework, derive the calculation formula for the yield acceleration coefficient under the non-circular arc assumption condition and solve for the yield acceleration coefficient. Specifically, when the slope stability coefficient is 1, according to the force and moment balance equations of the sliding body and the unconstrained boundary conditions, the calculation formula for the yield acceleration coefficient Is:

[0243]

[0244] wherein represents the shear force of the slip surface generated by the mechanical load other than the seismic action on soil strip i, represents the shear force generated by the horizontal seismic action, is an auxiliary variable, represents the number of soil strips.

[0245] The corresponding specific calculation formula is:

[0246]

[0247]

[0248]

[0249]

[0250] wherein respectively represent the equivalent cohesion and internal friction angle at the bottom of soil strip i, represents the resultant force of the hydrostatic pressure acting on soil strip i, represents the equivalent inclination angle of the slope of soil strip i, represents the length of the bottom side of the soil strip, represents the direction angle of the inter-strip force between soil strip i and soil strip i + 1. is the item related to seepage, is another auxiliary variable.

[0251] The corresponding specific expression is:

[0252]

[0253]

[0254] wherein respectively represent the gas phase pressure and pore pressure at the midpoint of the bottom side of soil strip i, represents the effective saturation at the midpoint of the bottom side of soil strip i. The effective saturation at the midpoint of the bottom side The calculation formula is as follows:

[0255]

[0256] According to the Morgenstern-Price hypothesis, the inter-strip direction angle is calculated from the inter-strip force coefficient and satisfies the following equation:

[0257]

[0258] wherein is the inter-strip force direction coefficient, To satisfy a smooth function with values at the sliding-in point and sliding-out point being the tangent values of the slope angles at the corresponding positions, it is a sine function with a period of twice the horizontal projection length of the sliding surface.

[0259] According to the moment equilibrium equations of each soil strip and the unconstrained conditions, the expression satisfied by the direction variable of the inter-strip force can be derived:

[0260]

[0261] where , represents the horizontal component of the inter-strip force , represents the width of the soil strip, represents the moment of the hydrostatic pressure about the center of the bottom of the soil strip, respectively represent the vertical coordinates of the midpoint of the bottom of soil strip i and the centroid of the soil strip, represents the self-weight of soil strip i. The recurrence calculation formula for the inter-strip force is:

[0262]

[0263] where is the auxiliary variable of soil strip i + 1.

[0264] Iteratively calculate the yield acceleration coefficient, the direction of the inter-strip force, and the direction coefficient of the inter-strip force alternately until the change amounts of the two in two adjacent iterations are less than the threshold so as to obtain the final result of the yield acceleration coefficient.

[0265] In this embodiment, in the fifth step, the following steps are specifically included:

[0266] First step, for the most dangerous sliding surface obtained in the fourth step, perform edge grid iterative encryption to eliminate the dependence of the yield acceleration result on the input discrete geometric dimensions. In each encryption process, calculate the distance from each unit node in the discrete geometry to the edge of the soil strip , and split the unit into two, and merge the unit with the adjacent unit, where is a model parameter with a value of 30% of the width of the soil strip.

[0267] Second step, use the grid after edge encryption and recalculate the yield acceleration coefficient by the method described in the fourth step. The convergence conditions are as follows:

[0268]

[0269] where are the yield acceleration coefficients obtained in the q-th and (q - 1)-th iterations respectively, is the convergence residual, with a value of 0.001. If it does not converge, return to the first step.

[0270] Finally, the yield acceleration coefficient is converted into the yield acceleration , and the final yield acceleration and the results of the related most dangerous slip surface are obtained. The calculation formula of the yield acceleration is as follows:

[0271]

[0272] where is the gravitational acceleration.

[0273] In this specific embodiment, the same number of soil strips as in step four is taken, that is, 10 soil strips. After 15 iterations of search, the most dangerous non-circular slip surface, that is, the most dangerous slope slip surface, is obtained. Subsequently, after 4 times of mesh refinement, the yield acceleration coefficient converges, and the final yield acceleration is obtained.

[0274] To further verify the calculation accuracy of the method proposed by the present invention, the present invention compares the calculation results of the specific embodiment with the results of the strength reduction method. The latter uses the finite element method with the same discrete geometric mesh model input. By gradually reducing the shear strength parameters of the soil layer material until a plastic penetration zone is formed, the calculation diverges; and then calculates the yield acceleration coefficient of the multi-layer soil structure under only seismic action. Its result can be regarded as an accurate solution, but it is less applied in the engineering field due to its low efficiency and complex calculation settings.

[0275] The calculation results of the yield acceleration coefficient of the two methods are shown in Table 3. It can be seen that compared with the yield acceleration coefficient calculated by the strength reduction method, the deviation of the yield acceleration coefficient in working condition 1 of the present invention is 1.2%, which proves that the direct solution method of the yield acceleration and the cross-layer equivalent shear strength algorithm of the present invention have high accuracy and are applicable to the rapid seismic stability analysis of complex multi-layer soil structures.

[0276] Table 3: Calculation results of yield acceleration

[0277]

[0278] At the same time, after comprehensively considering the hydrostatic pressure, seepage pore pressure and matrix suction of unsaturated soil in working condition 4 of the present invention, the calculated result of the yield acceleration coefficient is 0.215, that is, the yield acceleration is , which is significantly higher than the result of working condition 1 considering only seismic action, indicating that the hydrostatic pressure and matrix suction of the soil play a greater role in the seismic stability of the slope than the destabilizing effect of the seepage pore pressure, which is beneficial to the seismic stability of the slope. The yield acceleration coefficient results of the non-circular method in step four are similar to the results of the circular arc method in step three, forming a cross-validation. Furthermore, according to Figure 3The corresponding relationship between the equipotential lines of the steady seepage pore pressure field in [[]] and the most dangerous slip surface of the slope. Due to the existence of the weak interface between the middle and lower soil layers, it can be concluded that most of the slip surface is located near the interface, making the overall slip surface significantly non-circular, which reflects the necessity of Step 4.

[0279] In summary, the seismic yield acceleration simulation calculation method for multi-layer soil slopes based on discrete geometry proposed by the present invention is accurate and effective, and is applicable to the rapid seismic stability analysis of complex soil slopes.

[0280] An equipment embodiment applying the method of the present invention:

[0281] An electronic device, which includes:

[0282] One or more processors;

[0283] A storage device for storing one or more programs;

[0284] When the one or more programs are executed by the one or more processors, the one or more processors implement the above-mentioned grid-based multi-layer soil slope yield acceleration calculation method.

[0285] A computer medium embodiment applying the method of the present invention:

[0286] A computer-readable storage medium, on which a computer program is stored, and when the program is executed by a processor, it implements the above-mentioned grid-based multi-layer soil slope yield acceleration calculation method.

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

[0288] The present application is described according to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products of the embodiments of the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, and the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processors of general-purpose computers, special-purpose computers, embedded processors, or other programmable data processing devices to generate a machine, so that the instructions executed by the processors of the computer or other programmable data processing devices generate for implementing in the process Figure 1 a process or multiple processes or / and blocks Figure 1A device for the functions specified in one or more boxes.

[0289] 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 functions specified in one Figure 1 process or multiple processes or / and boxes Figure 1 or multiple boxes.

[0290] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to produce a computer-implemented process. Thus, the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one Figure 1 process or multiple processes or / and boxes Figure 1 or multiple boxes.

[0291] The unit in this application is an object that constitutes an objective description of the morphological structure by means of an entity or a virtual representation. The object is not equal to an object and is not limited to entities and virtuals. It can be a data processing function, a software program, a processing mode, a usage method, an operation method, a workflow, an application process, an electronic hardware, a circuit module, a processing system, a system imitation, or a simulation object.

[0292] 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 above embodiments, those of ordinary skill in the art can still modify or equivalently replace the specific implementation manners of the present invention. Any modification or equivalent replacement that does not depart from the spirit and scope of the present invention shall be covered by the protection scope of the claims of the present invention.

Claims

1. A grid-based multi-layer slope yield acceleration calculation method, characterized by: The following steps are involved: Step 1: Generate a geometric grid model through a pre-built slope digital simulation unit according to the geological data of a multi-soil layer slope and coupling hydraulic parameters and mechanical shear information; Step 2: Using the pre-built hydraulic coupling unit, the unsaturated and unsteady seepage algorithm, and the finite element method, the steady-state seepage field of the multi-layer slope before the earthquake is calculated, and the pore pressure results are extracted as the input conditions for the yield acceleration calculation; Step 3: Using the pre-built arc sliding surface calculation unit, determine the most dangerous arc sliding surface according to the pore pressure results; The method is as follows: The first step is to calculate the radius search upper and lower limits of the potential arc sliding surface based on the geometric mesh model using the mesh adaptive method to obtain the upper and lower limits of the arc radius; the method is as follows: Setting geometric control parameters of the potential arc sliding surface, including the sliding-in point coordinates, the sliding-out point coordinates and the arc radius; For any combination of sliding-in point and sliding-out point, the arc radius perpendicular to the slope surface at the sliding-in point is selected as the lower limit of the arc radius; Calculate the maximum point arc distance according to the lower limit of the arc radius and the distance between the sliding-in point and the sliding-out point; Based on the maximum point-arc distance and model parameters, the minimum point-arc distance is calculated; Calculate the upper limit of the sliding arc radius according to the minimum point arc distance; Based on the grid characteristics, the upper limit of the sliding arc radius is iteratively corrected so that the sliding surface generated by it keeps a distance from the slope surface, and the final upper limit of the arc radius is obtained; Sum up the upper limit and lower limit of the arc radius to obtain the upper and lower limits of the arc radius; The second step is to calculate the radius of each arc sliding surface according to the upper and lower limits of the arc radius, and generate the potential arc sliding surface of each combination of sliding-in point and sliding-out point; The third step is to calculate the corresponding yield acceleration coefficient of the potential arc sliding surface using the Bishop method according to the pore pressure results; and summarize multiple potential sliding surfaces and corresponding yield acceleration coefficients to form an array of arc sliding surfaces; The fourth step is to obtain the minimum value of the yield acceleration coefficient based on the circular arc sliding surface array, and search near the minimum value of the yield acceleration coefficient through local dichotomy to find the most dangerous circular arc sliding surface; Step 4: Using the pre-built slope sliding surface calculation unit, based on the most dangerous circular arc sliding surface, the improved Morgenstern-Price method is used to solve the yield acceleration coefficients of several slope sliding surfaces, and according to the yield acceleration coefficient, the most dangerous non-circular arc sliding surface generation and search algorithm based on discrete geometry is used to determine the most dangerous non-circular arc sliding surface, that is, the most dangerous slope sliding surface; Step 5: Based on the pre-built acceleration calculation unit, the local grid encryption algorithm is used to process the most dangerous slope sliding surface to obtain the yield acceleration data of a multi-soil layer slope.

2. The grid-based multi-layer slope yield acceleration calculation method according to claim 1, characterized in that: Step 1: The method for generating a geometric grid model by using a pre-built slope digital simulation unit, according to the geological data of a multi-layer slope, and coupling hydraulic parameters and mechanical shear information is as follows: Obtain geological data of a multi-layer slope, including structural material zoning characteristics, foundation stratification, and slope gradient changes; Determine the slope sliding position information based on the structural material zoning characteristics, foundation stratification and slope gradient changes; The slope sliding position information at least includes the position information of each level of the slope platform, the position information of the slope mutation, and the position information of the intersection of the layer and the slope; Based on the information of slope sliding position, the finite element strength reduction method is used to select the search area of ​​the sliding entry point and sliding exit point of the sliding surface; According to the search area, establish line elements and surface elements; The line elements and surface elements are encrypted to obtain the first-order discrete grid elements for the multi-layer slope; Obtain hydraulic parameters and mechanical shear information, including permeability coefficient, water retention curve, density, internal friction angle and cohesion; The water holding curve is used to describe the relationship between effective saturation and pore pressure; The hydraulic parameters and mechanical shear information are coupled with the first-order discrete grid units to form a geometric grid model.

3. The grid-based multi-layer slope yield acceleration calculation method according to claim 1, characterized in that: Step 2: Use the pre-built hydraulic coupling unit and the unsaturated and unsteady seepage algorithm to process the geometric grid model, obtain the steady-state seepage field, and extract the pore pressure results as follows: The unsaturated and unsteady seepage algorithm is used to calculate the mass flow rate of the liquid phase and the mass flow rate of the gas phase based on the permeability coefficient, the relative permeability coefficient of the liquid phase and the gas phase, the gravitational acceleration and the gradient function. According to the mass flow rate of the liquid phase and the mass flow rate of the gas phase, and based on the upstream and downstream water levels of the multi-layer slope when encountering an earthquake, the porous medium flow conservation equation is established; The porous media flow conservation equation is used to calculate the liquid phase density and gas phase density based on the porosity, time, saturation, mass flow rate of the liquid phase and the mass flow rate of the gas phase; Based on the water level and the liquid and gas phase densities, upstream and downstream water pressure boundary conditions are constructed; According to the upstream and downstream water pressure boundary conditions, the steady-state seepage field is obtained and the pore pressure results are extracted.

4. The grid-based multi-layer slope yield acceleration calculation method according to claim 1, characterized in that: The third step is to calculate the corresponding yield acceleration coefficient of the potential arc sliding surface based on the pore pressure results using the Bishop method as follows: Project the pore pressure results onto the midpoint of the bottom edge of each soil strip to obtain parameters related to seepage pressure; According to the parameters related to seepage pressure, the anti-sliding force generated by soil cohesion, the anti-sliding force generated by the normal stress of the sliding surface and the matrix suction generated by the unsaturated soil above the infiltration line, a calculation formula for the anti-sliding force of soil strips is established; Considering the multi-soil structure, the equivalent cohesion at the bottom of a soil strip and the equivalent internal friction angle of a soil strip are calculated, and the resultant force of the hydrostatic pressure on the soil strip, the equivalent inclination angle of the slope of the soil strip, the length of the bottom edge of the soil strip and the effective saturation of the soil at the midpoint of the bottom edge of the soil strip are input into the calculation formula of the soil strip anti-sliding force, and the soil strip anti-sliding force considering the seepage effect and the mechanical load on the soil strip is calculated; Based on the Bishop method, the horizontal earthquake peak acceleration coefficient, namely the yield acceleration coefficient, is obtained according to the anti-sliding force of the soil strip, the inclination angle of the bottom edge of a soil strip, the distance from the center of mass of a soil strip to the center of the sliding arc, and the radius of the sliding arc.

5. The grid-based multi-layer slope yield acceleration calculation method according to claim 1, characterized in that: Step 4: Use the pre-built slope sliding surface calculation unit to solve the yield acceleration coefficients of several slope sliding surfaces based on the most dangerous arc sliding surface, and obtain the method for determining the most dangerous slope sliding surface based on the yield acceleration coefficient as follows: Based on the most dangerous arc sliding surface, a mesh adaptive recursive algorithm is used to generate several sliding surfaces to ensure the integrity of the soil strips in the geometric mesh model. Using several sliding surfaces, the control parameters are biased by a random algorithm to generate multiple potential slope sliding surfaces for solution; For all potential slope sliding surfaces, the yield acceleration coefficients of several slope sliding surfaces are solved by alternately iterating the explicit soil strip force and moment equilibrium equations. Screening the yield acceleration coefficients of several slope sliding surfaces to obtain the minimum yield acceleration coefficient; The slope sliding surface corresponding to the minimum yield acceleration coefficient is taken as the most dangerous slope sliding surface.

6. The grid-based multi-layer slope yield acceleration calculation method according to claim 1, characterized in that: Step 5: Based on the pre-built acceleration calculation unit, the local grid encryption algorithm is used to process the most dangerous slope sliding surface, and the yield acceleration data of a multi-soil layer slope is obtained as follows: Step 51, for the most dangerous slope sliding surface, iteratively encrypt the edge grid to obtain an encrypted grid; Step 52, recalculating the yield acceleration coefficient according to the encrypted grid, and judging whether the yield acceleration coefficient converges. When the yield acceleration coefficient does not converge, executing step 51; when the yield acceleration coefficient converges, executing step 53; Step 53, converting the yield acceleration coefficient into yield acceleration to obtain yield acceleration data of a multi-soil layer slope.

7. The multi-layer slope yield acceleration calculation system based on hydraulic coupling is characterized by: It includes: one or more processors; A storage device for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the grid-based multi-soil layer slope yield acceleration calculation method as described in any one of claims 1-6.

Images