Risk simulation method for high-cold and high-altitude freeze-thaw slope based on region identification

Abstract

The present application relates to the technical field of cold region geotechnical engineering and geological disaster prevention, and particularly relates to a high-cold high-altitude freeze-thaw slope risk simulation method based on regional identification, which comprises the following steps: discretizing the slope and initializing the particle group by using the SPH method; calculating the physical field evolution index, and identifying the phase transition active area, the shear damage area and the stable inert area by threshold comparison; solving the heat conduction and water migration equation to obtain the temperature field, phase transition interface and water field data; performing adaptive adjustment on the particles according to the regional identification result, and calculating the stress response in combination with the dynamic freeze-thaw constitutive model; solving the momentum equation to update the motion state and correct the permeability coefficient, determining the time step according to the Courant condition and performing cyclic calculation. The present application realizes adaptive adjustment of the particle resolution through regional identification, takes into account the accuracy and calculation efficiency of the freeze-thaw interface tracking, and realizes fine simulation of the multi-field coupling process of the freeze-thaw slope.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering and geological disaster prevention technology in cold regions, and in particular to a method for simulating the risk of freeze-thaw slopes in high-altitude and cold regions based on regional identification. Background Technology

[0002] Freeze-thaw cycle-induced slope instability involves a strongly nonlinear coupling process of water migration, phase change frost heave, and mechanical damage. Existing numerical simulations mostly employ mesh-based methods such as the finite element method, which are prone to mesh distortion when dealing with large deformations caused by frost heave, and struggle to efficiently track moving freezing fronts. Global refinement also faces the problem of high computational costs. While smoothed particle hydrodynamics methods can overcome the mesh distortion problem, current techniques typically use uniform particle distributions, failing to adaptively adjust the resolution based on the evolution characteristics of the physical field, resulting in insufficient computational accuracy in critical areas such as phase change interfaces or shear zones. Furthermore, existing methods often neglect the permeability evolution caused by soil shear damage, making it difficult to accurately describe the water enrichment resulting from damage and its feedback mechanism that exacerbates frost heave.

[0003] Therefore, there is an urgent need for a simulation method that can adaptively adjust the calculation accuracy and couple the damage-seepage effect. Summary of the Invention

[0004] Therefore, this invention provides a risk simulation method for high-altitude and cold-climate freeze-thaw slopes based on regional identification, in order to solve the aforementioned problems existing in the prior art.

[0005] To achieve the above objectives, this invention provides a method for simulating the risk of freeze-thaw slopes in high-altitude and cold regions based on regional identification, comprising:

[0006] Step S1: Based on the slope geological survey data, the computational domain is discretized using the smooth particle hydrodynamics method, and the particle swarm carrying mass, position, temperature, volume of unfrozen water content, porosity and effective stress tensor is obtained through gravity initialization.

[0007] Step S2: Calculate the physical field evolution index in the particle neighborhood based on the current state variables of the particle swarm, and compare the physical field evolution index with a preset threshold to obtain the region identification result.

[0008] Step S3: Based on the region identification results and environmental boundary conditions, the heat conduction equation and moisture migration equation are solved using the SPH discrete scheme that considers the latent heat source of phase change, so as to obtain the temperature field update data, phase change interface location data and moisture field distribution data at the current moment.

[0009] Step S4: Based on the temperature field update data, the variable interface position data, and the region identification results, an adaptive adjustment strategy is executed on the particles of the particle swarm, and the freeze-thaw constitutive model is combined to calculate the freeze-thaw strain increment and stress response to obtain the current particle resolution, current stress tensor, and current plastic state.

[0010] Step S5: Solve the momentum conservation equation to obtain particle motion data based on the current stress tensor and the current plastic state, and correct the soil permeability coefficient based on the current plastic state.

[0011] Step S6: Based on the particle motion data, determine the time step of the next calculation step using the Courant condition, and repeat steps S2 to S5 until the preset freeze-thaw cycle is completed.

[0012] Furthermore, the process of step S2 includes:

[0013] Based on the temperature state variable carried by the particle, the temperature gradient modulus in the particle's neighborhood is calculated, and the region where the temperature gradient modulus is greater than a first threshold is identified as an active phase transition region.

[0014] Based on the effective stress tensor carried by the particle, the equivalent plastic strain increment of the particle is calculated, and the region where the particle with the equivalent plastic strain increment is greater than the second threshold is identified as the shear damage zone.

[0015] The regions containing particles that do not belong to the active phase transition region and the shear damage region are identified as stable inert regions.

[0016] Furthermore, the process of calculating the temperature gradient modulus within the particle's neighborhood includes:

[0017] The SPH kernel function approximation method is adopted. Based on the temperature, mass, density and kernel function gradient vector of the particle and its neighboring particles, the temperature gradient tensor is calculated by discrete summation. The norm of the tensor is taken as the temperature gradient modulus.

[0018] Furthermore, the process of step S3 includes:

[0019] Construct environmental boundary entity particles, and assign temperature time history or water potential state to the environmental boundary entity particles according to the environmental boundary conditions.

[0020] Based on the region identification results, the thermal conductivity and permeability coefficients of particles in different regions are dynamically assigned corresponding values.

