Multi-soil layer slope yield acceleration calculation method and system based on grids

Abstract

The invention discloses a grid-based multi-soil layer slope yield acceleration calculation method and system, and belongs to the technical field of geotechnical structure stability simulation. An existing geotechnical structure stability analysis scheme does not consider a multi-soil-layer structure, is not suitable for a multi-soil-layer side slope and cannot accurately calculate the earthquake yield acceleration. According to the grid-based multi-soil-layer slope yield acceleration calculation method, the most dangerous slope sliding surface of a certain multi-soil-layer slope and corresponding yield acceleration data are obtained by constructing a slope digital simulation unit, a hydraulic coupling unit, an arc sliding surface calculation unit, the most dangerous slope sliding surface and an acceleration calculation unit; therefore, the geological characteristics of the multi-soil-layer slope structure can be fully considered, and it is ensured that the obtained most dangerous sliding surface conforms to the actual situation. Meanwhile, the hydraulic coupling effect is considered, so that the anti-seismic slope stability of the water retaining building can be accurately analyzed, and accurate simulation calculation of the yield acceleration can be achieved.

Description

Technical Field

[0001] The invention relates to a grid-based multi-soil layer slope yield acceleration calculation method and system, belonging to the technical field of geotechnical structure stability simulation. Background Art

[0002] Earth-rock dams, levees, natural slopes and other geotechnical structures are generally composed of granular structures such as soil or sand and gravel, which are prone to sliding damage. External factors that cause soil slope sliding instability include increased soil pore pressure caused by seepage, top cracking, slope top ballast, sudden drop in water level, and earthquake action. However, due to the periodicity, short-term and randomness of earthquake 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 Code for Seismic Design of Hydraulic Structures in Hydropower Engineering NB / T35047-2015, it is necessary to calculate the seismic stability along the potential sliding surface and the permanent deformation of the dam body according to the effects of earthquake action. In this analysis process, before calculating the permanent deformation of the slope sliding body, it is necessary to accurately calculate the potential sliding surface of the slope (the most dangerous slope sliding surface) and the yield acceleration under earthquake conditions.

[0004] However, in the actual calculation process, since soil slopes are often composed of multiple soil layers and the physical and mechanical properties of each soil layer are discontinuous, the damage conditions of different soil layers are different, forming a complex sliding surface composed of multiple linear and circular arcs, which makes it difficult to accurately obtain the most dangerous sliding surface, and therefore it is impossible to accurately calculate the earthquake yield acceleration.

[0005] Furthermore, a Chinese patent application (publication number: CN116663274A) discloses a method for assessing the seismic stability of a prestressed anchor reinforced slope, comprising the following steps: Step 1, establishing a dynamic destruction mechanism of a prestressed anchor reinforced slope under earthquake action; Step 2, establishing an instantaneous displacement mode of a prestressed anchor reinforced slope under earthquake action, and obtaining a temporal recursive formula for anchor axial force, anchoring angle, slope height and slope angle; Step 3, deriving the change of yield acceleration over time based on the limit analysis upper bound method; Step 4, decoupling the coupling relationship between yield acceleration and seismic displacement based on a temporal analysis method, thereby obtaining the change of yield acceleration and seismic permanent displacement of the prestressed anchor reinforced slope; Step 5, using the obtained seismic permanent displacement to assess the seismic stability of the prestressed anchor reinforced slope.

[0006] The shortcomings of the above scheme are that the shape of the sliding surface is assumed to be a logarithmic spiral, and the multi-layer structure is not considered. It is not suitable for multi-layer slope scenarios. The most dangerous sliding surface obtained may deviate greatly from the actual situation, so it is still impossible to accurately calculate the seismic yield acceleration. In addition, the above scheme does not consider the hydraulic coupling effect, and has certain limitations in the seismic slope stability analysis of water retaining buildings.

[0007] The information disclosed in this Background Art is only for understanding the background of the inventive concept and therefore it may include information that does not constitute the prior art. Summary of the invention

[0008] In response to the above problem or one of the above problems, an object of the present invention is to provide a grid-based multi-soil layer slope yield acceleration calculation method and system. By constructing a slope digital simulation unit, a hydraulic coupling unit, an arc sliding surface calculation unit, the most dangerous slope sliding surface, and an acceleration calculation unit, the most dangerous slope sliding surface of a multi-soil layer slope and the corresponding yield acceleration data are obtained. 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, the hydraulic coupling effect is considered, so that the seismic slope stability of the water retaining building can be accurately analyzed, and thus accurate simulation calculation of the yield acceleration can be achieved.

[0009] In response to the above problem or one of the above problems, the second purpose of the present invention is to provide a grid-based multi-soil layer slope yield acceleration calculation method and system, which can realize the yield acceleration calculation of multi-soil layer structures under complex working conditions, and can accurately obtain the most dangerous slope sliding surface and corresponding yield acceleration data. It has a wide range of applications and high efficiency. It is suitable for quickly evaluating the seismic slope stability of multi-soil layer slope structures under hydraulic coupling conditions, and facilitates the integrated simulation design of multi-soil layer slopes and the rapid iteration of design schemes.

