A method, medium and equipment for evaluating the stability evolution of unsaturated soil slopes in cold regions

Abstract

The present invention proposes a method, medium, and equipment for evaluating the stability evolution of unsaturated soil slopes in cold regions. It establishes a framework that combines a soil-atmosphere interaction model and a water-vapor-heat coupled transmission model for unsaturated soil in cold regions. It 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. This solves the technical problem that when traditional numerical models simulate the freeze-thaw phenomenon of sandy soil, the latent heat release amount is easily lost due to the narrow phase change temperature range. Based on the upper limit theorem of limit analysis, a slope stability assessment method is established using discretization recursive technology, and an innovative iterative solution scheme combining a dichotomy method and a meta-heuristic optimization algorithm is proposed. A method for coupling a finite element transient analysis model with stability assessment is proposed, which can be used to evaluate the stability evolution of unsaturated soil slopes in cold regions driven by climate data.

Description

Technical Field

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

[0002] In recent years, with the frequent occurrence of abnormal global weather and the surge in human socioeconomic activities, a large number of landslides and slope failures have occurred in mountainous areas and along transportation infrastructure in cold regions. The evolutionary process of geological hazards (gestation, formation, movement, and disaster-causing) is closely related to the coupled interactions among various spheres of the Earth system. The gradual loss of strength of near-surface soils under the influence of seasonal climatic conditions triggers a series of natural disasters such as landslides and collapses.

[0003] However, few studies have directly linked dynamic climate change to slope stability assessment. Most studies on the impact of climate on slope stability are based on narrative development and empirical evaluation. Current analytical methods and related computational software typically rely on snapshots of the hydrological state of the soil at a specific time period, failing to reflect the evolution of slope stability under the influence of dynamic climate systems. 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, the stability assessment based on traditional methods will not match the actual state. Furthermore, most stability analyses are conducted only on slopes in areas without freeze-thaw activity. Studies on slopes in cold regions tend to use qualitative analysis, and relevant computational theories and methods are particularly scarce. The few quantitative assessments that have been conducted almost entirely use infinite slope models and assume the distribution of ground ice and hydrological conditions.

[0004] The above factors seriously restrict the assessment of long-term slope stability in cold areas under climate influence. 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 objectives, the technical solutions adopted in the embodiments of the present invention are 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 based on geometric control points and stratigraphic control lines, dividing the slope geometric model into finite element meshes using a mesh generator, and assigning basic physical and mechanical parameters to the finite element meshes and nodes;

[0009] S2, performing sub-daily weather data fitting, 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 generating hourly series weather data using a trigonometric function-based weather data fitting program;

[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 in 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 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 fluxes exchanged at the soil-atmosphere interaction interface obtained in S3 as upper boundary conditions to drive the operation of 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;

[0012] S5, obtain the evolution information of unsaturated soil slope stability 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 mesh and nodes established in S1 as the background, constructing the entire slope area into a Delaunay triangle diagram, using the discretized point-to-point recursive technology, generating potential failure mechanisms within 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 slip surface information and the energy ratio R corresponding to the potential slip 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 bisection method and the meta-heuristic optimization algorithm. At the beginning of each iteration, the bisection 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 bisection 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 slip 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, which 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 existing technology, the embodiment of the present invention provides a method, medium and equipment 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-thermal state response of unsaturated soil in cold regions under long-term 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 sand. 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 scheme combining the dichotomy method and 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 following briefly introduces the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any 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 according to an embodiment of the present invention;

[0026] FIG5( b ) is a diagram showing changes in evaporation on a slope surface according to 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 the 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 by 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 C-C′ provided by an embodiment of the present invention;

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

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

[0032] Figure 8 A schematic structural diagram 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] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions of the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only 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 herein 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 as claimed, but rather merely represents selected embodiments of the present invention. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort shall fall 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. Therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are used only to distinguish the description and should not be understood as indicating or implying relative importance.

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

[0038] In the description of the present invention, it should be noted that the terms "upper", "lower", "inside", "outside", etc. indicate orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, or are the orientations or positional relationships in which the inventive product 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 orientation, be constructed and operated in a specific orientation. Therefore, they should not be understood as limiting the present invention.

