Method, medium and equipment for evaluating stability evolution of unsaturated soil slope in cold region

Abstract

The invention provides a method, medium and equipment for evaluating stability evolution of an unsaturated soil slope in a cold region, and establishes a framework combining a soil-atmosphere interaction model and a cold region unsaturated soil water-steam-heat coupling transmission model. A finite element transient simulation program for simulating the water-heat state response of the unsaturated soil body in the cold region under the action of long-term and short-term climates is developed, and the technical problem that when a traditional numerical model simulates the sand freezing and thawing phenomenon, the latent heat release amount is prone to losing due to the narrow phase change temperature interval is solved. According to a limit analysis upper limit theorem, a slope stability evaluation method is established by adopting a discretization recursion technology, and a novel iterative solution scheme combining a dichotomy and a meta-heuristic optimization algorithm is innovatively provided. A finite element transient analysis model-stability evaluation coupling method is provided, so that the stability evolution evaluation of the unsaturated soil slope in the cold region driven by climate data can be carried out.

Description

Technical Field

[0001] The present invention relates to the field of environmental monitoring, and in particular to a method, a medium and a device for evaluating the stability evolution of an unsaturated soil slope in a cold region. Background Art

[0002] In recent years, with the frequent occurrence of abnormal weather around the world and the surge in human social and economic activities, a large number of landslide disasters and slope damage events have occurred in mountains and along transportation facilities in cold regions. The evolution process of geological disasters (gestation, formation, movement and disaster-causing) is closely related to the coupling interaction of various spheres of the earth system. The near-surface soil gradually loses its strength under the influence of seasonal climate and environmental conditions, inducing a series of natural disasters such as landslides and collapses.

[0003] However, there are few studies that directly link dynamically changing climate with slope stability assessment. Most studies on the impact of climate on slope stability are based on narrative development and empirical evaluation. Current analysis methods or related computational software are generally based on snapshots of the hydrological state of the field soil investigated during a specific period of time, and cannot reflect the evolution of slope stability under the influence of the dynamic climate system. Extending to engineering practice, design and assessment are almost always based on hydraulic profiles drawn during a short-term survey. However, due to seasonal climate change, when the actual hydraulic state of the slope deviates from the threshold hydraulic state considered in the initial assessment stage, the assessment of stability based on traditional methods will not match the actual state. In addition, most stability analyses are only conducted on slopes in areas without freeze-thaw activity, and studies on slopes in cold areas mostly use qualitative analysis, and relevant computational theories and methods are particularly scarce. The few quantitative assessments almost all use infinite slope models and assume the distribution of underground ice and hydrology.

[0004] The above factors seriously restrict the assessment of long-term stability of slopes in cold areas under the influence of climate. Summary of the invention

[0005] The purpose of the present invention is to provide a method, medium and equipment for evaluating the stability evolution of unsaturated soil slopes in cold regions, so as to improve the above-mentioned problems.

[0006] In order to achieve the above purpose, the technical solution adopted by the embodiment of the present invention is as follows:

[0007] In a first aspect, an embodiment of the present invention provides a method for evaluating the stability evolution of an unsaturated soil slope in a cold region, the method comprising:

[0008] S1, finite element modeling of unsaturated soil slope, including: generating a slope geometric model according to geometric control points and stratum control lines, dividing the slope geometric model into finite element meshes using a mesh generator, and assigning basic physical and mechanical parameters of the finite element meshes and nodes;

[0009] S2, fitting the sub-daily weather data, including: extracting daily data of target meteorological variables from the meteorological data set, wherein the target meteorological variables include any one or more of air temperature, relative humidity, wind speed, sunshine hours and air pressure, and using a trigonometric function type weather data fitting program to generate hourly sequence weather data;

[0010] S3, determining the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface, including: constructing a balance equation and a parameterized model based on the energy balance and mass balance of the soil-atmosphere interface, inputting the hourly series weather data obtained by S2 into the parameterized model, and solving the balance equation to obtain the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface;

[0011] S4, water-vapor-heat multi-field coupled finite element simulation of unsaturated soil slope in cold region, including: constructing a water-vapor-heat multi-field coupled finite element simulation model of unsaturated soil slope in cold region, taking the heat and moisture flux exchanged at the soil-atmosphere interaction interface obtained from S3 as the upper boundary conditions, driving the water-vapor-heat multi-field coupled finite element simulation model of unsaturated soil slope in cold region to run, and obtaining the continuous changes of the internal hydrological and thermal states of unsaturated soil slope in cold region under the forcing of meteorological data after solving;

[0012] S5, obtain the evolution information of the stability of unsaturated soil slopes in cold regions, including: S501 and S502;

[0013] S501, establish a discretized failure mechanism generator based on the upper limit method of limit analysis, including: using the finite element grid and nodes established in S1 as the background, constructing the entire soil slope area as a Delaunay triangle diagram, using the discretized point-to-point recursive technology, generating potential failure mechanisms in the Delaunay triangle area, interpolating the hydrological and thermal state variable values ​​of the nodes obtained in step S4 to determine the hydrological and thermal state variable values ​​at the discrete points, and calculating the internal energy dissipation rate and external force power based on this, and then solving and outputting the sliding surface information and the energy ratio R corresponding to the potential sliding surface a , where the energy ratio R a is the ratio of internal energy dissipation rate to external force power;

[0014] S502, construct an iterative solution for slope stability, use the strength reduction method to define the safety factor, and build a new iterative solution that combines the dichotomy method and the meta-heuristic optimization algorithm. At the beginning of each iteration, the dichotomy method is used to update the safety factor. At this time, the soil strength parameter is the soil strength parameter after the safety factor is reduced. The energy ratio R in S501 is used. a As the objective function, a meta-heuristic optimization algorithm is used to solve the optimal sliding surface and the corresponding energy ratio R under a given safety factor. a , until the safety factor determined by the dichotomy method reaches the minimum value, and the energy ratio Ra When the difference between the value and 1 meets the threshold error, the iterative calculation is terminated and the safety factor and critical slip surface information are output;

[0015] S503: Repeat the method described in S502 to calculate the slope safety factor and critical sliding surface information over time, thereby obtaining the evolution information of the slope stability.

[0016] In a second aspect, an embodiment of the present invention provides a storage medium having a computer program stored thereon, and the computer program implements the above method when executed by a processor.

[0017] In a third aspect, an embodiment of the present invention provides an electronic device, comprising: a processor and a memory, wherein the memory is used to store one or more programs; when the one or more programs are executed by the processor, the above method is implemented.

[0018] Compared with the prior art, the embodiment of the present invention provides a method, medium and device for evaluating the stability evolution of unsaturated soil slopes in cold regions, establishes a framework combining the soil-atmosphere interaction model and the water-vapor-heat coupling transmission model of unsaturated soil in cold regions, and develops a finite element transient simulation program that simulates the water-heat state response of unsaturated soil in cold regions under long- and short-term climate effects, solving the technical problem that the traditional numerical model is prone to lose the amount of latent heat release due to the narrow phase change temperature range when simulating the freeze-thaw phenomenon of sandy soil. According to the upper limit theorem of limit analysis, a slope stability evaluation method is established using discretization recursive technology, and an innovative iterative solution combining the dichotomy method and the meta-heuristic optimization algorithm is proposed. A method of coupling the finite element transient analysis model and stability evaluation is proposed, which can be used to evaluate the stability evolution of unsaturated soil slopes in cold regions driven by climate data.

[0019] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments will be briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without paying creative work.