[0010] To achieve one of the above purposes, the first technical solution of the present invention is: The grid-based multi-layer slope yield acceleration calculation method includes the following steps: 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 and the unsaturated and unsteady seepage algorithm, the geometric grid model is processed to obtain the steady-state seepage field and extract the pore pressure results; Step 3: Using the pre-built arc sliding surface calculation unit, determine the most dangerous arc sliding surface according to the pore pressure results; Step 4: using the pre-built slope sliding surface calculation unit, based on the most dangerous arc sliding surface, solve the yield acceleration coefficients of several slope sliding surfaces, and determine the most dangerous slope sliding surface according to the yield acceleration coefficients; 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.

[0011] The present invention obtains the most dangerous slope sliding surface of a multi-soil layer slope and corresponding yield acceleration data by constructing a slope digital simulation unit, a hydraulic coupling unit, an arc sliding surface calculation unit, a most dangerous slope sliding surface, and an acceleration calculation unit. 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, the hydraulic coupling effect is considered, so that the seismic slope stability of the water retaining building can be accurately analyzed, and thus accurate simulation calculation of the yield acceleration can be achieved. The scheme is scientific, reasonable, and feasible.

[0012] Furthermore, the present invention can realize the yield acceleration calculation of multi-soil layer structures under complex working conditions, can accurately obtain the most dangerous slope sliding surface and the corresponding yield acceleration data, has a wide range of applications and high efficiency, is suitable for rapid evaluation of the seismic slope stability of multi-soil layer slope structures under hydraulic coupling conditions, and facilitates the integrated simulation design of multi-soil layer slopes and the rapid iteration of design schemes.

[0013] As the preferred technical measures: 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, water retention curve, density, internal friction angle and cohesion; the water retention 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.

[0014] As the preferred technical measures: 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, the 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.

[0015] As the preferred technical measures: Step 3: Using the pre-built arc sliding surface calculation unit, the method for determining the most dangerous arc sliding surface according to the pore pressure results 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, so as 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 array of arc sliding surfaces, and search near the minimum value of the yield acceleration coefficient through local dichotomy to find the most dangerous arc sliding surface.

[0016] As the preferred technical measures: The first step is to calculate the radius search upper and lower limits of the potential arc sliding surface according to the geometric mesh model. The method to obtain the upper and lower limits of the arc radius 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; The upper and lower limits of the arc radius are summed up to obtain the upper and lower limits of the arc radius.

[0017] As the preferred technical measures: 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.

[0018] As the preferred technical measures: 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.

[0019] As the preferred technical measures: 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.

[0020] To achieve one of the above purposes, the second technical solution of the present invention is: The grid-based multi-layer slope yield acceleration calculation method includes the following: According to the geological data of a multi-layer slope, the geometric grid model is generated by coupling the hydraulic parameters and mechanical shear information; The geometrical grid model is processed to obtain the steady-state seepage field; According to the steady-state seepage field, the most dangerous arc sliding surface is determined; Based on the most dangerous arc sliding surface, the most dangerous slope sliding surface is obtained; The most dangerous slope sliding surface is processed to obtain the yield acceleration data of a multi-layer slope.

[0021] To achieve one of the above purposes, the third technical solution of the present invention is: The multi-layer slope yield acceleration calculation system based on hydraulic coupling 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 above-mentioned grid-based multi-soil layer slope yield acceleration calculation method.

[0022] Compared with the prior art solutions, the present invention has the following beneficial effects: The present invention obtains the most dangerous slope sliding surface of a multi-soil layer slope and corresponding yield acceleration data by constructing a slope digital simulation unit, a hydraulic coupling unit, an arc sliding surface calculation unit, a most dangerous slope sliding surface, and an acceleration calculation unit. 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, the hydraulic coupling effect is considered, so that the seismic slope stability of the water retaining building can be accurately analyzed, and thus accurate simulation calculation of the yield acceleration can be achieved. The scheme is scientific, reasonable, and feasible.

[0023] Furthermore, the present invention can realize the yield acceleration calculation of multi-soil layer structures under complex working conditions, can accurately obtain the most dangerous slope sliding surface and the corresponding yield acceleration data, has a wide range of applications and high efficiency, is suitable for rapid evaluation of the seismic slope stability of multi-soil layer slope structures under hydraulic coupling conditions, and facilitates the integrated simulation design of multi-soil layer slopes and the rapid iteration of design schemes. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 A schematic diagram of a flow chart of a multi-soil layer slope yield acceleration simulation calculation method of the present invention; Figure 2 A schematic diagram of a relationship curve between the yield acceleration coefficient and the sliding arc radius of the present invention; Figure 3 This is a schematic diagram of the corresponding relationship between the contour lines of the stable seepage pore pressure field and the most dangerous sliding surface of the slope during an earthquake. DETAILED DESCRIPTION

[0025] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with 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 intended to limit the present invention.

[0026] On the contrary, the present invention covers any substitution, modification, equivalent method and scheme made on the essence and scope of the present invention as defined by the claims. Further, in order to make the public have a better understanding of 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.