[0021] The heat flux and moisture flux between a particle and its neighboring particles and the boundary particles of the environment are calculated using the discrete scheme of the SPH kernel function. The temperature field update data and moisture field distribution data at the current moment are obtained by solving the problem.

[0022] Furthermore, the process of calculating the heat flux and moisture flux between a particle and its neighboring particles and environmental boundary particles based on the SPH kernel function discrete scheme includes:

[0023] In the heat conduction calculation, a discrete equation is constructed based on the temperature difference between the calculated particle and its neighboring particles, thermal conductivity, mass, density, and kernel function gradient. A latent heat source term for phase change is introduced to solve for the rate of temperature change. The thermal conductivity is dynamically assigned based on the particle position and region identification results.

[0024] In the calculation of water migration, a discrete equation is constructed based on the pore pressure difference between the calculated particle and its neighboring particles, the permeability coefficient, and the kernel function gradient. An anti-singularity coefficient is introduced to correct the calculation singularity when the particle spacing is too close, and the water flux is obtained by solving the equation.

[0025] Furthermore, the process of implementing an adaptive adjustment strategy for the particles in the particle swarm includes:

[0026] For particles identified as active phase transition regions by the region identification results, when their temperature gradient modulus exceeds the first threshold, a particle splitting operation is performed to decompose a single parent particle into a preset number of sub-particles. Based on the principles of mass conservation and momentum conservation, the sub-particles inherit the mass, velocity, temperature, and effective stress tensor of the parent particle. At the same time, the smooth length of the sub-particles is reduced accordingly to improve the local computational resolution.

[0027] For particles identified as stable inert regions by the region identification results, when the particle spacing is less than a preset merging threshold and the difference in physical field variables is lower than a limit value, a particle merging operation is performed to aggregate neighboring particles and reconstruct them into a single coarse particle.

[0028] Furthermore, the process of calculating the freeze-thaw strain increment and stress response in conjunction with the dynamic freeze-thaw constitutive model to obtain the current particle resolution, current stress tensor, and current plastic state includes:

[0029] Based on the temperature change and phase transition interface location data in the temperature field update data, combined with the changes in soil porosity and unfrozen water content, the frost heave strain increment is calculated using the frost heave coefficient.

[0030] The total strain increment is calculated based on particle motion data, and the thermal strain increment and the frost heave strain increment are subtracted from the total strain increment to obtain the effective mechanical strain increment.

[0031] The effective mechanical strain increment is calculated using the elastic constitutive matrix to obtain the test stress tensor, and the cohesion and internal friction angle parameters in the yield criterion are dynamically updated based on the current temperature.

[0032] The radial back mapping algorithm is used to plastically correct the experimental stress tensor, and the stress state is mapped back to the updated yield surface to obtain the current stress tensor and the current plastic state.

[0033] Furthermore, the process of step S5 includes:

[0034] Based on the current stress tensor, the stress divergence term in the particle's neighborhood is calculated using the SPH kernel function approximation method. The momentum conservation equation is then constructed by combining the gravitational volume force term, and the particle acceleration is obtained by solving it.

[0035] The particle acceleration is integrated using an explicit time integration algorithm to update the particle's velocity and position data;

[0036] Extract the cumulative equivalent plastic strain in the current plastic state, substitute it into the preset permeability coefficient evolution model, and dynamically correct the soil permeability coefficient.

[0037] Furthermore, the process of dynamically correcting the soil permeability coefficient includes:

[0038] Damage variables are introduced to describe the evolution of soil micropore structure, and a mapping relationship between cumulative equivalent plastic strain and damage variables is established.

[0039] The change in porosity is calculated based on the damage variables, and the permeability coefficient correction factor at the current moment is calculated using the exponential damage permeability formula.

[0040] The initial permeability coefficient is multiplied by the permeability correction factor to obtain the corrected soil permeability coefficient, which is used for water migration calculation at the next time step.

[0041] Furthermore, the process of determining the time step size of the next calculation step using the Courant condition includes:

[0042] The critical step size for mechanical stability was calculated based on the interparticle spacing and the current sound velocity, the critical step size for thermal stability was calculated based on the thermal diffusivity, and the critical step size for water migration stability was calculated based on the corrected permeability coefficient.

[0043] The minimum value among the mechanical stability critical step size, thermal stability critical step size, and moisture migration stability critical step size is selected and multiplied by a preset safety factor to serve as the time step size for the next calculation step.

[0044] Determine whether the current accumulated time has reached the preset freeze-thaw cycle period. If it has, stop the calculation and output the slope displacement field, stress field and plastic zone distribution data.

[0045] Compared with the prior art, the beneficial effects of the present invention are as follows: the present invention accurately identifies the active phase transition region and the shear damage region by calculating the temperature gradient modulus and the equivalent plastic strain increment, and drives the particle resolution to adaptively adjust according to the intensity of the physical field change, so that the computational node density is accurately matched with the spatial characteristics of the thermodynamic gradient and mechanical damage evolution, effectively solving the contradiction between computational accuracy and efficiency under fixed resolution; at the same time, it uses the cumulative plastic strain to describe the pore structure damage and corrects the permeability coefficient in real time, and establishes a dynamic feedback mechanism between mechanical damage and water migration, thereby realizing a high-fidelity simulation of the multi-field coupled evolution process of freeze-thaw slopes. Attached Figure Description