[0021] Figure 1 A schematic flow chart of a method for evaluating the stability evolution of unsaturated soil slopes in cold regions provided by an embodiment of the present invention;

[0022] Figure 2A schematic diagram of the discretization failure mechanism of an unsaturated soil slope provided by an embodiment of the present invention;

[0023] Figure 3 The geometric model, finite element mesh and boundary conditions provided for the embodiments of the present invention;

[0024] Figure 4 A fitted sub-daily weather data graph provided by an embodiment of the present invention;

[0025] FIG5( a) is a diagram showing changes in energy flux on a slope surface provided by an embodiment of the present invention;

[0026] FIG5( b ) is a diagram showing changes in evaporation on a slope surface provided by an embodiment of the present invention;

[0027] FIG6( a ) is a diagram showing changes in temperature, unfrozen water (liquid water) content and ice content of a monitoring section AA′ provided by an embodiment of the present invention;

[0028] FIG6( b ) is a diagram showing changes in temperature, unfrozen water (liquid water) content and ice content of the monitoring section BB′ provided in an embodiment of the present invention;

[0029] FIG6( c ) is a diagram showing changes in temperature, unfrozen water (liquid water) content and ice content of the monitoring section CC′ provided by an embodiment of the present invention;

[0030] FIG. 7( a ) is a calculation result of the safety factor FoS provided by an embodiment of the present invention;

[0031] FIG. 7( b ) is a partial enlarged view of the melting period safety factor FoS provided by an embodiment of the present invention;

[0032] Figure 8 A schematic diagram of the structure of an electronic device provided by an embodiment of the present invention.

[0033] In the figure: 10 - processor; 11 - memory; 12 - bus; 13 - communication interface. DETAILED DESCRIPTION

[0034] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

[0035] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0036] It should be noted that similar reference numerals and letters represent similar items in the following drawings, so once an item is defined in one drawing, it does not need to be further defined and explained in the subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are only used to distinguish the description and cannot be understood as indicating or implying relative importance.

[0037] It should be noted that, in this article, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "comprise a ..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

[0038] In the description of the present invention, it should be noted that the terms "upper", "lower", "inside", "outside", etc. indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings, or are directions or positional relationships in which the product of the invention is usually placed when in use. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore should not be understood as a limitation on the present invention.

[0039] In the description of the present invention, it is also necessary to explain that, unless otherwise clearly specified and limited, the terms "disposed" and "connected" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be an indirect connection through an intermediate medium, or it can be the internal communication of two elements. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0040] Some embodiments of the present invention are described in detail below in conjunction with the accompanying drawings. In the absence of conflict, the following embodiments and features in the embodiments can be combined with each other.

[0041] In the case of the above-mentioned problems, it is necessary to develop a model based on physical foundations to promote the stability assessment of mountain slopes and engineering facility slopes in cold regions under long-term climate effects, especially to develop physical models that can handle dynamic processes and threshold crossings. In view of the lack of methods in this regard, an embodiment of the present invention provides a method for evaluating the evolution of stability of unsaturated soil slopes in cold regions, so as to evaluate the evolution of stability of unsaturated soil slopes in cold regions over time under the influence of climate, and innovatively couples the climate data-driven water-vapor-heat coupling transmission finite element analysis model of unsaturated soil slopes in cold regions with the slope stability calculation model, providing a guiding modeling method for the deductive process of evaluating the long-term stability of near-surface unsaturated soil slopes in cold regions under the action of multiple meteorological factors, which can be used to evaluate the long-term service performance and potential damage risks of unsaturated slopes in cold regions, and provide scientific advice for the design, construction and long-term operation and maintenance of slopes in cold regions.

[0042] The method for evaluating the stability evolution of unsaturated soil slopes in cold regions provided by the embodiment of the present invention can be applied to, but not limited to, the electronic devices described below. For the specific process, please refer to Figure 1 , the methods for evaluating the stability evolution of unsaturated soil slopes in cold regions include:

[0043] S1, finite element modeling of unsaturated soil slope, including: generating slope geometric model according to geometric control points and stratum control lines, dividing the slope geometric model into finite element meshes using a mesh generator, and assigning basic physical and mechanical parameters of finite element meshes and nodes.

[0044] Among them, the nodes belong to the finite element grid. Considering that the surface soil is strongly affected by climate, the grid corresponding to the surface soil is appropriately encrypted, and a basic grid information database is established to input and store the basic physical and mechanical parameters of each finite element grid and node.

[0045] S2, fitting of sub-daily weather data, including: extracting daily data of target meteorological variables from the meteorological data set, the target meteorological variables including any one or more of air temperature, relative humidity, wind speed, sunshine hours and air pressure, and using a trigonometric function type weather data fitting program to generate hourly series weather data.

[0046] Since the daily data of the target meteorological variables obtained are all daily values, they need to be converted into sub-daily scale data. The sub-daily weather data fitting program (trigonometric function weather data fitting program) is used to fit the daily meteorological data into hourly series weather data. The converted weather data and the longitude, latitude and altitude information of the slope location are saved as input files.

[0047] S3, determining the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface, including: constructing a balance equation and a parameterized model based on the energy balance and mass balance of the soil-atmosphere interface, inputting the hourly series weather data obtained by S2 into the parameterized model, and solving the balance equation to obtain the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface.

[0048] S4, water-vapor-heat multi-field coupled finite element simulation of unsaturated soil slopes in cold regions, including: constructing a water-vapor-heat multi-field coupled finite element simulation model for unsaturated soil slopes in cold regions, using the heat and moisture flux exchanged at the soil-atmosphere interaction interface obtained from S3 as the upper boundary conditions, driving the water-vapor-heat multi-field coupled finite element simulation model for unsaturated soil slopes in cold regions, and obtaining the continuous changes in the internal hydrological and thermal states of unsaturated soil slopes in cold regions under the forcing of meteorological data after solving the problem.

[0049] Optionally, the continuous changes in the internal hydrological and thermal states of the unsaturated soil slope in the cold region under the forcing of meteorological data include the distribution of soil temperature, water content, ice content and pore fluid pressure in the cold region over time.

[0050] S5, obtain the evolution information of the stability of unsaturated soil slopes in cold regions, including: S501 and S502;

[0051] S501, establish a discretized failure mechanism generator based on the upper limit method of limit analysis, including: using the finite element grid and nodes established in S1 as the background, constructing the entire soil slope area as a Delaunay triangle diagram, using the discretized point-to-point recursive technology, generating potential failure mechanisms in the Delaunay triangle area, interpolating the hydrological and thermal state variable values ​​of the nodes obtained in step S4 to determine the hydrological and thermal state variable values ​​at the discrete points, and calculating the internal energy dissipation rate and external force power based on this, and then solving and outputting the sliding surface information and the energy ratio R corresponding to the potential sliding surface a , where the energy ratio R a is the ratio of internal energy dissipation rate to external force power;

[0052] S502, construct an iterative solution for slope stability, use the strength reduction method to define the safety factor, and build a new iterative solution that combines the dichotomy method and the meta-heuristic optimization algorithm. At the beginning of each iteration, the dichotomy method is used to update the safety factor. At this time, the soil strength parameter is the soil strength parameter after the safety factor is reduced. The energy ratio R in S501 is used. a As the objective function, a meta-heuristic optimization algorithm is used to solve the optimal sliding surface and the corresponding energy ratio R under a given safety factor. a , until the safety factor determined by the dichotomy method reaches the minimum value, and the energy ratio R a When the difference between the value and 1 meets the threshold error, the iterative calculation is terminated and the safety factor and critical slip surface information are output;