[0027] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those commonly understood by those skilled in the art to which the present invention pertains. 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" as used herein includes any and all combinations of one or more of the associated listed items.

[0028] like Figure 1 As shown, the first specific embodiment of the grid-based multi-soil slope yield acceleration calculation method of the present invention is: The grid-based multi-layer slope yield acceleration calculation method includes the following steps: 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 and the unsaturated and unsteady seepage algorithm, the geometric grid model is processed to obtain the steady-state seepage field and extract the pore pressure results; Step 3: Using the pre-built arc sliding surface calculation unit, determine the most dangerous arc sliding surface according to the pore pressure results; Step 4: using the pre-built slope sliding surface calculation unit, based on the most dangerous arc sliding surface, solve the yield acceleration coefficients of several slope sliding surfaces, and determine the most dangerous slope sliding surface according to the yield acceleration coefficients; 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.

[0029] A second specific embodiment of the grid-based multi-layer slope yield acceleration calculation method of the present invention: The grid-based multi-layer slope yield acceleration calculation method includes the following: According to the geological data of a multi-layer slope, the geometric grid model is generated by coupling the hydraulic parameters and mechanical shear information; The unsaturated and unsteady seepage algorithm is used to process the geometric grid model, obtain the steady-state seepage field, and extract the pore pressure results; According to the pore pressure results, the most dangerous arc sliding surface is determined and the corresponding yield acceleration coefficient is obtained; Based on the most dangerous arc sliding surface and the yield acceleration coefficient, the most dangerous slope sliding surface is obtained, and the corresponding yield acceleration is solved; The local grid encryption algorithm is used to process the most dangerous slope sliding surface and the corresponding yield acceleration, and the yield acceleration data of a multi-layer slope is obtained.

[0030] A specific embodiment of applying the present invention to simulate a natural water retaining slope in the west: 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, the upper soil layer has weak cohesion, and the friction angle of the middle layer is small. The multi-soil slope consists of two sections with slope ratios of 2:1 and 2.5:1 from bottom to top, with a wide platform of 8 meters in the middle. The two sections are 20m and 15m high, with a total height of 35m. The interfaces of the soil layers are approximately parallel and slightly inclined, with a slope ratio of 22:1.

[0031] The soil slope needs to be analyzed for seismic slope stability using dynamic methods to calculate the permanent deformation of the slope sliding body. According to the discrete geometry-based multi-soil layer hydraulic coupling seismic yield acceleration solution method proposed in the present invention, the yield acceleration and its corresponding most dangerous sliding body are calculated for dynamic analysis, which specifically includes the following steps: Step 1: According to the geotechnical structure design drawings and geological survey data, a grid model containing the main material partitions is established, the hydraulic parameters and mechanical shear properties of each area are defined, and the sliding surface endpoint search area is selected.

[0032] Step 2: Use the unsaturated and unsteady seepage algorithm and the finite element method to calculate the steady-state seepage field of the geotechnical structure before the earthquake, and extract the pore pressure results as the input conditions for the yield acceleration calculation.

[0033] Step 3: Using the improved simplified Bishop method and the circular arc sliding surface generation and search algorithm based on discrete geometry, the most dangerous circular arc sliding surface is determined, and the earthquake yield acceleration coefficient is solved. The circular arc sliding surface is used as the initial condition for step 4.

[0034] Step 4: Using the improved Morgenstern-Price method and the non-circular sliding surface generation and search algorithm based on discrete geometry, the most dangerous non-circular sliding surface, that is, the most dangerous slope sliding surface, is determined, and its corresponding seismic yield acceleration is solved.

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

[0036] In this embodiment, the present invention proposes a sliding surface generation algorithm that takes into account the size of discrete boundary units, which satisfies the kinematic rationality of the sliding surface while also satisfying the soil strip integrity condition, thereby improving the efficiency and robustness of the sliding surface search algorithm. For complex soil layer structures of slopes, the present invention proposes a layer intersection identification algorithm suitable for discrete geometry, which is used to calculate the equivalent shear strength parameters of the bottom of the cross-soil strip, thereby improving the calculation accuracy of the yield acceleration. In terms of yield acceleration solution, the present invention derives and proposes a calculation format for directly solving the yield acceleration coefficient in one go based on the simplified Bishop method and Spencer method calculation frameworks, comprehensively utilizing the high solution efficiency of the arc method and the calculation accuracy of the non-arc method, and fully considering complex working conditions such as seepage pore pressure, unsaturated matrix suction and external pressure loads, thereby improving the solution efficiency. The method of the present invention has a wide range of applications and high efficiency, and is suitable for rapid evaluation of seismic slope stability of multi-soil layer slope structures under hydraulic coupling conditions.

[0037] In this embodiment, the step 1 specifically includes the following steps: The first step is to obtain the main structure material zoning characteristics and foundation stratification from the geotechnical structure design drawings and geological survey reports. After removing the detailed structure, the first-order discrete grid units are divided for each material area, requiring the maximum grid size to be less than 5% of the characteristic size of the area, and the line units and surface units are encrypted around the dam slope to obtain a discrete geometric grid model suitable for seepage and other necessary finite element simulation analysis. The discrete geometric grid model contains a total of 4605 nodes, 552 line units and 8890 first-order triangular plane units. The maximum unit characteristic size of the slope line unit and its adjacent area is 0.5m.