[0046] Figure 1 A flowchart illustrating the risk simulation method for high-altitude and cold-climate freeze-thaw slopes based on region identification provided by this invention;

[0047] Figure 2 A flowchart illustrating step S2 in the method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification provided by the present invention.

[0048] Figure 3 This is a flowchart illustrating step S3 in the method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification provided by the present invention. Detailed Implementation

[0049] To make the objectives and advantages of the present invention clearer, the present invention will be further described below with reference to embodiments; it should be understood that the specific embodiments described herein are merely for explaining the present invention and are not intended to limit the present invention.

[0050] Preferred embodiments of the present invention will now be described with reference to the accompanying drawings. Those skilled in the art should understand that these embodiments are merely illustrative of the technical principles of the present invention and are not intended to limit the scope of protection of the present invention.

[0051] It should be noted that in the description of this invention, the terms "upper", "lower", "left", "right", "inner", "outer", etc., which indicate directions or positional relationships, are based on the directions or positional relationships shown in the accompanying drawings. This is only for the convenience of description and is not intended to indicate or imply that the device or element must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, it should not be construed as a limitation of this invention.

[0052] Furthermore, it should be noted that, in the description of this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0053] Please see Figure 1 As shown, this invention provides a method for simulating the risk of freeze-thaw slopes in high-altitude and cold regions based on regional identification, including:

[0054] Step S1: Based on the slope geological survey data, the computational domain is discretized using the smooth particle hydrodynamics method, and the particle swarm carrying mass, position, temperature, volume of unfrozen water content, porosity and effective stress tensor is obtained through gravity initialization.

[0055] Specifically, based on the slope geological survey data, the geometric boundary of the computational domain and the soil layer distribution information of the slope are determined, and the physical and mechanical parameters of each soil layer are assigned values. These physical and mechanical parameters include density, elastic modulus, Poisson's ratio, cohesion, internal friction angle, thermal conductivity, specific heat capacity, latent heat of phase change, and initial moisture content. The initial particle spacing Δx is set according to the required computational accuracy. Based on the Smooth Particle Hydrodynamics (SPH) method, the continuous computational domain is discretized into a series of independent particles carrying physical information. The smoothing length h is determined, typically taken as h = 1.2Δx ~ 1.3Δx to ensure sufficient supporting particles in the neighborhood. The mass m of each particle is calculated based on the soil density ρ and the particle spacing, using the following formula: Where d is the computational dimension; each particle is assigned initial position coordinates and its initial velocity is set to zero; based on the ambient temperature and ground temperature distribution data during the exploration, the particles are assigned an initial temperature T0, and the corresponding initial unfrozen water content is calculated based on the unfrozen water content function; based on the initial porosity of the soil obtained from the geological exploration... Calculate and assign initial porosity to particles : Based on the soil pore characteristics, the particles are assigned an initial permeability coefficient k0; gravitational acceleration is applied to the particle swarm, and a damped dynamic relaxation algorithm is used to iteratively solve the momentum conservation equation to simulate the settlement and consolidation process of the slope under its own weight; the maximum velocity or maximum unbalanced force of the particle swarm is monitored, and when it is less than the preset gravity equilibrium convergence threshold... When the slope reaches its initial static equilibrium state, the stress tensor of each particle at this time is saved as the initial effective stress tensor for subsequent calculations.

[0056] Specifically, the gravity equilibrium convergence threshold is used to determine whether the initial geostress field of a slope under gravity-only conditions has reached a stable state. In SPH calculations, particles undergo acceleration and motion under gravity until they come to rest under resistance. This threshold defines the critical condition for particle motion to cease. Those skilled in the art can select this threshold according to different simulation scales. Generally, a maximum particle velocity threshold of 10 is chosen. −5 m / s to 10 −8 m / s, or take the maximum unbalanced force ratio as 10 −5 If the threshold is too large, the initial stress field will not be fully balanced, resulting in false deformations in the early stages of subsequent freeze-thaw calculations; if the threshold is too small, it will significantly increase the calculation time in the initialization stage and reduce efficiency.

[0057] Step S2: Calculate the physical field evolution index in the particle neighborhood based on the current state variables of the particle swarm, and compare the physical field evolution index with a preset threshold to obtain the region identification result.

[0058] Specifically, step S2 includes the following process:

[0059] Based on the temperature state variable carried by the particle, the temperature gradient modulus in the particle's neighborhood is calculated, and the region where the temperature gradient modulus is greater than a first threshold is identified as an active phase transition region.

[0060] Specifically, the process of calculating the temperature gradient modulus within the particle's neighborhood includes:

[0061] The SPH kernel function approximation method is adopted. Based on the temperature, mass, density and kernel function gradient vector of the particle and its neighboring particles, the temperature gradient tensor is calculated by discrete summation. The norm of the tensor is taken as the temperature gradient modulus.

[0062] Specifically, the first threshold This refers to the temperature gradient modulus threshold, which defines the boundary of the region in soil where a dramatic phase transition (water-ice transition) occurs. Near a freezing front, the temperature gradient is significantly higher than in the frozen or unfrozen areas. By setting this threshold, the location of the phase transition interface can be accurately captured, thereby triggering particle density operations to precisely simulate the growth of ice lenses. The value of this threshold is related to the soil type and freezing rate. Generally, The range of values ​​is If the value is too large, it may result in a narrow identification range of the phase transition region, missing key details of frost heave; if the value is too small, it may result in an overly wide identification range, causing unnecessary waste of computing resources.