[0053] S503: Repeat the S502 method to calculate the slope safety factor and critical sliding surface information over time, thereby obtaining the evolution information of the slope stability.

[0054] The method for evaluating the stability evolution of unsaturated soil slopes in cold regions provided by the embodiment of the present invention establishes a framework combining the soil-atmosphere interaction model and the water-vapor-heat coupling transmission model of unsaturated soil in cold regions, and develops a finite element transient simulation program that simulates the water-heat state response of unsaturated soil in cold regions under long- and short-term climate effects, solving the technical problem that the latent heat release amount is easily lost due to the narrow phase change temperature range when simulating the freeze-thaw phenomenon of sandy soil with traditional numerical models. According to the upper limit theorem of limit analysis, a soil slope stability evaluation method is established using discretization recursive technology, and an innovative iterative solution combining the dichotomy method and the meta-heuristic optimization algorithm is proposed. A method of coupling the finite element transient analysis model and the stability evaluation is proposed, which can be used to evaluate the stability evolution of unsaturated soil slopes in cold regions driven by climate data.

[0055] Air temperature, relative humidity and wind speed generally show trigonometric fluctuations on a sub-daily time scale, so trigonometric functions can be used to fit their sub-daily variations. Optionally, S2, the step of fitting sub-daily weather data, includes: fitting the daily weather data of the target meteorological variable obtained from the meteorological data set to the sub-daily fluctuations of the weather data through equations (1) to (3), wherein the sub-daily fluctuations are hourly series of weather data:

[0056]

[0057] Where t represents the current time, T a is the temperature corresponding to the current time t, T m is the daily average temperature of the day at the current time t, T max is the highest temperature of the day at the current time t, T min is the lowest temperature of the day at the current time t, RH a is the relative humidity of the air corresponding to the current time t, RH a,m is the daily average relative humidity of the air at the current time t, RH a,max is the maximum relative humidity of the day at the current time t, V a is the wind speed corresponding to the current time t, V m The daily average wind speed for the day of the current time t, V max The maximum wind speed of the day at the current time t, t max,T is the time when the daily maximum temperature occurs, t max,RHa is the time when the maximum relative humidity of the air occurs each day, t max,VThe time when the maximum wind speed occurs each day.

[0058] Optionally, S3, determining the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface, and calculating the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface by equations (4) and (5), respectively:

[0059] G=R n -LE-H (4)

[0060] q net =PRE / ρ l (5)

[0061] Where G is the heat exchanged at the soil-atmosphere interaction interface, R n is the net radiation flux, LE is the latent heat flux, L is the latent heat of phase change from liquid water to gaseous water, E is the evaporation, H is the sensible heat flux, q net is the moisture flux exchanged at the soil-atmosphere interaction interface, P is the precipitation, R is the runoff, and ρ l is the density of liquid water.

[0062] Net radiation flux is obtained by adding short-wave and long-wave radiation. Short-wave radiation can be calculated by theoretical formulas based on local latitude, longitude, calculation time, surface soil moisture, sunshine duration, etc. Long-wave radiation is calculated by Stefan-Boltzmann radiation formula (requires air temperature and surface soil temperature) based on air emissivity (calculated by relative air humidity and air temperature), surface emissivity (surface soil moisture), etc. Other energy fluxes such as LE and H are determined by Monin-Obukhov similarity theory.

[0063] Optionally, S4, water-vapor-heat multi-field coupled finite element simulation of unsaturated soil slope in cold regions, comprises the following steps:

[0064] S401, establish the control equation of water-vapor-heat multi-field coupling of unsaturated soil: Taking unit volume soil as the research object, based on the principle of water mass conservation and energy conservation, establish the strong form control equation of the coupling system:

[0065]

[0066] Among them, ρ l is the density of liquid water, ρ i is the density of ice, ρ v is the density of gaseous water, θ l is the volume content of liquid water, θ i is the volume content of ice, θ a is the volume content of air, θ v is the volume content of water vapor, is the gradient operator, Klh is the isothermal liquid water migration coefficient, K lT is the non-isothermal liquid water migration coefficient, K vh is the isothermal water vapor diffusion coefficient, K vT is the non-isothermal water vapor diffusion coefficient, h is the suction head, y is the vertical spatial coordinate, q l is the amount of liquid water transferred, q v is the amount of water vapor transferred, C is the volume heat capacity of the soil, C l is the volumetric heat capacity of liquid water, C v is the volume heat capacity of water vapor, L0 is the volume latent heat of vaporization of liquid water, L f is the latent heat of melting of ice, T is the temperature, and λ is the effective thermal conductivity of the soil;

[0067] It should be noted that, according to the principle of conservation of mass, the total mass conservation equation of water is established as shown in formula (6), including the mass change rate of liquid water, ice and water vapor, as well as the isothermal and non-isothermal liquid water migration mass and the isothermal and non-isothermal water vapor migration mass per unit time. According to the principle of conservation of energy, the energy conservation equation is established as shown in formula (7), including the internal energy change rate, the water-ice phase change latent heat release rate, the water-gaseous water phase change latent heat release rate, the conductive heat flux, the convective heat flux carried by the liquid water and water vapor migration, and the vaporization latent heat flux carried by the water vapor migration.

[0068] The number of unknown variables contained in equations (6) and (7) exceeds the number of equations, and other relationships need to be introduced to close the equations, so S402 needs to be executed.

[0069] S402, Degenerate treatment of the strong form governing equations of coupled systems:

[0070] The van-Genuchten soil-water characteristic curve (SWCC, describing the constitutive relationship between soil matrix suction head and water content - soil-water characteristic curve) is expressed as formula (8), and its derivative relationship is formula (9):

[0071]

[0072] Among them, θ s is the volumetric water content at saturation, θ r is the volumetric water content in the residual state, α v 、n v and m v is the fitting parameter, S e is the effective saturation, and the derivative value C' is the specific water capacity;

[0073] Secondly, once the soil temperature is lower than the freezing temperature, the driving force for water migration will be provided by the low-temperature suction head. Therefore, the relationship between temperature and low-temperature suction head needs to be introduced:

[0074]

[0075] Wherein, h0 is the pressure head before the soil is frozen (specifically the matrix suction head before the unsaturated soil is frozen); HS is the Heaviside function; ω is the ratio of the water-ice interface free energy to the water-air interface free energy, which is 1 for colloidal soil particles and 2.2 for non-colloidal soil particles; T' is the freezing temperature;

[0076] Combining equations (8) and (10), we can establish the relationship between unfrozen liquid water and negative temperature (11), which is called the soil freezing characteristic curve (SFCC):

[0077]

[0078] Let θ = θ l +ρ i θ i / ρ l , then the pore air content is expressed as the difference between the porosity n and θ: θ a =n-θ, substituting the above relationship into equation (6), we get the degenerated total water mass conservation equation:

[0079]

[0080] in, K HH =ρ l (K lh +K vh ), K HT =ρ l (K lT +K vT );

[0081] Based on equation (7), the same order and the same operator terms are merged to obtain the simplified form of the energy conservation equation:

[0082]

[0083] in, L ii =L f ρ i , K TT =(λ+L0K vT ), K TH =L0K vh and

[0084] The coupled system composed of equations (12) and (13) can be described as follows:

[0085]

[0086]