[0038] The second step is to obtain the regional hydraulic parameters, soil mechanics and shear strength parameters of the main structure and foundation materials from the geotechnical structure design drawings or safety appraisal reports, including permeability coefficient, water retention curve, density, internal friction angle and cohesion, etc. The water retention curve is used to describe the relationship between effective saturation and pore pressure, and the van Gruchteng model is used, which is expressed as follows:

[0039]

[0040]

[0041] in represents the effective saturation, , are model parameters, is the pore pressure, is the air pressure, For water pressure.

[0042] By taking soil samples for indoor tests, the parameters of the approximate Van Gruchteng model of each soil layer were calibrated, and the specific values ​​are as follows:

[0043]

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

[0045] Table 1: Hydraulic and mechanical material performance parameters of each soil layer in multi-layer slope

[0046] The third step is to consider the slope gradient changes and structural stratification. At various levels of slope platforms, slope mutation points, and intersections between layers and slopes, based on expert experience or the results of finite element strength reduction method calculations under the action of seismic inertial force, select the search areas for sliding entry and exit points of the sliding surface as input for yield acceleration calculation.

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

[0048] In this embodiment, in step 2, according to the upstream and downstream water levels of the retaining geotechnical structure when encountering an earthquake, the Newton-Raphson iteration is used to solve the porous medium flow conservation equation based on the finite element method to calculate and solve the porous medium flow conservation equation. The specific formula is as follows:

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

[0050] The mass flow rates of the liquid and gas phases are calculated as follows:

[0051]

[0052] in represents the permeability coefficient, represent the relative permeability coefficients of liquid and gas phases, respectively. represents the acceleration due to gravity, represents the gradient function, .

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

[0054] in Indicates the water level, represents the vertical coordinate, Indicates taking the minimum value.

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

[0056] In this embodiment, the step three specifically includes the following steps: The first step is to calculate the search upper and lower limits of the potential arc sliding surface radius of each combination of sliding-in point and sliding-out point based on the discrete geometric input, which includes the following: The arc sliding surface can be determined by three geometric parameters: the coordinates of the sliding-in point, the coordinates of the sliding-out point and the arc radius. Therefore, these three geometric parameters are defined as the arc sliding surface control parameters.

[0057] For any combination of sliding-in point and sliding-out point, the method for selecting the lower limit of the arc radius is: select the arc radius at the sliding-in point where the tangent line of the arc is perpendicular to the slope surface as the lower limit of the radius, recorded as If there are other intersections between this arc and the boundary of the computational domain except the slide-out point, the radius of the arc tangent to the relevant boundary is selected as the lower limit.

[0058] For any combination of slide-in and slide-out points, the method for selecting the upper limit of the arc radius is as follows: definition It is the distance between the straight line passing through the sliding-in point and the sliding-out point and the arc tangent line parallel to it, which can replace the arc radius as the arc sliding surface control parameter. Minimum value of The formula is as follows:

[0059] in express The maximum value of Calculated, is a model parameter, which indicates the number of potential sliding surfaces generated by each combination of sliding-in point and sliding-out point.

[0060] Calculating distance The maximum value of The formula is:

[0061] in Indicates half the distance between the slide-in point and the slide-out point.

[0062] Upper limit of sliding arc radius Depend on The conversion formula is as follows:

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

[0064] in and They represent the upper limits of the radius for the n+1th and nth iterations respectively, Indicates the minimum unit size of discrete geometry. Secondly, check whether the vertical distance between the bottom end of each soil strip and the slope surface is greater than the discrete geometry unit size at the corresponding slope surface. If it is not greater than the discrete geometry unit, correct the ordinate for:

[0065] in denote the right ordinate of the bottom edge of the ith soil strip in the n+1th and nth iterations, respectively. Indicates the maximum size of discrete geometric units on the slope.

[0066] At the same time, the slide-in point and slide-out point, as well as the point The arc radius is used as the upper limit of the arc radius after the n+1th iteration correction. Indicates the horizontal coordinate of the slide-in point. Indicates the width of the soil strip, is any integer, Indicates the ordinate of the middle node.

[0067] The second step is to generate the potential arc sliding surface of each combination of sliding-in point and sliding-out point, which includes the following contents: First, calculate the distance based on the corrected upper limit of the arc radius Minimum value of The formula is as follows:

[0068] in Indicates the upper limit of the arc radius.

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

[0070]

[0071] The third step is to calculate the yield acceleration coefficient for each potential sliding surface of the combination of sliding-in point and sliding-out point using the improved simplified Bishop method, which includes the following: First, calculate the deadweight, hydrostatic pressure, and moment about the center of the sliding arc of each soil strip of the sliding body. According to the material properties of each soil layer area input in step 1, perform numerical integration within the soil strip to calculate the deadweight of the soil strip. The calculation formula is as follows:

[0072] in represents the self-weight of soil strip i, and represents the density and area of ​​unit k, Represents the acceleration due to gravity.

[0073] By integrating over the one-dimensional unit of the soil strip slope, the hydrostatic pressure on each soil strip and its moment about the center of the sliding arc are calculated. , and its calculation formula is as follows:

[0074]

[0075] in represents the resultant force vector of the hydrostatic pressure on soil strip i, represents the moment of the hydrostatic pressure on soil strip i about the center of the sliding arc, represents the slope normal vector, represents the hydrostatic pressure distribution function, They represent the coordinate vectors of the hydrostatic pressure point and the sliding arc center, respectively. Respectively represent the distances from the starting point of the slope of strip i and strip i-1, is the integration variable, is the two-dimensional vector product operator.

[0076] Secondly, since the input discrete geometric grid model only contains the spatial discrete information of the soil layer area but does not contain the geometric information of the soil layer interface, a recursive algorithm is required to identify the intersection of the bottom edge of the cross-soil strip and the soil layer interface for calculating the equivalent shear strength parameters.

[0077] For the two ends of the bottom of soil strip i, determine whether the shear strength parameters of the discrete geometric units are the same, that is, determine 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 determine whether it belongs to the same soil layer as the left and right endpoints. Divide the bottom edge of the soil strip step by step according to the above method until the distance between the two points involved in the judgment is less than the predefined threshold, that is, the following inequality is satisfied in the nth recursive layer:

[0078] in Indicates the two coordinate points located at the bottom of the soil strip that recursively participate in the soil layer parameter judgment of the nth layer, Indicates the distance threshold, which is 1 / 20 of the bottom length of the soil strip.

[0079] After the recursion is completed, The intersection point of the bottom edge of a soil strip and the soil layer interface, Indicates the number of soil layers spanned by the bottom edge of the soil strip. Then calculate the equivalent shear strength and cohesion and the equivalent internal friction angle , and its calculation formula is as follows:

[0080]

[0081] in denote the material cohesion and internal friction angle of the soil layer k spanned by soil strip i, represents the central angle of soil layer k within the bottom edge of soil strip i. Under the approximate assumption that the vertical normal stress distribution at the bottom edge of the soil strip is uniform, the above equivalent shear strength parameter algorithm makes the simplified numerical calculation results of the anti-sliding force and moment on the final sliding surface close to the exact integration results, thus better reflecting the true limit equilibrium state of the sliding body.

[0082] The parameters related to seepage pressure are obtained by projecting the pore pressure field and gas phase pressure field in the finite element unsaturated seepage simulation calculation results to the midpoint of the bottom edge of each soil strip. Specifically, for the midpoint of the bottom edge of each soil strip, the finite element where it is located is searched in the input discrete geometric grid model, and the pore pressure at this point is calculated using the element shape function. Gas phase pressure , and its calculation formula is as follows:

[0083]

[0084] in Represents the shape function of Gaussian point k, which satisfies the value of 1 at Gaussian point k and the value of 0 at other Gaussian points. They respectively represent the gas phase pressure and pore pressure at the unit Gaussian point k in the seepage finite element simulation results.

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

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

[0087] in denote the equivalent cohesion and internal friction angle at the bottom of soil strip i, respectively. represents the resultant hydrostatic pressure on soil strip i, represents the equivalent inclination angle of the slope surface of soil strip i, Indicates the length of the bottom side of the soil strip, is the effective saturation of soil at the midpoint of the bottom edge of soil strip i, Indicates the self-weight of soil bar i.

[0088] In the expression, the first term is the anti-sliding force generated by the cohesion of the soil, the second term is the anti-sliding force generated by the normal stress of the sliding surface, and the third part is the matrix suction generated by the unsaturated soil above the infiltration line.

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

[0090] The fourth step is to accurately search for the most dangerous arc sliding surface and its yield acceleration coefficient through the local binary search method. According to the combination of the sliding-in point and the sliding-out point where the minimum yield acceleration coefficient is obtained in the third step, the line connecting the sliding-in and sliding-out points, and the arc tangent distance parallel to the line connecting the sliding-in and sliding-out points , respectively calculated according to the yield acceleration coefficient solution method described in the third step The yield acceleration coefficient of the corresponding sliding arc, where:

[0091] If a new minimum yield acceleration coefficient is generated, update And continue to iterate The yield acceleration coefficient of the corresponding sliding arc is calculated until no new minimum value is generated, where represents the number of iterations, The most dangerous arc sliding control parameters updated for the k-1th iteration.

[0092] The coordinates of the sliding entry point, the sliding exit point and the sliding arc radius of the most dangerous sliding surface are used as the initial conditions of the non-circular sliding surface search algorithm in step 4.

[0093] In this specific embodiment, the number of soil strips is divided into 10, and the most dangerous sliding arc radius and yield acceleration coefficient of the following four working conditions are calculated: Condition 1: earthquake action only; Working condition 2: Based on working condition 1, the effect of seepage pore pressure is considered; Working condition 3: Based on working condition 2, the hydrostatic pressure is considered; Working condition 4: Based on working condition 3, the unsaturated matrix suction of soil is considered.