[0039] In the description of the present invention, it should also be noted that, unless otherwise expressly specified or limited, the terms "disposed" and "connected" should be understood in a broad sense. For example, they can refer to fixed connections, detachable connections, or integral connections; mechanical connections, or electrical connections; direct connections, indirect connections through an intermediate medium, or internal connections between two components. Those skilled in the art will understand the specific meanings of the above terms in the present invention based on the specific circumstances.

[0040] The following embodiments of the present invention are described in detail with reference to the accompanying drawings. In the absence of conflict, the following embodiments and features in the embodiments may 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 stability evolution of unsaturated soil slopes in cold regions, so as to evaluate the temporal evolution of the stability of unsaturated soil slopes in cold regions under the influence of climate. It innovatively couples the climate data-driven water-vapor-heat coupled 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 embodiment of the present invention provides a method for evaluating the stability evolution of unsaturated soil slopes in cold regions, which 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 stratigraphic control lines, dividing the slope geometric model into finite element meshes using a mesh generator, and assigning basic physical and mechanical parameters to the finite element meshes and nodes.

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

[0045] S2, fitting 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 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 information such as the longitude, latitude, and altitude 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 in 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 of 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 of unsaturated soil slopes in cold regions to run, and after solving, obtaining the continuous changes in the internal hydrological and thermal states of unsaturated soil slopes in cold regions under the forcing of meteorological data.

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

[0050] S5, obtain the evolution information of unsaturated soil slope stability 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 mesh and nodes established in S1 as the background, constructing the entire slope area into a Delaunay triangle diagram, using the discretized point-to-point recursive technology, generating potential failure mechanisms within 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 slip surface information and the energy ratio R corresponding to the potential slip 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 bisection method and the meta-heuristic optimization algorithm. At the beginning of each iteration, the bisection 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 bisection 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 slip 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 coupled transmission model of unsaturated soil in cold regions, and develops a finite element transient simulation program that simulates the water-thermal 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 sand. 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 scheme combining the dichotomy method and the meta-heuristic optimization algorithm is proposed. A method of coupling the finite element transient analysis model with 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 exhibit 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 dataset to the sub-daily fluctuations of the weather data using equations (1) to (3), where the sub-daily fluctuations are hourly series 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 maximum 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 The daily average relative humidity of the air at the current time t, RH a,max is the maximum relative humidity of the air on the day of the current time t, V a is the wind speed corresponding to the current time t, V m The daily average wind speed at 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 determines the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface, and calculates 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 amount, H is the sensible heat flux, q net is the water flux exchanged at the soil-atmosphere interaction interface, P is precipitation, R is runoff, and ρ l is the density of liquid water.

[0062] Net radiation flux is obtained by summing shortwave and longwave radiation. Shortwave radiation can be calculated using theoretical formulas based on local latitude and longitude, calculation time, surface soil moisture, and sunshine duration. Longwave radiation is calculated using the Stefan-Boltzmann radiation formula (which requires air and surface soil temperatures) based on air emissivity (calculated from relative humidity and air temperature) and surface emissivity (surface soil moisture). Other energy fluxes, such as L and H, are determined using the Monin-Obukhov similarity theory.

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

[0064] S401, establish the governing equations for the water-vapor-heat multi-field coupling of unsaturated soil: Taking unit volume soil as the research object, based on the principles of water mass conservation and energy conservation, establish the strong form governing equations for the coupled 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 fraction 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 migration of liquid water, q v is the amount of water vapor migration, C is the volume heat capacity of the soil, and C l is the volumetric heat capacity of liquid water, C v is the volumetric heat capacity of water vapor, L0 is the volumetric 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-gas 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 Equation (8), and its derivative relationship is Equation (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 below 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] Where h0 is the pressure head before the soil freezes (specifically the matrix suction head before the unsaturated soil freezes); 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 formula (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 formula (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 consisting of equations (12) and (13) can be described as the following mathematical problem:

[0085]

[0086]

[0087] Formula (14) represents the equation of the D-dimensional Euclidean space R with 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 system of partial differential equations 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 vector normal to 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, perform finite element spatial discretization on Equation (14), and integrate Equation (14) using the Galerkin weighted residual method to obtain the equivalent integral weak form of the partial differential equations for total water mass and energy conservation:

[0093]

[0094] Among them, L w is the residual of the equivalent integral equation of the total mass of water, 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), first solve the unit integral to obtain the unit stiffness equation, and then assemble it to obtain 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), first solve the integral in the unit to get the unit stiffness equation, and then assemble it to get 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 time discretized to obtain their conventional iterative formats. To overcome the problems of the conventional format simulating the freeze-thaw phenomenon of sand soil, which is prone to losing latent heat information and difficult to ensure the conservation of total water mass due to its narrow temperature range, 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 and j~j+1:

[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 terms on the right of Equations (24) and (25) unchanged, and add the first terms 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), and combining the backward Euler difference method, we obtain the final iterative format of the discretized equation:

[0124]

[0125]

[0126] Among them, the continuous changes of hydrological and thermal states inside unsaturated soil slopes in cold regions 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 ramp 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 solving the two sub-steps sequentially, iterate the calculation cycle until the convergence condition is met, output and store the results, and the results include the solved T t+Δt,j+1 、h t+Δt,j+1 、 as well as in, represents the ice volume content of the slope node obtained after the j+1th iteration at the current time t+Δt, Indicates the volume content of liquid water at the slope node at the current time t+Δt, obtained after the j+1th iteration, By changing T t+Δt,j+1 Substitute into formula (11) to calculate and determine, By h t+Δt,j+1 After conversion to water content Determine after subtraction.

[0128] Optionally, S501, a discretized failure mechanism generator based on the upper limit method of limit analysis is established. The prerequisite is that under given strength parameters, the 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 used as the starting point and the rotation radius r0 and rotation angle θ0 are used as the starting conditions. A series of discrete points are generated in sequence using the recursive method and the ideal plastic orthogonality condition. The ideal plastic orthogonality condition means that the angle between the velocity direction of the discrete point and the slip surface is the friction angle at the discrete point. The node variable values ​​include the continuous changes in the internal hydrological and thermal state 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 , rotation radius r i , then the recursively generated point P i+1 The rotation angle β i+1 , rotation radius 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 at the toe of the slope, 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 graph 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 internal friction angle values ​​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 as follows: i-1 P i 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 slip surface, n is the total number of discretized segments of the slip surface, dL is the length of the micro-segment of the slip surface, and c is the length of the micro-segment of the slip surface. 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 using formula (39) by summing the micro-segment power:

[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 velocity at discrete point i on the slip surface;

[0142] The soil in different areas of the sliding body will have different densities and specific gravity due to differences in composition. 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 specific gravity γ 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 i-th small triangle, v i is the speed of the i-th 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 a bisection method and a metaheuristic optimization algorithm is used to solve the safety factor and the optimal sliding surface of the slope. The implementation steps are as follows: a strength reduction method is used to define the safety factor (FoS). After each iteration, a new safety factor is obtained by bisection method. Then, the original strength parameters are reduced by the safety factor to obtain a new set of cohesion and internal friction angles. The new cohesion, internal friction angles and other state variable values ​​are input into the energy ratio R. a It is the optimization procedure of the objective function. When the energy ratio Ra When the minimum value under the given strength parameter is reached, the current iteration ends, and a smaller safety factor is calculated by bisection method, and the next iteration begins. When the minimum value of the safety factor is found, the safety factor value at that 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 The variables of the optimization program are the starting rotation radius r0 and the starting rotation angle θ0 with the upper and lower limit ranges. The variable seeds are randomly arranged according to the optimization algorithm. Under each set 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] Alternatively, in a novel iterative solution scheme combining bisection and metaheuristic optimization algorithms, the reduced strength parameters include cohesion c' and safety factor Calculate using formula (43):

[0154]

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