[0063] Specifically, a particle neighborhood search mechanism is constructed. Based on the smoothness length *h* of the current particle, its influence radius is determined. All particles within the computational domain are traversed, and neighboring particles within the influence radius of the current particle are selected to construct a neighborhood particle list. Active phase transition regions are identified. Based on the temperature state variable carried by the particle, the temperature gradient modulus within the particle's neighborhood is calculated using the SPH kernel function approximation method. The specific calculation formula is as follows:

[0064]

[0065]

[0066] Where i is the current particle, and j is a neighboring particle. Let m be the set of neighboring particles, ρ be the density, and T be the temperature. The kernel function gradient; the calculated temperature gradient modulus Second threshold Perform a comparison, if If so, the particle is marked as a particle in the phase transition active region.

[0067] Based on the effective stress tensor carried by the particle, the equivalent plastic strain increment of the particle is calculated, and the region where the particle with the equivalent plastic strain increment is greater than the second threshold is identified as the shear damage zone.

[0068] Specifically, the second threshold This threshold represents the equivalent plastic strain increment, used to define the area where significant plastic shear failure occurs in the soil. During slope instability, the plastic strain increment accumulates significantly at the slip zone. This threshold filters out minor numerical disturbances and elastic deformations, focusing only on the actual damage evolution area. This threshold is typically set to a minimum value to capture the damage initiation. Generally, The value range is 1.0 × 10 −5 ~1.0×10 -3 If the value is too large, it may cause a lag in shear zone identification, making it impossible to capture the microscopic failure characteristics of the slope initiation stage.

[0069] Specifically, based on the effective stress tensor carried by the particle, the plastic strain increment tensor at the current moment is calculated according to the plastic flow law. Further calculation of the equivalent plastic strain increment The calculation formula is:

[0070]

[0071] in, The components of the plastic strain increment tensor are represented using Einstein's summation convention.

[0072] The calculated equivalent plastic strain increment is compared with the second threshold. If... If so, the particle is marked as a particle in the shear damage zone.

[0073] The regions containing particles that do not belong to the active phase transition region and the shear damage region are identified as stable inert regions.

[0074] Step S3: Based on the region identification results and environmental boundary conditions, the heat conduction equation and moisture migration equation are solved using the SPH discrete scheme that considers the latent heat source of phase change, so as to obtain the temperature field update data, phase change interface location data and moisture field distribution data at the current moment.

[0075] Specifically, step S3 includes the following process:

[0076] Construct environmental boundary entity particles, and assign temperature time history or water potential state to the environmental boundary entity particles according to the environmental boundary conditions.

[0077] Specifically, environmental boundary entity particles are constructed. A layer of virtual environmental boundary entity particles is generated outside the boundary of the computational domain, and these particles are assigned a time-varying temperature time-history function based on environmental monitoring data. Or a constant water potential state It is used to simulate the exchange of heat and moisture between the atmospheric environment and the slope surface.

[0078] Based on the region identification results, the thermal conductivity and permeability coefficients of particles in different regions are dynamically assigned corresponding values.

[0079] Specifically, for particles in the phase change active region, the equivalent thermal conductivity and equivalent heat capacity are dynamically calculated and assigned based on the current temperature and unfrozen water content function, and the permeability coefficient in the transition region is assigned; for particles in the shear damage region, the thermal conductivity after damage and the corrected permeability coefficient are assigned; for particles in the stable inert region, the fixed thermal conductivity and initial permeability coefficient based on the initial physical state are assigned.

[0080] The heat flux and moisture flux between a particle and its neighboring particles and the boundary particles of the environment are calculated using the discrete scheme of the SPH kernel function. The temperature field update data and moisture field distribution data at the current moment are obtained by solving the problem.

[0081] Specifically, the process of calculating the heat flux and moisture flux between a particle and its neighboring particles and environmental boundary particles based on the SPH kernel function discrete scheme includes:

[0082] In the heat conduction calculation, a discrete equation is constructed based on the temperature difference between the calculated particle and its neighboring particles, thermal conductivity, mass, density, and kernel function gradient. A latent heat source term for phase change is introduced to solve for the rate of temperature change. The thermal conductivity is dynamically assigned based on the particle position and region identification results.

[0083] Specifically, in SPH calculations, when the distance between adjacent particles approaches zero, the product of the kernel function gradient and the reciprocal of the distance becomes numerically singular (the denominator is zero). The anti-singularity coefficient is a very small positive number used to smooth the kernel function gradient calculation and prevent numerical divergence. It is typically set between 0.1h and 0.2h (where h is the smoothing length). This value effectively eliminates singularities without significantly interfering with the physical reality of the inter-particle interactions.

[0084] Specifically, the latent heat source term of phase change represents the enormous heat absorbed or released during the phase change of water-ice, and is the core element that distinguishes freeze-thaw simulation from ordinary heat conduction calculations. By introducing the equivalent heat capacity method or the sensible heat capacity method, the latent heat term is integrated into the energy equation, enabling temperature field calculations to accurately capture the latent heat release effect at the phase change interface. When calculating the rate of temperature change, it is necessary to monitor the rate of change of the unfrozen water content in real time. The node temperature is then corrected based on the latent heat of phase change.