[0094] The calculation results are shown in Table 2. It can be seen that if only the seepage pore pressure effect is considered in working condition 2, the yield acceleration coefficient is significantly reduced compared with other working conditions due to the high water level and small weight of the sliding body, and the risk of slope seismic instability is higher. After considering the hydrostatic pressure load in working condition 3, the yield acceleration coefficient is slightly higher than that in working condition 1, indicating that the stabilizing effect of water pressure on the slope is stronger than the instability effect of the current water level pore pressure, which is conducive to slope stability. After further considering the matrix suction of the unsaturated soil in working condition 4, the yield acceleration coefficient is further improved, reflecting the strengthening effect of the soil matrix suction on the shear strength, which is consistent with the actual situation.

[0095] Table 2: Calculation results of yield acceleration coefficient and most dangerous sliding arc radius under various working conditions

[0096] In working condition 4, the equivalent shear strength algorithm of the cross-soil strip proposed in the present invention is compared with the non-equivalent shear strength calculation method in the prior art. The comparison results can be seen in Figure 2 .according to Figure 2It can be seen from the relationship curve shown that the earthquake yield acceleration coefficient and the sliding arc radius calculated directly using the midpoint of the bottom edge of the soil strip without an equivalent shear strength calculation method result in a relationship curve between the sliding arc radius and the yield acceleration coefficient that is not smooth enough; 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 accurate identification of the most dangerous sliding surface in a complex soil layer structure; the most dangerous sliding arc radius finally obtained is 160.4m, which is approximately tangent to the boundary line between the middle and lower soil layers, indicating that during an earthquake, shear failure is prone to occur at the interface between the middle and lower soil layers, leading to slope sliding instability.

[0097] Based on the above calculation results, the yield acceleration solution method based on discrete geometry proposed in the present invention is suitable for searching the most dangerous arc sliding surface and calculating earthquake yield acceleration. The calculation results are reasonable and applicable to complex geotechnical structures with multiple soil layers. The most dangerous arc sliding surface obtained in working conditions 1, 3, and 4 is the same, so the sliding surface can be used as the input and initial condition of step 4.

[0098] In this embodiment, the step 4 specifically includes the following steps: The first step is to determine whether the minimum value of the yield acceleration coefficient of the current iteration step converges. The criteria are as follows:

[0099] in They represent the minimum values ​​of the yield acceleration coefficients of the kth and k-1th iterations, respectively. Indicates the convergence threshold, and its value is 0.002.

[0100] If convergence occurs, the most dangerous sliding surface control parameter corresponding to the minimum yield acceleration coefficient is input into step five, otherwise continue with the current step.

[0101] Step 2: Generate or filter Potential non-circular sliding surfaces include the following: Select the number of soil strips , define the control parameters of the non-circular sliding surface The expression is as follows:

[0102] in and is the horizontal coordinate of the left and right endpoints of the sliding surface, is the ordinate of the middle node of the sliding surface.

[0103] The strip dividing line is vertical, and the horizontal width of each strip is The calculation formula is as follows: , According to the horizontal width of each soil strip , calculate the horizontal coordinates of each intermediate node , and its calculation formula is as follows:

[0104] If the current iteration is the first one, an initial set of non-circular sliding surfaces is generated, in which one potential sliding surface is the circular sliding surface obtained in step 3, and the other sliding surfaces are generated using a recursive algorithm to ensure the kinematic rationality of the sliding surface and the integrity of the soil strips in the discrete geometric grid model.

[0105] Specifically, firstly, a pair of slide-in point coordinates and slide-out point coordinates are generated in the slide-in point and slide-out point search range according to a uniform distribution. Starting from the slide-out point, the vertical coordinates of the intermediate nodes are calculated in sequence. Upper and lower limits, select lower limit Make If the sliding surface is located above the extension line of the bottom edge of the previous soil strip and the angle between the left end of the sliding surface and the slope is less than 45 degrees, select the upper limit. Make lie in On the line connected to the slide-in point, ,in and They respectively represent the ordinate of the slope corresponding to the middle node i and the longitudinal unit size of the discrete geometric grid model at that location.

[0106] Check whether the upper and lower limits of the generated intermediate nodes are reasonable. Then return to the previous intermediate node and re-randomly generate , otherwise it is randomly generated according to a uniform distribution, and its expression is as follows:

[0107] If the current iteration is the second or subsequent iteration, the roulette algorithm is used to select all potential sliding surfaces in the previous iteration. sliding surfaces enter the next iteration, and the probability of each sliding surface being selected for:

[0108] in Respectively represent the maximum and minimum values ​​of the yield acceleration coefficient in the previous iteration, represents the yield acceleration coefficient of the kth sliding surface, Represents the probability of sliding surface k being selected in the roulette algorithm.

[0109] In the third step, based on each sliding surface generated in the second step, multiple potential sliding surfaces for solution are generated by randomly biasing the control parameters. , random bias amplitude for:

[0110] in are model parameters, with values ​​of 120m, 100% and 10% of the strip width respectively. is the maximum number of iterations, the value is 50, is the current iteration number, represents the yield acceleration coefficient of sliding surface i, represents the yield acceleration factor of the sliding surface k, represents the minimum yield acceleration factor, Represents the maximum value function.

[0111] Then, the number of sliding surfaces generated for solution for:

[0112] in represents the maximum yield acceleration coefficient, represents the yield acceleration factor of the sliding surface k, Represents the yield acceleration coefficient of sliding surface i.

[0113] in is a model parameter and its value is 40.

[0114] The calculation formula of the randomly selected sliding surface control parameters is:

[0115] in represents the control parameter after bias, Represents the control parameter before biasing. Similarly, after the i-th control parameter is randomly biased, 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.

[0116] In the fourth step, for all potential sliding surfaces generated in the third step, the yield acceleration coefficient is solved based on the Morgenstern-Price calculation framework.

[0117] First, calculate the deadweight 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 hydrostatic pressure on each soil strip about the bottom edge of the soil strip is:

[0118] in represents the moment of the hydrostatic pressure on soil strip i about the center of the bottom edge of the soil strip, represents the slope normal vector, represents the hydrostatic pressure distribution function, denote the hydrostatic pressure point and the coordinate vector of the center of the bottom edge of soil strip i, respectively. Respectively represent the distances from the starting point of the slope of strip i and strip i-1, is the integration variable, is the two-dimensional vector product operator.

[0119] Similarly, the equivalent shear strength parameters of each soil strip are calculated according to the method described in the third step of step three, and the pore pressure and gas phase pressure parameters are obtained from the finite element seepage simulation results through discrete geometric unit projection.

[0120] Finally, based on the Morgenstern-Price calculation framework, the yield acceleration coefficient calculation formula under the non-circular arc assumption is derived and the yield acceleration coefficient is solved. Specifically, when the slope stability coefficient is 1, the yield acceleration coefficient can be derived according to the sliding body force and moment equilibrium equation and the unconstrained boundary condition: The calculation formula is:

[0121] in represents the shear force on the sliding surface of soil strip i caused by mechanical loads other than earthquake action, represents the shear force generated by horizontal earthquake action, is an auxiliary variable, Indicates the number of soil strips.

[0122] The corresponding specific calculation formula is:

[0123]

[0124]

[0125]

[0126] in denote the equivalent cohesion and internal friction angle at the bottom of soil strip i, respectively. represents the resultant hydrostatic pressure on soil strip i, represents the equivalent inclination angle of the slope surface of soil strip i, Indicates the length of the bottom side of the soil strip, It represents the direction angle of the inter-strip force between soil strip i and soil strip i+1. is the seepage related term, is another auxiliary variable.

[0127] The corresponding specific expression is:

[0128]

[0129] in denote the gas phase pressure and pore pressure at the midpoint of the bottom edge of soil strip i, respectively. Indicates the effective saturation of the midpoint of the bottom edge of soil strip i. The calculation formula is as follows:

[0130] According to the Morgenstern-Price hypothesis, the angle between the strips It is calculated from the inter-bar force coefficient and satisfies the following equation:

[0131] in is the inter-strip force direction coefficient, In order to satisfy the requirement that the values ​​of the sliding-in point and sliding-out point are smooth functions of the tangent value of the slope angle at the corresponding position, It is a sinusoidal function with a period of twice the length of the horizontal projection of the sliding surface.

[0132] According to the moment equilibrium equation of each soil strip and the unconstrained condition, the inter-strip force direction variable can be derived The expression satisfied is:

[0133] in , Indicates inter-strike force The horizontal component of Indicates the width of the soil strip, represents the moment of hydrostatic pressure about the center of the bottom edge of the soil strip, They represent the ordinates of the midpoint of the bottom edge of soil strip i and the centroid of soil strip i, respectively. represents the deadweight of soil strip i. The recursive calculation formula of inter-strip force is:

[0134] in It is the auxiliary variable of soil strip i+1.

[0135] Alternately iterate and calculate the inter-strip force direction and inter-strip force direction coefficient of the yield acceleration coefficient until the change between the two adjacent iterations is less than the threshold Thus, the final yield acceleration coefficient result is obtained.

[0136] In this embodiment, the step 5 specifically includes the following steps: In the first step, the edge mesh is iteratively refined for the most dangerous sliding surface obtained in step 4 to eliminate the dependence of the yield acceleration result on the input discrete geometry size. During each step of the refinement process, the distance from each unit node in the discrete geometry to the edge of the soil strip is calculated. ,Will The unit is split into two, The cells of are merged with the adjacent cells, where is a model parameter, and its value is 30% of the width of the soil strip.

[0137] In the second step, the edge-encrypted mesh is used to recalculate the yield acceleration coefficient using the method described in step 4. The convergence criteria are as follows:

[0138] in are the yield acceleration coefficients obtained from the qth and q-1th iterations, respectively. is the convergence residual, and its value is 0.001. If it does not converge, return to the first step.

[0139] The yield acceleration coefficient Convert to yield acceleration , and obtain the final yield acceleration and the most dangerous sliding surface results. The calculation formula is as follows:

[0140] in is the acceleration due to gravity.

[0141] In this specific embodiment, the same number of soil strips as in step 4, i.e., 10 soil strips, is used. After 15 iterations of searching, the most dangerous non-circular sliding surface, i.e., the most dangerous slope sliding surface, is obtained. After 4 mesh densifications, the yield acceleration coefficient converges to obtain the final yield acceleration.

[0142] In order to further verify the calculation accuracy of the method proposed in 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 and the same discrete geometric grid model input. The shear strength parameters of the soil material are gradually reduced until a plastic penetration zone is formed, which leads to calculation divergence; then the yield acceleration coefficient of the multi-soil layer structure under earthquake action alone is calculated. The result can be regarded as an exact solution, but it is rarely used in the engineering field due to its low efficiency and complex calculation settings.

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

[0144] Table 3: Yield acceleration calculation results

[0145] At the same time, after comprehensively considering the hydrostatic pressure, seepage pore pressure and unsaturated soil matrix suction, the yield acceleration coefficient of working condition 4 of the present invention is calculated to be 0.215, that is, the yield acceleration is , which is significantly higher than the result of condition 1 considering only the earthquake effect, indicating that the hydrostatic pressure and soil matrix suction have a greater effect on the seismic stability of the slope than the destabilizing effect of the seepage pore pressure, which is beneficial to the stability of the seismic slope. The yield acceleration coefficient result of the non-circular arc method in step 4 is similar to the result of the circular arc method in step 3, forming a cross-validation. Then, according to Figure 3 The correspondence between the stable seepage pore pressure field contour lines and the most dangerous sliding surface of the slope, and due to the existence of the weak interface between the middle and lower soil layers, it can be concluded that most of the sliding surface is located near the interface, making the overall sliding surface significantly non-circular, reflecting the necessity of step four.

[0146] In summary, the discrete geometry-based multi-soil slope seismic yield acceleration simulation calculation method proposed in the present invention is accurate and effective, and is suitable for rapid seismic stability analysis of complex soil layer slopes.

[0147] An embodiment of a device applying the method of the present invention: An electronic device comprising: 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 above-mentioned grid-based multi-soil layer slope yield acceleration calculation method.

[0148] A computer medium embodiment using the method of the present invention: A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the above-mentioned grid-based multi-soil layer slope yield acceleration calculation method.

[0149] 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 take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application 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, optical storage, etc.) containing computer-usable program codes.

[0150] The present application is described by flowcharts or / and 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 or / and block in the flowchart or / and block diagram and the combination of the processes or / and blocks in the flowchart or / and block diagram can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart or / and block diagram. Figure 1 Process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0151] These computer program instructions may also be stored in a computer readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture including an instruction device, which implements the process Figure 1 Process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

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

[0153] The unit in this application is an object that objectively describes the morphological structure with the help of physical or virtual expressions. The object is not equal to the object and is not limited to physical and virtual. It can be a data processing function, software program, processing mode, usage method, operation method, workflow, application process, electronic hardware, circuit module, processing system, system imitation or simulation object.

[0154] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, ordinary technicians in the relevant field can still modify or replace the specific implementation methods of the present invention with equivalents. Any modification or equivalent replacement that does not depart from the spirit and scope of the present invention should be included in the scope of protection 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 and the unsaturated and unsteady seepage algorithm, the geometric grid model is processed to obtain the steady-state seepage field and extract the pore pressure results; Step 3: Using the pre-built arc sliding surface calculation unit, determine the most dangerous arc sliding surface according to the pore pressure results; Step 4: using the pre-built slope sliding surface calculation unit, based on the most dangerous arc sliding surface, solve the yield acceleration coefficients of several slope sliding surfaces, and determine the most dangerous slope sliding surface according to the yield acceleration coefficients; 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: Step 3: Using the pre-built arc sliding surface calculation unit, the method for determining the most dangerous arc sliding surface according to the pore pressure results 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, so as 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 array of arc sliding surfaces, and search near the minimum value of the yield acceleration coefficient through local dichotomy to find the most dangerous arc sliding surface.

5. The grid-based multi-layer slope yield acceleration calculation method according to claim 4, characterized in that: The first step is to calculate the radius search upper and lower limits of the potential arc sliding surface according to the geometric mesh model. The method to obtain the upper and lower limits of the arc radius 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; The upper and lower limits of the arc radius are summed up to obtain the upper and lower limits of the arc radius.

6. The grid-based multi-layer slope yield acceleration calculation method according to claim 5, 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.

7. 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.

8. 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.

9. A grid-based method for calculating yield acceleration of multi-layer slopes, characterized in that: Includes the following: According to the geological data of a multi-layer slope, the geometric grid model is generated by coupling the hydraulic parameters and mechanical shear information; The geometrical grid model is processed to obtain the steady-state seepage field; According to the steady-state seepage field, the most dangerous arc sliding surface is determined; Based on the most dangerous arc sliding surface, the most dangerous slope sliding surface is obtained; The most dangerous slope sliding surface is processed to obtain the yield acceleration data of a multi-layer slope.

10. 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-9.

Images