[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: Establish a finite element model of unsaturated soil slope;

[0160] The embodiment is a high-speed railway embankment slope in a certain area, with a slope ratio of 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 also 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 with a setting of 0.03W / m 2A constant geothermal flux is maintained, and the lateral boundaries are insulated and impermeable. Monitoring profiles were established 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 using step S2 are shown. 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 steps S1 and S2 are loaded into the finite element transient analysis model in file format. This drives the finite element program, which first executes step S3 to calculate the heat and moisture fluxes exchanged between the soil and the atmosphere. Figure 5(a) shows the diurnal variation series of various energy fluxes on the slope surface. Given the low rainfall and runoff in the area, Figure 5(b) only shows the diurnal variation series of evaporation.

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

[0170] Then, in the finite element transient model, the flux calculated in step S3 is used as the upper boundary condition. The driver program continues to run, and according to the two-step solution described in step S4, the hydrological state change of the underlying unsaturated soil is first calculated, followed by the thermal state change. This cycle iterates until the current time step converges. The converged value is stored, the time is updated, and the calculation of the next time step begins. This process is repeated until the end of the time period is reached, and the calculation ends. Please refer to Figures 6(a), 6(b), and 6(c). Figure 6 shows the temperature contours 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 in 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 conditions 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 slip surface must be the friction angle at this discrete point. For unfrozen soil, the friction angle remains unchanged. For frozen soil, the frozen soil strength model proposed by Nishiruma and Wang (2019) is used. The friction angle and cohesion are related to the low-temperature suction and can be determined using 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 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 failure mechanism is established, step S502 is executed to calculate the stability of the slope. The time interval for the 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] An embodiment of the present invention further provides a storage medium storing computer instructions and programs that, when read and executed, execute the method for evaluating the stability evolution of unsaturated soil slopes in cold regions described in the above embodiment. The storage medium may include memory, flash memory, registers, 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. The processor 10 is used to execute executable modules stored in the memory 11, such as computer programs.

[0181] The processor 10 can be an integrated circuit chip with signal processing capabilities. During 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 Only one bidirectional arrow is used in the figure, but it does not mean that there is only one bus 12 or one type of bus 12.

[0184] The memory 11 is used to store programs, such as a program corresponding to an apparatus for evaluating the stability evolution of unsaturated soil slopes in cold regions. The apparatus for evaluating the stability evolution of unsaturated soil slopes in cold regions includes at least one software functional module that can be stored in the memory 11 in the form of software or firmware or embedded in the operating system (OS) of the electronic device. Upon receiving an execution instruction, the processor 10 executes the program to implement the method for evaluating the stability evolution 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 than 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 equipment for evaluating the stability evolution 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-thermal 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 sand. 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 scheme combining the dichotomy method and the meta-heuristic optimization algorithm is proposed. A method of coupling the finite element transient analysis model with stability evaluation is proposed, which can be used to evaluate the stability evolution of unsaturated soil slopes in cold regions driven by climate data.

[0188] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Those skilled in the art will readily appreciate that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention are intended to be within the scope of protection of the present invention.