[0087] Formula (14) represents the equation in a D-dimensional Euclidean space R with a closed boundary Γ. D In the physical domain Ω, at time t greater than or equal to t0, the state variables h, T and θ i Satisfied by The partial differential equation system composed of , in addition, on the boundary Γ, satisfies the following conditions:

[0088] Dirichlet frontier: Newman Boundary:

[0089] And at time t0, the initial variable value is:

[0090]

[0091] Where n is the unit normal vector outside the boundary, q and Q are the water flow rate and heat flow rate on the boundary, q w is the water flux exchanged at the soil-atmosphere interaction interface, q h heat exchanged at the soil-atmosphere interaction interface;

[0092] S403, finite element space discretization is performed on equation (14), and the Galerkin weighted residual method is used to integrate equation (14), and the equivalent integral weak form of the partial differential equation of total water mass and energy conservation can be obtained:

[0093]

[0094] Among them, L w is the residual of the equivalent integral equation of the total water mass, L h is the residual of the energy equivalent integral equation, N is the shape function, N T is the transpose of the shape function;

[0095] For equation (17), firstly solve it in unit integration to get the unit stiffness equation, and then assemble it to get the spatial discretization of the moisture equation:

[0096] form:

[0097]

[0098] Among them, the expressions of each coefficient matrix are as follows:

[0099]

[0100]

[0101] Among them, nel is the total number of units, n b is the number of nodes on the unit boundary to which the Newman boundary condition belongs;

[0102] Similarly, for equation (18), firstly solve it in the unit integral to obtain the unit stiffness equation, and then assemble it to obtain the space equation of the energy equation:

[0103] Discretized form:

[0104]

[0105] Among them, the expressions of each coefficient matrix are as follows:

[0106]

[0107]

[0108]

[0109]

[0110] Equations (19) and (21) are processed by time discretization to obtain their conventional iterative formats. In order to overcome the problem that the conventional format is prone to lose latent heat information and it is difficult to ensure the conservation of total water mass due to its narrow temperature range when simulating the freezing and thawing phenomenon of sandy soil, an iterative level operator splitting method is introduced to transform the time term. Specifically, assuming that t is the previous time step, t+Δt is the current time step, the current iteration level is marked as j+1, the first iteration is 1, and the value of the first iteration level is inherited from the convergence value at time t. In the iterative calculation of [t, t+Δt], the time term is split into two terms, 1~j terms and j~j+1 terms:

[0111]

[0112] First, for the second term on the far right of equation (23), restore its original meaning, that is, the rate of change of the mass of liquid water and ice:

[0113]

[0114] Among them, M ww and M ii The expression is as follows:

[0115]

[0116]

[0117] Next, keep the second term on the right of equation (24) and equation (25) unchanged, and add the first term to obtain:

[0118]

[0119] Among them, M Ti The expression is:

[0120]

[0121] Coefficient M Ti is the sensible heat capacity, and its expression is:

[0122]

[0123] Substituting equations (23) and (26) into equation (19), and substituting equations (24), (25) and (28) into equation (21), combined with the backward Euler difference method, the final iterative format of the discretized equation is obtained:

[0124]

[0125]

[0126] Among them, the continuous changes of the hydrological and thermal states of the unsaturated soil slopes in the cold region under the forcing of meteorological data include T t +Δt,j+1 and h t+Δt,j+1 ; T t+Δt,j represents the temperature value of the ramp node obtained after the jth iteration at the current time t+Δt, T t +Δt,j+1 represents the temperature value of the slope node obtained after the j+1th iteration at the current time t+Δt, h t+Δt,j represents the water head value of the slope node obtained after the jth iteration at the current time t+Δt, h t+Δt,j+1 represents the water head value of the slope node obtained after the j+1th iteration at the current time t+Δt, T t represents the temperature value of the ramp node after the previous time step t converges, θ t represents the total volume content of liquid water and ice at the slope node after convergence at the previous time step t, represents the ice volume content of the slope node obtained after the jth iteration at the current time t+Δt, represents the volume content of ice at the slope node after convergence at the previous time step t;

[0127] S404, compile a finite element solution program, propose and adopt a single time step sequential solution scheme: that is, divide the current time step into two sub-steps, in the first sub-step, assume that the temperature is constant to solve the moisture discretization equation, update the pressure head and total water content; in the second sub-step, assume that the pressure head and total water content are constant to solve the energy discretization equation, update the temperature and ice content; after the two sub-steps are solved sequentially, iterate and loop the calculation until the convergence condition is met, output and store the results, and the results include the T t+Δt,j+1 、h t+Δt,j+1 , as well as in, represents the ice volume content of the slope node at the current time t+Δt, after the j+1th iteration, represents the volume content of liquid water at the slope node obtained after the j+1th iteration at the current time t+Δt, By changing T t+Δt,j+1 Substitute into formula (11) to calculate and determine: By changing h t+Δt,j+1 After conversion to water content Determine after subtraction.

[0128] Optionally, S501, a discretized failure mechanism generator based on the limit analysis upper limit method is established, and its premise is that under given strength parameters, its establishment includes the following three steps:

[0129] (1) The Delaunay triangle diagram is constructed with the finite element mesh and nodes as the background. The node variable values ​​and soil strength parameters calculated in step S4 are assigned to the vertices of the Delaunay triangle. The slope toe point is taken as the starting point and the rotation radius r0 and rotation angle θ0 are taken as the starting conditions. A series of discrete points are generated in sequence by using the recursive method and following the ideal plastic orthogonal condition. The ideal plastic orthogonal condition indicates that the angle between the velocity direction of the discrete point and the sliding surface is the friction angle at the discrete point. The node variable values ​​include the continuous changes of the internal hydrological and thermal states of the unsaturated soil slope in the cold region under the forcing of meteorological data. The soil strength parameters include the internal friction angle and cohesion. It is known that the previous discrete point P i The corresponding rotation angle β i , Radius of rotation r i , then the recursively generated point P i+1 The rotation angle β i+1 , Radius of rotation r i+1 and coordinates (x i+1 ,y i+1 ) are calculated using formulas (33) to (35) respectively:

[0130] β i+1 =β i +δ (33)

[0131]

[0132] Among them, (x A ,y A ) is the coordinate of the starting point A of the slope toe, which is determined by the slope geometry and is known; δ is the angle between adjacent rotation radii, which is specified by the user; is a discrete point P i The internal friction angle at P is obtained by first finding i The position of the point in the Delaunay triangle diagram returns the number of the triangle and the vertex (node) number, and then uses the barycentric interpolation method to determine P i Angle of internal friction at point:

[0133]

[0134] Among them, κ, ξ and ζ are interpolation coefficients that vary with the position of the discrete point i in the triangle. and are the values ​​of the internal friction angles at the triangle vertices a, b, and c respectively;

[0135] When the newly generated discrete point P i+1 When the slope line is exceeded, the recursive procedure will terminate and the last discrete point will be corrected to the slope line. Then, all discrete points are connected to form a potential slip surface. The slip lines between adjacent discrete points are shown in P i-1 P i The length L i Determined as:

[0136]

[0137] (2) It is also necessary to calculate the internal energy dissipation rate occurring on the slip surface and the power of all external forces. The internal energy dissipation rate only occurs on the slip surface. The total dissipation rate is calculated by summing the micro-segment dissipation rates:

[0138]

[0139] Where c is the cohesion, v is the velocity of the sliding surface, n is the total number of discretized segments of the sliding surface, dL is the length of the micro-segment of the sliding surface, and c i is the cohesive force at the discrete point. The pore fluid pressure also does work on the shear zone of the slip surface. The power of the pore fluid pressure is calculated by the micro-segment power summation method using formula (39):