[0085] Specifically, based on the calculated temperature difference between the particle and its neighboring particles Thermal conductivity ,quality ,density The discrete equations for particle temperature change rate are constructed using the kernel function gradient, and a latent heat source term for phase change is introduced to handle the phase change process. The resulting discrete equations for heat conduction SPH are as follows:

[0086]

[0087] in, Let be the equivalent heat capacity of particle i. The latent heat of the ice-water phase transition. Ice content, To prevent singular coefficients;

[0088] If the neighboring particles contain environmental boundary entity particles, then replace j with the boundary particle index and use the environmental temperature. It participates in the calculation to apply the heat flux boundary conditions.

[0089] In the calculation of water migration, a discrete equation is constructed based on the pore pressure difference between the calculated particle and its neighboring particles, the permeability coefficient, and the kernel function gradient. An anti-singularity coefficient is introduced to correct the calculation singularity when the particle spacing is too close, and the water flux is obtained by solving the equation.

[0090] Specifically, based on the calculated pore pressure difference between the particle and its neighboring particles... Permeability coefficient The discrete equations for water migration (SPH) are constructed using kernel function gradients and anti-singularity coefficients are introduced to correct computational singularities when particle spacing is too close. The constructed SPH discrete equations are as follows:

[0091]

[0092] in, Where is the unfrozen water content, and p is the pore water pressure. The water source and sink terms caused by the phase change are determined by updating the unfrozen water content and ice content to determine the number of phase change interface locations.

[0093] Specifically, the thermal conductivity of frozen soil changes significantly with variations in ice content and unfrozen water content. In regions with active phase transitions, dynamic interpolation formulas (such as the weighted average method) are used to calculate transient thermal conductivity to reflect the abrupt changes in the thermophysical properties of the soil during phase transitions and improve simulation accuracy.

[0094] Step S4: Based on the temperature field update data, the variable interface position data, and the region identification results, an adaptive adjustment strategy is executed on the particles of the particle swarm, and the freeze-thaw constitutive model is combined to calculate the freeze-thaw strain increment and stress response to obtain the current particle resolution, current stress tensor, and current plastic state.

[0095] Specifically, the process of implementing an adaptive adjustment strategy for the particles in the particle swarm includes:

[0096] For particles identified as active phase transition regions by the region identification results, when their temperature gradient modulus exceeds the first threshold, a particle splitting operation is performed to decompose a single parent particle into a preset number of sub-particles. Based on the principles of mass conservation and momentum conservation, the sub-particles inherit the mass, velocity, temperature, and effective stress tensor of the parent particle. At the same time, the smooth length of the sub-particles is reduced accordingly to improve the local computational resolution.

[0097] Specifically, for particles identified as being in the active phase transition region, their temperature gradient modulus is calculated. With the first threshold The difference ratio; when Significantly greater than Furthermore, if the current time step is greater than the preset minimum time step, a splitting decision is triggered; the individual parent particle is decomposed into... Individual particles, in two-dimensional computation In three-dimensional calculation The position offset of the child particle relative to the parent particle is set as Based on the principle of mass conservation, the mass of the subparticle Set as the mass of the parent particle of Based on the principle of momentum conservation, the daughter particle inherits the velocity vector of the parent particle; the temperature, unfrozen water content, and effective stress tensor of the daughter particle directly inherit the values ​​of the parent particle; the smooth length of the daughter particle... Reduced to the smooth length of the parent particle of To match the new particle spacing and improve local computational resolution.

[0098] For particles identified as stable inert regions by the region identification results, when the particle spacing is less than a preset merging threshold and the difference in physical field variables is lower than a limit value, a particle merging operation is performed to aggregate neighboring particles and reconstruct them into a single coarse particle.

[0099] Specifically, in particle merging determination, the variable difference threshold ensures that the merged coarse particles accurately represent the state of the microscopic physical field before merging, preventing the forced merging of particles with significant temperature or stress differences, which could lead to the loss of local physical field characteristics (such as temperature peaks or stress concentration points). For the temperature difference threshold, a value of 0.5 is typically used. ∘ C∼1.0 ∘ C; The limit value for the difference in stress tensor norm is usually taken as 10 kPa to 50 kPa.

[0100] Specifically, for a particle identified as a stable inert region, search for neighboring inert region particles within a radius of 2h centered on that particle; when the number of neighboring particles reaches a preset merging quantity... (Usually consistent with the number of splits), and the differences in physical field variables between adjacent particles (including temperature differences and stress tensor norm differences) are all less than the preset variable difference limit. At that time, perform the merge operation; Neighboring particles aggregate and reconstruct into a single coarse particle; the mass of the coarse particle is... The sum of the masses of all merged particles is used; the position and velocity of the coarse particles are calculated based on a mass-weighted average to ensure the conservation of the total momentum and center of mass position of the system before and after merging; the smooth length of the coarse particles is increased to match the coarsened particle spacing and reduce local computational overhead.

[0101] Specifically, the process of calculating the freeze-thaw strain increment and stress response in conjunction with the dynamic freeze-thaw constitutive model to obtain the current particle resolution, current stress tensor, and current plastic state includes:

[0102] Based on the temperature change and phase transition interface location data in the temperature field update data, combined with the changes in soil porosity and unfrozen water content, the frost heave strain increment is calculated using the frost heave coefficient.