[0189] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above and that the invention can be embodied in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered in all respects as illustrative and non-restrictive, and the scope of the invention is defined by the appended claims, not the foregoing description, and all variations within the meaning and range of equivalents of the claims are intended to be included therein. Any reference sign in a claim should not be construed 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 based on 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 to the finite element meshes and nodes; S2, performing sub-daily weather data fitting, 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 generating hourly series weather data using a trigonometric function-based weather data fitting program; 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 in 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 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 fluxes exchanged at the soil-atmosphere interaction interface obtained in S3 as upper boundary conditions to drive the operation of 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; S5, obtain the evolution information of unsaturated soil slope stability 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 mesh and nodes established in S1 as the background, constructing the entire slope area into a Delaunay triangle diagram, using the discretized point-to-point recursive technology, generating potential failure mechanisms within 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 slip surface information and the energy ratio R corresponding to the potential slip 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 bisection method and the meta-heuristic optimization algorithm. At the beginning of each iteration, the bisection 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 bisection 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, includes: fitting the daily weather data of the target meteorological variable obtained from the meteorological data set into sub-daily fluctuations of the weather data by using equations (1) to (3), wherein the sub-daily fluctuations are hourly series weather data: 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 maximum 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 The daily average relative humidity of the air at the current time t, RH a,max is the maximum relative humidity of the air on the day of the current time t, V a is the wind speed corresponding to the current time t, V m The daily average wind speed at 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: The S3 determines the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface, and the heat and moisture fluxes exchanged at the soil-atmosphere interaction interface are calculated 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 amount, H is the sensible heat flux, q net is the water flux exchanged at the soil-atmosphere interaction interface, P is precipitation, R is 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 slopes in cold regions, includes the following steps: S401, establish the governing equations for the water-vapor-heat multi-field coupling of unsaturated soil: Taking unit volume soil as the research object, based on the principles of water mass conservation and energy conservation, establish the strong form governing equations for the coupled 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 fraction 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 migration of liquid water, q v is the amount of water vapor migration, C is the volume heat capacity of the soil, and C l is the volumetric heat capacity of liquid water, C v is the volumetric heat capacity of water vapor, L0 is the volumetric 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 is: Where 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 can 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 formula (6), we get the degenerated total water mass conservation equation: in, Based on formula (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 of the D-dimensional Euclidean space R with 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 system of partial differential equations 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 vector normal to 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 integrate Equation (14) using the Galerkin weighted residual method. 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 mass of water, 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 it in the unit integration to get the unit stiffness equation, and then assemble it to get 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), first solve it in the unit integration to get the unit stiffness equation, and then assemble it to get the spatial discretization form of the energy equation: Among them, the expressions of each coefficient matrix are: Equations (19) and (21) are time discretized to obtain their conventional iterative formats. To overcome the problems of the conventional format simulating the freeze-thaw phenomenon of sand soil, which is prone to losing latent heat information and difficult to ensure the conservation of total water mass due to its narrow temperature range, 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 and j~j+1: 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 terms on the right of Equations (24) and (25) unchanged, and add the first terms 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), and combining the backward Euler difference method, we obtain the final iterative format of the discretized equation: Among them, the continuous changes of hydrological and thermal states inside unsaturated soil slopes in cold regions 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 ramp 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, and 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, and update the temperature and ice content; after solving the two sub-steps sequentially, iterate and loop the calculation until the convergence condition is met, output and store the results, which 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 obtained after the j+1th iteration at the current time t+Δt, Indicates the volume content of liquid water at the slope node at the current time t+Δt, obtained after the j+1th iteration, By putting 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 above S501 is to establish a discretized failure mechanism generator based on the upper limit method of limit analysis. The premise is that under given strength parameters, the establishment includes the following three steps: (1) Construct a Delaunay triangle diagram with the finite element grid and nodes as the background, assign the node variable values ​​and soil strength parameters calculated in step S4 to the vertices of the Delaunay triangle, take the slope toe point as the starting point and the rotation radius r0 and rotation angle θ0 as the starting conditions, and use the recursive method and the ideal plastic orthogonality condition to generate a series of discrete points in sequence. The ideal plastic orthogonality condition means 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 state 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 , rotation radius r i , then the recursively generated point P i+1 The rotation angle β i+1 , rotation radius 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 internal friction angle values ​​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 slip surface, n is the total number of discretized segments of the slip surface, dL is the length of the micro-segment of the slip surface, and c is the length of the micro-segment of the slip surface. 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 using formula (39) by summing the micro-segment power: 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 velocity at discrete point i on the slip surface; The soil in different areas of the sliding body will have different densities and specific gravity due to differences in composition. 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 specific gravity γ 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 i-th small triangle, v i is the velocity 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 combining the bisection method and the metaheuristic 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 bisection method. Then, the original strength parameters are reduced by the safety factor to obtain a new set of 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. When the energy ratio R a When the minimum value under the given strength parameter is reached, the current iteration ends, and a smaller safety factor is calculated by bisection method, and the next iteration begins. When the minimum value of the safety factor is found, the safety factor value at that 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: A meta-heuristic optimization algorithm is used to determine the minimum energy ratio R under a given safety factor. a The variables of the optimization program are the starting rotation radius r0 and the starting rotation angle θ0 with the upper and lower limit ranges. The variable seeds are randomly arranged according to the optimization algorithm. Under each set 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 novel iterative solution scheme combining the bisection method and the metaheuristic optimization algorithm, the reduced strength parameters include the cohesion c' and the safety factor φ', which are calculated using Equation (43): The safety factor after each iteration is updated using the bisection method shown in formula (44): 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 configured 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