[0140]

[0141] Among them, W u is the pore pressure power, u is the pore fluid pressure, u i is the pore fluid pressure at discrete point i on the slip surface, v iis the movement velocity at discrete point i on the slip surface;

[0142] The soil in different areas of the sliding body will have different densities and weights due to differences in components. Therefore, the closed area surrounded by two adjacent rotation radii, slope lines and sliding lines is divided into many small quadrilaterals. Each quadrilateral is further divided into a series of smaller paired triangles. The weight γ of the soil at the centroid of each small triangle is determined by the barycentric interpolation method:

[0143] γ=ρ l gθ l +ρ s g(1-n)+ρ a gθ a +ρ i gθ i (40)

[0144] Where g is the acceleration due to gravity;

[0145] After calculating the gravity power of all small triangles in the sliding area, the total gravity power is calculated by summing the micro-units:

[0146]

[0147] Among them, W g is the total gravitational power, N is the total number of small triangles, γ i is the weight at the centroid of the i-th small triangle, S i is the area of ​​the ith small triangle, v i is the speed of the ith small triangle (the direction is perpendicular to the line connecting the centroid and the rotation center), β i ′ is the angle between the line connecting the centroid and the rotation center and the x-axis;

[0148] (3) Define the energy ratio as:

[0149]

[0150] For a given strength parameter, when the calculated energy ratio reaches the minimum value, the corresponding slip surface is the optimal slip surface under the current conditions.

[0151] Optionally, in S502, a novel iterative solution combining the dichotomy method and the meta-heuristic optimization algorithm is used to solve the safety factor and the optimal sliding surface of the slope, and the implementation steps are: the safety factor (FoS) is defined by the strength reduction method, and after each iteration, a new safety factor is obtained by the dichotomy method, and then the original strength parameter is reduced by the safety factor to obtain a set of new cohesion and internal friction angle, and the new cohesion, internal friction angle and other state variable values ​​are input into the energy ratio R a It is the optimization procedure of the objective function. After iterative calculation, the energy ratio Ra When it is the minimum value under the given strength parameter, the iteration ends, and a smaller safety factor is calculated by dichotomy, and the next iteration starts until the minimum value of the safety factor is found. The safety factor value at this time and the critical slip surface information under the safety factor value are output. Among them, other state variable values ​​include temperature, fluid pressure, water content and ice content.

[0152] Optionally, a meta-heuristic optimization algorithm is used to determine the minimum energy ratio R under a given safety factor. a and the corresponding slip surface. The variables of the optimization program are the starting rotation radius r0 and the starting rotation angle θ0 with specified upper and lower limits. The variable seeds are randomly arranged according to the optimization algorithm. Under each group of variable combinations, the S501 program is called to generate potential failure mechanisms, thereby generating different potential slip surfaces. The optimal slip surface under the current safety factor is obtained through optimization.

[0153] Optionally, in a novel iterative solution scheme combining the bisection method and the meta-heuristic optimization algorithm, the reduced strength parameters include the cohesion c' and the safety factor Calculate using formula (43):

[0154]

[0155] The bisection method shown in formula (44) is used to update the safety factor after each iteration:

[0156]

[0157] Among them, FoS ub and FoS lb It is a repository of safety factors. When the iterative FoS is greater than or equal to 1, it is stored in the database FoS ub When the iteratively obtained FoS is less than 1, it is stored in the database FoS lb middle.

[0158] The embodiment of the present invention also provides an optional implementation method, please refer to the following.

[0159] Step S1: Establishing a finite element model of an unsaturated soil slope;

[0160] The embodiment is a high-speed railway embankment slope in a certain area, and the slope ratio is 1:2. According to step S1, a geometric model and a finite element mesh of the slope are established, such as Figure 3 The figure shows the boundary conditions of the finite element model. The surface mesh is refined to accurately capture the exchange flux between the soil and the atmosphere. The top of the model is the atmospheric boundary. The bottom boundary is a fixed water head and is set to 0.03W / m 2The constant geothermal flux is 1.3°, and the side boundaries are insulated and impermeable. Monitoring profiles were set up at the toe, middle and top of the slope to record the continuous changes in temperature, water content and ice content during the simulation period. Table 1 shows the hydraulic parameters and structural conditions of the sandy silt that constitutes the slope, and Table 2 shows the physical and mechanical parameters of the sandy silt.

[0161] Table 1 Hydraulic parameters and structural conditions of sandy silt

[0162]

[0163] Table 2 Physical and mechanical parameters of sandy silt

[0164]

[0165] Step S2: fitting of sub-daily weather data;

[0166] The daily weather data of a certain area’s weather station from July 2015 to July 2019 were extracted from the meteorological data daily value dataset, such as Figure 4 shown. Figure 4 The sub-daily weather changes from May 2 to May 17, 2017, fitted by step S2, are used. The weather data from July 1, 2015 to June 30, 2016 are used as the forcing data to obtain the initial conditions, and the meteorological data from July 2016 to October 2019 are used for the formal simulation.

[0167] Step S3: determining the exchange flux at the soil-atmosphere interaction interface;

[0168] The basic model information and meteorological data obtained in step S1 and step S2 are loaded into the finite element transient analysis model in the form of a file to drive the finite element program to run. Step S3 will be executed first to calculate the heat and moisture flux exchanged between the soil and the atmosphere. Figure 5(a) shows the daily variation sequence of various energy fluxes on the slope surface. Considering that the local rainfall is not significant and the runoff is small, Figure 5(b) only shows the daily variation sequence of evaporation.

[0169] Step S4: calculation of soil water-vapor-heat coupling model;

[0170] Then, in the finite element transient model, the flux calculated by step S3 is used as the upper boundary condition, and the driver continues to run. According to the two-step solution described in S4, the hydrological state change of the underlying unsaturated soil is first calculated, and then the thermal state change is calculated. The cycle is iterated until the current time step converges. The convergence value is stored, the time is updated and the calculation of the next time step is started. The process is repeated until the end of the time period is reached and the calculation is terminated. Please refer to Figure 6(a), Figure 6(b) and Figure 6(c). Figure 6 is the temperature cloud map of the three monitoring sections extracted after the finite element model calculation is completed, as well as the changes in water content and ice content at different depths of the three monitoring sections.

[0171] Step S5: Calculate the evolution of slope stability over time;

[0172] After obtaining the simulation results, the hydrological and thermal states at a certain time are extracted according to the specified time step to construct the potential failure mechanism of the slope. As mentioned above, the potential failure mechanism relies on the recursive algorithm described in step S501. After all the information of the previous discrete point is known, the position of the next discrete point needs to meet the orthogonality mandatory condition specified by the ideal plastic material, that is, Figure 2 The angle between the vertical direction (movement speed) of the rotation radius shown and the micro-segment line of the sliding surface must be the friction angle at the discrete point. For non-frozen soil, the friction angle remains unchanged; for frozen soil, the frozen soil strength model proposed by Nishiruma and Wang (2019) is adopted, and its friction angle and cohesion are related to the low-temperature suction, which can be determined by equations (45) to (48):

[0173] M f =M a +a M s c (45)

[0174] q f =a q s c (46)

[0175]

[0176] Among them, q f is the shear strength enhancement caused by pore ice, a q 、a M is a constant in the model, s c (MPa) is the low temperature suction, M f is the slope of the critical state line of frozen soil, M ais the slope of the critical state line of the melting soil. In step S501, when the soil at the discrete point is in a frozen state, the friction angle determined by formula (47) is used to recursively obtain the next discrete point. After the potential failure mechanism of the slope is constructed, the internal energy dissipation rate and the external force power are calculated, and then the energy ratio R of the corresponding failure mechanism is calculated. a .