[0103] Specifically, based on the updated temperature field data at the current moment, the temperature drop ΔT and the change in unfrozen water content within the current time step are calculated. Introducing the coefficient of frost heave The volume expansion caused by in-situ freezing and segregation freezing is calculated to obtain the frost heave strain increment. :

[0104]

[0105] Where n is the porosity, For the increase in ice content, The value is typically taken from 0.01 to 0.05, depending on the soil type.

[0106] Specifically, the frost heave coefficient is a coefficient that characterizes the ability of soil to expand in volume during freezing, and it is related to the soil particle composition, pore structure, and external moisture supply conditions. This coefficient is a key coupling parameter connecting temperature field changes and mechanical deformation fields.

[0107] The total strain increment is calculated based on particle motion data, and the thermal strain increment and the frost heave strain increment are subtracted from the total strain increment to obtain the effective mechanical strain increment.

[0108] Specifically, the formula for calculating the effective mechanical strain increment is:

[0109]

[0110] in, For effective mechanical strain increment, This represents the total strain increment. is the coefficient of thermal expansion, and I is the unit tensor.

[0111] The effective mechanical strain increment is calculated using the elastic constitutive matrix to obtain the test stress tensor, and the cohesion and internal friction angle parameters in the yield criterion are dynamically updated based on the current temperature.

[0112] Specifically, the calculation process for the test stress tensor is as follows:

[0113]

[0114] Specifically, the yield criterion parameters are dynamically updated based on the current temperature T, and the cohesion at the current moment is calculated using a temperature dependence function. and friction angle .

[0115] The radial back mapping algorithm is used to plastically correct the experimental stress tensor, and the stress state is mapped back to the updated yield surface to obtain the current stress tensor and the current plastic state.

[0116] Specifically, a yield function F(σ,c(T),ϕ(T)) based on the updated parameters is constructed. If the test stress exceeds the yield surface, the plasticity correction factor Δλ is calculated to pull the stress state back to the yield surface along the steepest descent direction, thus obtaining the current stress tensor σ_new and the current plastic state (including cumulative equivalent plastic strain).

[0117] Step S5: Solve the momentum conservation equation to obtain particle motion data based on the current stress tensor and the current plastic state, and correct the soil permeability coefficient based on the current plastic state.

[0118] Specifically, step S5 includes the following process:

[0119] Based on the current stress tensor, the stress divergence term in the particle's neighborhood is calculated using the SPH kernel function approximation method. The momentum conservation equation is then constructed by combining the gravitational volume force term, and the particle acceleration is obtained by solving it.

[0120] Specifically, the stress divergence term is derived in a symmetrical form to ensure momentum conservation, and the acceleration formula for particle i is as follows:

[0121]

[0122] in, Let be the acceleration of particle i. , where are the current stress components of the particle and its neighboring particles, respectively, and g is the gravitational acceleration vector. This is an artificial viscosity term, used to dissipate high-frequency numerical oscillations and prevent particles from penetrating non-physically.

[0123] The particle acceleration is integrated using an explicit time integration algorithm to update the particle's velocity and position data;

[0124] Specifically, the Verlet integral algorithm is used to first update the particle velocity half-step using the current acceleration, then calculate the position full-step using the updated velocity, and finally calculate the acceleration at the next moment using the force at the new position.

[0125]

[0126]

[0127]

[0128] in, The time step is represented by the superscripts n and n+1, which indicate the time step nodes.

[0129] Extract the cumulative equivalent plastic strain in the current plastic state, substitute it into the preset permeability coefficient evolution model, and dynamically correct the soil permeability coefficient.

[0130] Specifically, the process of dynamically correcting the soil permeability coefficient includes:

[0131] Damage variables are introduced to describe the evolution of soil micropore structure, and a mapping relationship between cumulative equivalent plastic strain and damage variables is established.

[0132] Specifically, the mapping relationship is expressed as: Where D is the damage variable, The damage initiation threshold is when When D=0.

[0133] The change in porosity is calculated based on the damage variables, and the permeability coefficient correction factor at the current moment is calculated using the exponential damage permeability formula.

[0134] Specifically, the exponential damage permeability formula is expressed as follows: ,in, This is the permeability sensitivity coefficient.

[0135] The initial permeability coefficient is multiplied by the permeability correction factor to obtain the corrected soil permeability coefficient, which is used for water migration calculation at the next time step.

[0136] Specifically, the initial permeability coefficient Multiplying this by the permeability correction factor yields the corrected soil permeability coefficient. : .

[0137] Step S6: Based on the particle motion data, determine the time step of the next calculation step using the Courant condition, and repeat steps S2 to S5 until the preset freeze-thaw cycle is completed.

[0138] Specifically, the process of determining the time step size of the next calculation time step using the Courant condition includes:

[0139] The critical step size for mechanical stability was calculated based on the interparticle spacing and the current sound velocity, the critical step size for thermal stability was calculated based on the thermal diffusivity, and the critical step size for water migration stability was calculated based on the corrected permeability coefficient.

[0140] Specifically, the formula for calculating the critical step size Δt_mech for mechanical stability is:

[0141]