[0177] Figure 2 middle, is the inclination angle of the slope surface at the toe, α is the inclination angle of the slope surface at the top of the slope, P i ' is the radius of rotation r i The intersection point with the slope line, P′ i+1 is the radius of rotation r i+1 The intersection of the slope line, H is the slope height, B is the slope vertex, and area P i-1 P′ i-1 P i ′P i is a small quadrilateral, triangle D i,j D i+1,j D i+1,j+1 and D i,j D i,j+1 D i+1,j+1 It is one of a series of paired triangles formed after splitting the small quadrilateral.

[0178] After the discretized destruction mechanism is established, step S502 is executed to calculate the stability of the slope. The time interval for stability calculation is set. For example, in this embodiment, the stability of the slope is calculated every 7 days. Each calculation retrieves the transient simulation results of the corresponding time, such as node temperature, water content, effective saturation, ice content, and pressure head. The calculation of the stability of the slope is gradually advanced with the set time interval, and finally the evolution of the slope stability over time can be obtained. Figure 7(a) shows the change of the safety factor of the slope in this embodiment. Obviously, the strength of the soil increases after freezing in the cold season, resulting in an increase in the safety factor in the cold period, which then decreases in the melting season. Figure 7(b) shows a more detailed change in the safety factor in the melting season.

[0179] The embodiment of the present invention further provides a storage medium, which stores computer instructions and programs, and when the computer instructions and programs are read and run, the method for evaluating the stability evolution of unsaturated soil slopes in cold regions of the above embodiment is executed. The storage medium may include a memory, a flash memory, a register, or a combination thereof.

[0180] An embodiment of the present invention provides an electronic device, which may be a computer device or a server device. Figure 8, a schematic diagram of the structure of an electronic device. The electronic device includes a processor 10, a memory 11, and a bus 12. The processor 10 and the memory 11 are connected via the bus 12, and the processor 10 is used to execute an executable module stored in the memory 11, such as a computer program.

[0181] The processor 10 can be an integrated circuit chip with signal processing capabilities. In the implementation process, each step of the method for evaluating the stability evolution of unsaturated soil slopes in cold regions can be completed by hardware integrated logic circuits in the processor 10 or software instructions. The above-mentioned processor 10 can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components.

[0182] The memory 11 may include a high-speed random access memory (RAM), and may also include a non-volatile memory, such as at least one disk memory.

[0183] The bus 12 may be an ISA (Industry Standard Architecture) bus, a PCI (Peripheral Component Interconnect) bus, or an EISA (Extended Industry Standard Architecture) bus. Figure 8 Although only one bidirectional arrow is used in the figure, it does not mean that there is only one bus 12 or only one type of bus 12 .

[0184] The memory 11 is used to store programs, such as programs corresponding to the device for evaluating the evolution of stability of unsaturated soil slopes in cold regions. The device for evaluating the evolution of stability of unsaturated soil slopes in cold regions includes at least one software function module that can be stored in the memory 11 in the form of software or firmware or solidified in the operating system (OS) of the electronic device. After receiving the execution instruction, the processor 10 executes the program to implement the method for evaluating the evolution of stability of unsaturated soil slopes in cold regions.

[0185] Possibly, the electronic device provided by the embodiment of the present invention further includes a communication interface 13. The communication interface 13 is connected to the processor 10 via a bus.

[0186] It should be understood that Figure 8 The structure shown is only a schematic diagram of a portion of the electronic device. The electronic device may also include Figure 8 More or fewer components as shown, or with Figure 8 Different configurations shown. Figure 8 Each component shown in the figure can be implemented by hardware, software or a combination thereof.

[0187] In summary, the embodiments of the present invention provide a method, medium and device for evaluating the evolution of stability of unsaturated soil slopes in cold regions, establish a framework combining the soil-atmosphere interaction model and the water-vapor-heat coupling transmission model of unsaturated soil in cold regions, and develop a finite element transient simulation program that simulates the water-heat state response of unsaturated soil in cold regions under long- and short-term climate effects, which solves the technical problem that the latent heat release amount is easily lost due to the narrow phase change temperature range when simulating the freeze-thaw phenomenon of sandy soil with traditional numerical models. According to the upper limit theorem of limit analysis, a slope stability evaluation method is established using discretization recursive technology, and an innovative iterative solution combining the dichotomy method and the meta-heuristic optimization algorithm is proposed. A method of coupling the finite element transient analysis model and stability evaluation is proposed, which can be used to evaluate the evolution of stability of unsaturated soil slopes in cold regions driven by climate data.

[0188] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0189] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other specific forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations falling within the meaning and scope of the equivalent elements of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

Claims

1. A method for evaluating the stability evolution of unsaturated soil slopes in cold regions, characterized in that: The method comprises: S1, finite element modeling of unsaturated soil slope, including: generating a slope geometric model according to geometric control points and stratum control lines, dividing the slope geometric model into finite element meshes using a mesh generator, and assigning basic physical and mechanical parameters of the finite element meshes and nodes; S2, fitting the sub-daily weather data, including: extracting daily data of target meteorological variables from the meteorological data set, wherein the target meteorological variables include any one or more of air temperature, relative humidity, wind speed, sunshine hours and air pressure, and using a trigonometric function type weather data fitting program to generate hourly sequence weather data; S3, determining the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface, including: constructing a balance equation and a parameterized model based on the energy balance and mass balance of the soil-atmosphere interface, inputting the hourly series weather data obtained by S2 into the parameterized model, and solving the balance equation to obtain the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface; S4, water-vapor-heat multi-field coupled finite element simulation of unsaturated soil slope in cold region, including: constructing a water-vapor-heat multi-field coupled finite element simulation model of unsaturated soil slope in cold region, taking the heat and moisture flux exchanged at the soil-atmosphere interaction interface obtained from S3 as the upper boundary conditions, driving the water-vapor-heat multi-field coupled finite element simulation model of unsaturated soil slope in cold region to run, and obtaining the continuous changes of the internal hydrological and thermal states of unsaturated soil slope in cold region under the forcing of meteorological data after solving; S5, obtain the evolution information of the stability of unsaturated soil slopes in cold regions, including: S501 and S502; S501, establish a discretized failure mechanism generator based on the upper limit method of limit analysis, including: using the finite element grid and nodes established in S1 as the background, constructing the entire soil slope area as a Delaunay triangle diagram, using the discretized point-to-point recursive technology, generating potential failure mechanisms in the Delaunay triangle area, interpolating the hydrological and thermal state variable values ​​of the nodes obtained in step S4 to determine the hydrological and thermal state variable values ​​at the discrete points, and calculating the internal energy dissipation rate and external force power based on this, and then solving and outputting the sliding surface information and the energy ratio R corresponding to the potential sliding surface a , where the energy ratio R a is the ratio of internal energy dissipation rate to external force power; S502, construct an iterative solution for slope stability, use the strength reduction method to define the safety factor, and build a new iterative solution that combines the dichotomy method and the meta-heuristic optimization algorithm. At the beginning of each iteration, the dichotomy method is used to update the safety factor. At this time, the soil strength parameter is the soil strength parameter after the safety factor is reduced. The energy ratio R in S501 is used. a As the objective function, a meta-heuristic optimization algorithm is used to solve the optimal sliding surface and the corresponding energy ratio R under a given safety factor. a , until the safety factor determined by the dichotomy method reaches the minimum value, and the energy ratio R a When the difference between the value and 1 meets the threshold error, the iterative calculation is terminated and the safety factor and critical slip surface information are output; S503: Repeat the method described in S502 to calculate the slope safety factor and critical slip surface information over time, thereby obtaining the evolution information of the slope stability.

2. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 1, characterized in that: S2, the step of fitting sub-daily weather data, comprising: fitting the daily weather data of the target meteorological variable obtained from the meteorological data set into the sub-daily fluctuation of the weather data by using equations (1) to (3), wherein the sub-daily fluctuation is the weather data of the hourly sequence: Where t represents the current time, T a is the temperature corresponding to the current time t, T m is the daily average temperature of the day at the current time t, T max is the highest temperature of the day at the current time t, T min is the lowest temperature of the day at the current time t, RH a is the relative humidity of the air corresponding to the current time t, RH a,m is the daily average relative humidity of the air at the current time t, RH a,max is the maximum relative humidity of the day at the current time t, V a is the wind speed corresponding to the current time t, V m The daily average wind speed for the day of the current time t, V max The maximum wind speed of the day at the current time t, t max,T is the time when the daily maximum temperature occurs, t max,RHa is the time when the maximum relative humidity of the air occurs each day, t max,V The time when the maximum wind speed occurs each day.

3. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 1, characterized in that: S3 determines the heat and moisture flux exchanged at the soil-atmosphere interaction interface, and calculates the heat and moisture flux exchanged at the soil-atmosphere interaction interface by equations (4) and (5), respectively: G=R n -LE-H (4) q net =PRE / ρ l (5) Where G is the heat exchanged at the soil-atmosphere interaction interface, R n is the net radiation flux, LE is the latent heat flux, L is the latent heat of phase change from liquid water to gaseous water, E is the evaporation, H is the sensible heat flux, q net is the moisture flux exchanged at the soil-atmosphere interaction interface, P is the precipitation, R is the runoff, and ρ l is the density of liquid water.

4. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 1, characterized in that: The S4, water-vapor-heat multi-field coupled finite element simulation of unsaturated soil slope in cold regions, comprises the following steps: S401, establish the control equation of water-vapor-heat multi-field coupling of unsaturated soil: Taking unit volume soil as the research object, based on the principle of water mass conservation and energy conservation, establish the strong form control equation of the coupling system: Among them, ρ l is the density of liquid water, ρ i is the density of ice, ρ v is the density of gaseous water, θ l is the volume content of liquid water, θ i is the volume content of ice, θ a is the volume content of air, θ v is the volume content of water vapor, is the gradient operator, K lh is the isothermal liquid water migration coefficient, K lT is the non-isothermal liquid water migration coefficient, K vh is the isothermal water vapor diffusion coefficient, K vT is the non-isothermal water vapor diffusion coefficient, h is the suction head, y is the vertical spatial coordinate, q l is the amount of liquid water transferred, q v is the amount of water vapor transferred, C is the volume heat capacity of the soil, C l is the volumetric heat capacity of liquid water, C v is the volume heat capacity of water vapor, L0 is the volume latent heat of vaporization of liquid water, L f is the latent heat of melting of ice, T is the temperature, and λ is the effective thermal conductivity of the soil; S402, Degenerate treatment of the strong form governing equations of coupled systems: The van-Genuchten soil-water characteristic curve (SWCC) is expressed as formula (8), and its derivative relationship is formula (9): Among them, θ s is the volumetric water content at saturation, θ r is the volumetric water content in the residual state, α v 、n v and m v is the fitting parameter, S e is the effective saturation, and the derivative value C' is the specific water capacity; Secondly, the relationship between the introduction temperature and the low temperature suction head: Wherein, h0 is the pressure head when the soil is not frozen; HS is the Heaviside function; ω is the ratio of the water-ice interface free energy to the water-air interface free energy, which is 1 for colloidal soil particles and 2.2 for non-colloidal soil particles; T' is the freezing temperature; Combining equations (8) and (10), we establish the relationship between unfrozen liquid water and negative temperature (11), which is called the soil freezing characteristic curve: Let θ = θ l +ρ i θ i / ρ l , then the pore air content is expressed as the difference between the porosity n and θ: θ a =n-θ, substituting the above relationship into equation (6), we get the degenerated total water mass conservation equation: in, Based on equation (7), the same order and the same operator terms are merged to obtain the simplified form of the energy conservation equation: in, and The coupled system consisting of equations (12) and (13) can be described as a mathematical problem: Formula (14) represents the equation in a D-dimensional Euclidean space R with a closed boundary Γ. D In the physical domain Ω, at time t greater than or equal to t0, the state variables h, T and θ i Satisfied by The partial differential equation system composed of , in addition, on the boundary Γ, satisfies the following conditions: Dirichlet frontier: Newman Boundary: And at time t0, the initial variable value is: Where n is the unit normal vector outside the boundary, q and Q are the water flow rate and heat flow rate on the boundary, q w is the water flux exchanged at the soil-atmosphere interaction interface, q h heat exchanged at the soil-atmosphere interaction interface; S403, perform finite element spatial discretization on equation (14), and use Galerkin weighted residual method to integrate equation (14): The equivalent integral weak form of the partial differential equation for conservation of total water mass and energy can be obtained: Among them, L w is the residual of the equivalent integral equation of the total water mass, L h is the residual of the energy equivalent integral equation, N is the shape function, N T is the transpose of the shape function; For equation (17), firstly solve the unit integral to obtain the unit stiffness equation, and then assemble it to obtain the spatial discretization form of the moisture equation: Among them, the expressions of each coefficient matrix are: Among them, nel is the total number of units, n b is the number of nodes on the unit boundary to which the Newman boundary condition belongs; For equation (18), firstly solve the unit integral to obtain the unit stiffness equation, and then assemble it to obtain the spatial discretization form of the energy equation: Among them, the expressions of each coefficient matrix are: Equations (19) and (21) are processed by time discretization to obtain their conventional iterative formats. In order to overcome the problem that the conventional format is prone to lose latent heat information and it is difficult to ensure the conservation of total water mass due to its narrow temperature range when simulating the freezing and thawing phenomenon of sandy soil, an iterative level operator splitting method is introduced to transform the time term. Specifically, assuming that t is the previous time step, t+Δt is the current time step, the current iteration level is marked as j+1, the first iteration is 1, and the value of the first iteration level is inherited from the convergence value at time t. In the iterative calculation of [t, t+Δt], the time term is split into two terms, 1~j terms and j~j+1 terms: First, for the second term on the far right of equation (23), restore its original meaning, that is, the rate of change of the mass of liquid water and ice: Among them, M ww and M ii The expression is: Next, keep the second term on the right of equation (24) and equation (25) unchanged, and add the first term to obtain: Among them, M Ti The expression is: Coefficient M Ti is the sensible heat capacity, and its expression is: Substituting equations (23) and (26) into equation (19), and substituting equations (24), (25) and (28) into equation (21), combined with the backward Euler difference method, the final iterative format of the discretized equation is obtained: Among them, the continuous changes of the hydrological and thermal states of the unsaturated soil slopes in the cold region under the forcing of meteorological data include T t+ Δt ,j+1 and h t+Δt,j+1 ; T t+Δt,j represents the temperature value of the ramp node obtained after the jth iteration at the current time t+Δt, T t+Δt,j+1 represents the temperature value of the slope node obtained after the j+1th iteration at the current time t+Δt, h t+Δt,j represents the water head value of the slope node obtained after the jth iteration at the current time t+Δt, h t+Δt , j+1 represents the water head value of the slope node obtained after the j+1th iteration at the current time t+Δt, T t represents the temperature value of the ramp node after the previous time step t converges, θ t represents the total volume content of liquid water and ice at the slope node after convergence at the previous time step t, represents the ice volume content of the slope node obtained after the jth iteration at the current time t+Δt, represents the volume content of ice at the slope node after convergence at the previous time step t; S404, compile a finite element solution program, propose and adopt a single time step sequential solution scheme: that is, divide the current time step into two sub-steps, in the first sub-step, assume that the temperature is constant to solve the moisture discretization equation, update the pressure head and total water content; in the second sub-step, assume that the pressure head and total water content are constant to solve the energy discretization equation, update the temperature and ice content; after the two sub-steps are solved sequentially, iterate and loop the calculation until the convergence condition is met, output and store the results, the results include the T t+Δt,j+1 、h t+Δt,j+1 , as well as in, represents the ice volume content of the slope node at the current time t+Δt, after the j+1th iteration, represents the volume content of liquid water at the slope node obtained after the j+1th iteration at the current time t+Δt, By adding T t+Δt ,j +1 Substitute into formula (11) to calculate and determine: By adding h t+Δt,j+1 After conversion to water content Determine after subtraction.

5. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 1, characterized in that: The S501 is to establish a discretized failure mechanism generator based on the limit analysis upper limit method. The premise is that under given strength parameters, the establishment includes the following three steps: (1) A Delaunay triangle diagram is constructed with the finite element grid and nodes as the background. The node variable values ​​and soil strength parameters calculated in step S4 are assigned to the vertices of the Delaunay triangle. The slope toe point is taken as the starting point and the rotation radius r0 and the rotation angle θ0 are taken as the starting conditions. A series of discrete points are generated in sequence by using a recursive method and following the ideal plastic orthogonal condition. The ideal plastic orthogonal condition indicates that the angle between the velocity direction of the discrete point and the sliding surface is the friction angle at the discrete point. The node variable values ​​include the continuous changes of the internal hydrological and thermal states of the unsaturated soil slope in the cold region under the forcing of meteorological data. The soil strength parameters include the internal friction angle and cohesion. It is known that the previous discrete point P i The corresponding rotation angle β i , Radius of rotation r i , then the recursively generated point P i+1 The rotation angle β i+1 , Radius of rotation r i+1 and coordinates (x i+1 ,y i+1 ) are calculated using formulas (33) to (35) respectively: b i+1 =b i +d (33) Among them, (x A ,y A ) is the coordinate of the starting point A of the slope toe, which is known; δ is the angle between adjacent rotation radii, which is specified by the user; is a discrete point P i The internal friction angle at P is obtained by first finding i The position of the point in the Delaunay triangle graph returns the number of the triangle and the vertex number, and then uses the barycentric interpolation method to determine P i Angle of internal friction at point: Among them, κ, ξ and ζ are interpolation coefficients that vary with the position of the discrete point i in the triangle. and are the values ​​of the internal friction angles at the triangle vertices a, b, and c respectively; When the newly generated discrete point P i+1 When the slope line is exceeded, the recursive procedure will terminate and the last discrete point will be corrected to the slope line. Then, all discrete points are connected to form a potential slip surface. The slip line P between adjacent discrete points is i-1 P i Length L i Determined as: (2) It is also necessary to calculate the internal energy dissipation rate occurring on the slip surface and the power of all external forces. The internal energy dissipation rate only occurs on the slip surface. The total dissipation rate is calculated by summing the micro-segment dissipation rates: Where c is the cohesion, v is the velocity of the sliding surface, n is the total number of discretized segments of the sliding surface, dL is the length of the micro-segment of the sliding surface, and c i is the cohesive force at the discrete point. The pore fluid pressure also does work on the shear zone of the slip surface. The power of the pore fluid pressure is calculated by the micro-segment power summation method using formula (39): Among them, W u is the pore pressure power, u is the pore fluid pressure, u i is the pore fluid pressure at discrete points on the slip surface, v i is the movement velocity at discrete point i on the slip surface; The soil in different areas of the sliding body will have different densities and weights due to differences in components. Therefore, the closed area surrounded by two adjacent rotation radii, slope lines and sliding lines is divided into many small quadrilaterals. Each quadrilateral is divided into a series of smaller paired triangles. The weight γ of the soil at the centroid of each small triangle is determined by the barycentric interpolation method: c = p l gth l +r s g(1-n)+ρ a gth a +r i gth i (40) Where g is the acceleration due to gravity; After calculating the gravity power of all small triangles in the sliding area, the total gravity power is calculated by summing the micro-units: Among them, W g is the total gravitational power, N is the total number of small triangles, γ i is the weight at the centroid of the i-th small triangle, S i is the area of ​​the ith small triangle, v i is the speed of the ith small triangle (the direction is perpendicular to the line connecting the centroid and the rotation center), β′ i is the angle between the line connecting the centroid and the rotation center and the x-axis; (3) Define the energy ratio as: For a given strength parameter, when the calculated energy ratio reaches the minimum value, the corresponding slip surface is the optimal slip surface under the current conditions.

6. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 1, characterized in that: In S502, a novel iterative solution scheme combining the dichotomy method and the meta-heuristic optimization algorithm is used to solve the safety factor and the optimal sliding surface of the slope. The implementation steps are as follows: the safety factor (FoS) is defined by the strength reduction method. After each iteration, a new safety factor is obtained by the dichotomy method. Then, the original strength parameters are reduced by the safety factor to obtain a set of new cohesion and internal friction angles. The new cohesion, internal friction angle and other state variable values ​​are input into the energy ratio R a It is the optimization procedure of the objective function. After iterative calculation, the energy ratio R a When it is the minimum value under the given strength parameter, the iteration ends, and a smaller safety factor is calculated by dichotomy, and the next iteration starts until the minimum value of the safety factor is found. The safety factor value at this time and the critical slip surface information under the safety factor value are output. Among them, other state variable values ​​include temperature, fluid pressure, water content and ice content.

7. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 6, characterized in that: The meta-heuristic optimization algorithm is used to determine the minimum energy ratio R under a given safety factor. a and the corresponding slip surface. The variables of the optimization program are the starting rotation radius r0 and the starting rotation angle θ0 with specified upper and lower limits. The variable seeds are randomly arranged according to the optimization algorithm. Under each group of variable combinations, the S501 program is called to generate potential failure mechanisms, thereby generating different potential slip surfaces. The optimal slip surface under the current safety factor is obtained through optimization.

8. The method for evaluating the stability evolution of unsaturated soil slopes in cold regions according to claim 6, characterized in that: In the new iterative solution scheme combining the bisection method and the meta-heuristic optimization algorithm, the reduced strength parameters include the cohesion c' and the safety factor φ', which are calculated using formula (43): The bisection method shown in formula (44) is used to update the safety factor after each iteration: Among them, FoS ub and FoS lb It is a repository of safety factors. When the iterative FoS is greater than or equal to 1, it is stored in the database FoS ub When the iteratively obtained FoS is less than 1, it is stored in the database FoS lb middle.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 8 is implemented.

10. An electronic device, characterized in that: include: A processor and a memory, the memory being used to store one or more programs; When the one or more programs are executed by the processor, the method according to any one of claims 1 to 8 is implemented.

Images