[0142] Where h is the smooth length of the current particle. The velocity of sound in the soil (calculated from the elastic modulus and density of the soil). Let be the maximum velocity modulus in the current particle swarm; this formula ensures that the propagation distance of the stress wave within one time step does not exceed the particle spacing.

[0143] The formula for calculating the thermal stability critical step size Δt_ther is:

[0144]

[0145] in, The thermal diffusivity of the soil is . ,in, Thermal conductivity, The formula, based on specific heat capacity, limits the diffusion range of heat flow within a time step, preventing the temperature field calculation from diverging.

[0146] The formula for calculating the steady-state critical step size Δt_water for water migration is:

[0147]

[0148] in, The effective diffusion coefficient for water migration is related to the corrected permeability coefficient and the parameters of the soil-water characteristic curve. This formula ensures the stability of the explicit solution of the water field, and the step size limit is particularly critical when the permeability coefficient increases sharply due to damage.

[0149] The minimum value among the mechanical stability critical step size, thermal stability critical step size, and moisture migration stability critical step size is selected and multiplied by a preset safety factor to serve as the time step size for the next calculation step.

[0150] Specifically, the time step size for the next calculation step is expressed as follows: ,in For safety, a value of 0.2-0.8 is used. The current accumulated time is then set. Updated to .

[0151] Specifically, the safety factor represents the theoretical upper limit of the explicit computational stability provided by the Courant conditions. In practical calculations, due to the truncation error of the SPH kernel function, the uncertainty of boundary treatment, and the influence of the nonlinear constitutive model, directly using the theoretical critical step size still carries the risk of numerical divergence. The safety factor is used to further reduce the time step, ensuring the absolute stability of the computation.

[0152] Determine whether the current accumulated time has reached the preset freeze-thaw cycle period. If it has, stop the calculation and output the slope displacement field, stress field and plastic zone distribution data.

[0153] Specifically, determine the current cumulative time after the update. Has the preset freeze-thaw cycle been reached? :

[0154] when If the time condition is met, return to step S2 and use the updated particle position, velocity, temperature, stress, and permeability coefficient to calculate the physical field evolution for the next time step; when When the calculation is complete, the loop stops, and the displacement field, stress field, plastic zone distribution data, and particle resolution evolution data of the entire slope process are output.

[0155] Specifically, this invention accurately identifies active phase transition regions and shear damage regions by calculating the temperature gradient modulus and equivalent plastic strain increment. It drives the particle resolution to adaptively adjust according to the intensity of physical field changes, so that the computational node density is precisely matched with the spatial characteristics of thermodynamic gradient and mechanical damage evolution, effectively solving the contradiction between computational accuracy and efficiency under fixed resolution. At the same time, it uses cumulative plastic strain to describe pore structure damage and corrects the permeability coefficient in real time, establishing a dynamic feedback mechanism between mechanical damage and water migration, thereby realizing high-fidelity simulation of the multi-field coupled evolution process of freeze-thaw slopes.

[0156] The technical solution of the present invention has been described above with reference to the preferred embodiments shown in the accompanying drawings. However, it will be readily understood by those skilled in the art that the scope of protection of the present invention is obviously not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions after these changes or substitutions will all fall within the scope of protection of the present invention.

[0157] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for simulating the risk of freeze-thaw slopes in high-altitude and cold regions based on regional identification, characterized in that... include: Step S1: Based on the slope geological survey data, the computational domain is discretized using the smooth particle hydrodynamics method, and the particle swarm carrying mass, position, temperature, volume of unfrozen water content, porosity and effective stress tensor is obtained through gravity initialization. Step S2: Calculate the physical field evolution index in the particle neighborhood based on the current state variables of the particle swarm, and compare the physical field evolution index with a preset threshold to obtain the region identification result. Step S3: Based on the region identification results and environmental boundary conditions, the heat conduction equation and moisture migration equation are solved using the SPH discrete scheme that considers the latent heat source of phase change, so as to obtain the temperature field update data, phase change interface location data and moisture field distribution data at the current moment. Step S4: Based on the temperature field update data, the variable interface position data, and the region identification results, an adaptive adjustment strategy is executed on the particles of the particle swarm, and the freeze-thaw constitutive model is combined to calculate the freeze-thaw strain increment and stress response to obtain the current particle resolution, current stress tensor, and current plastic state. Step S5: Solve the momentum conservation equation to obtain particle motion data based on the current stress tensor and the current plastic state, and correct the soil permeability coefficient based on the current plastic state. Step S6: Based on the particle motion data, determine the time step of the next calculation step using the Courant condition, and repeat steps S2 to S5 until the preset freeze-thaw cycle is completed. The process of implementing an adaptive adjustment strategy for the particles in the particle swarm includes: For particles identified as active phase transition regions by the region identification results, when their temperature gradient modulus exceeds the first threshold, a particle splitting operation is performed to decompose a single parent particle into a preset number of sub-particles. Based on the principles of mass conservation and momentum conservation, the sub-particles inherit the mass, velocity, temperature, and effective stress tensor of the parent particle. At the same time, the smooth length of the sub-particles is reduced accordingly to improve the local computational resolution. For particles identified as stable inert regions by the region identification results, when the particle spacing is less than a preset merging threshold and the difference in physical field variables is lower than a limit value, a particle merging operation is performed to aggregate neighboring particles and reconstruct them into a single coarse particle.

2. The method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification as described in claim 1, characterized in that, The process of step S2 includes: Based on the temperature state variable carried by the particle, the temperature gradient modulus in the particle's neighborhood is calculated, and the region where the temperature gradient modulus is greater than a first threshold is identified as an active phase transition region. Based on the effective stress tensor carried by the particle, the equivalent plastic strain increment of the particle is calculated, and the region where the particle with the equivalent plastic strain increment is greater than the second threshold is identified as the shear damage zone. The regions containing particles that do not belong to the active phase transition region and the shear damage region are identified as stable inert regions.

3. The method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification according to claim 2, characterized in that, The process of calculating the temperature gradient modulus in the particle's neighborhood includes: The SPH kernel function approximation method is adopted. Based on the temperature, mass, density and kernel function gradient vector of the particle and its neighboring particles, the temperature gradient tensor is calculated by discrete summation. The norm of the tensor is taken as the temperature gradient modulus.

4. The method for simulating the risk of high-altitude, cold-climate freeze-thaw slopes based on region identification according to claim 3, characterized in that, The process of step S3 includes: Construct environmental boundary entity particles, and assign temperature time history or water potential state to the environmental boundary entity particles according to the environmental boundary conditions. Based on the region identification results, the thermal conductivity and permeability coefficients of particles in different regions are dynamically assigned corresponding values. The heat flux and moisture flux between a particle and its neighboring particles and the boundary particles of the environment are calculated using the SPH kernel function discrete scheme. The current temperature field update data and moisture field distribution data are obtained by solving the problem.

5. The method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification according to claim 4, characterized in that, The process of calculating the heat flux and water flux between a particle and its neighboring particles and environmental boundary particles based on the SPH kernel function discrete scheme includes: In the heat conduction calculation, a discrete equation is constructed based on the temperature difference between the calculated particle and its neighboring particles, thermal conductivity, mass, density, and kernel function gradient. A latent heat source term for phase change is introduced to solve for the rate of temperature change. The thermal conductivity is dynamically assigned based on the particle position and region identification results. In the calculation of water migration, a discrete equation is constructed based on the pore pressure difference between the calculated particle and its neighboring particles, the permeability coefficient, and the kernel function gradient. An anti-singularity coefficient is introduced to correct the calculation singularity when the particle spacing is too close, and the water flux is obtained by solving the equation.

6. The method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification as described in claim 5, is characterized in that... The process of calculating the freeze-thaw strain increment and stress response in conjunction with the dynamic freeze-thaw constitutive model to obtain the current particle resolution, current stress tensor, and current plastic state includes: Based on the temperature change and phase transition interface location data in the temperature field update data, combined with the changes in soil porosity and unfrozen water content, the frost heave strain increment is calculated using the frost heave coefficient. The total strain increment is calculated based on particle motion data, and the thermal strain increment and the frost heave strain increment are subtracted from the total strain increment to obtain the effective mechanical strain increment. The effective mechanical strain increment is calculated using the elastic constitutive matrix to obtain the test stress tensor, and the cohesion and internal friction angle parameters in the yield criterion are dynamically updated based on the current temperature. The radial back mapping algorithm is used to plastically correct the experimental stress tensor, and the stress state is mapped back to the updated yield surface to obtain the current stress tensor and the current plastic state.

7. The method for simulating the risk of freeze-thaw slopes in high-altitude and cold regions based on regional identification as described in claim 6, is characterized in that... The process of step S5 includes: Based on the current stress tensor, the stress divergence term in the particle's neighborhood is calculated using the SPH kernel function approximation method. The momentum conservation equation is then constructed by combining the gravitational volume force term, and the particle acceleration is obtained by solving it. The particle acceleration is integrated using an explicit time integration algorithm to update the particle's velocity and position data; Extract the cumulative equivalent plastic strain in the current plastic state, substitute it into the preset permeability coefficient evolution model, and dynamically correct the soil permeability coefficient.

8. The method for simulating the risk of freeze-thaw slopes in high-altitude and cold regions based on regional identification according to claim 7, characterized in that, The process of dynamically correcting the soil permeability coefficient includes: Damage variables are introduced to describe the evolution of soil micropore structure, and a mapping relationship between cumulative equivalent plastic strain and damage variables is established. The change in porosity is calculated based on the damage variables, and the permeability coefficient correction factor at the current moment is calculated using the exponential damage permeability formula. The initial permeability coefficient is multiplied by the permeability correction factor to obtain the corrected soil permeability coefficient, which is used for water migration calculation at the next time step.

9. The method for simulating the risk of freeze-thaw slopes in cold and high-altitude regions based on regional identification as described in claim 8, characterized in that, The process of determining the time step size of the next calculation time step using the Courant condition includes: The critical step size for mechanical stability was calculated based on the interparticle spacing and the current sound velocity, the critical step size for thermal stability was calculated based on the thermal diffusivity, and the critical step size for water migration stability was calculated based on the corrected permeability coefficient. The minimum value among the mechanical stability critical step size, thermal stability critical step size, and moisture migration stability critical step size is selected and multiplied by a preset safety factor to serve as the time step size for the next calculation step. Determine whether the current accumulated time has reached the preset freeze-thaw cycle period. If it has, stop the calculation and output the slope displacement field, stress field and plastic zone distribution data.

Images