A Multi-Field Coupled Analysis Method for Loess Slopes Considering the Influence of Porosity

Abstract

This invention discloses a multi-field coupled analysis method for loess slopes considering the influence of porosity, comprising: geometric modeling using finite element software; establishing a constitutive model related to porosity; introducing a porosity random field and atmospheric boundary conditions; calculating the slope safety factor using an improved strength reduction method; and obtaining the safety factor. The advantage of this invention is that the improved strength reduction method can be better used for calculating the slope safety factor and better simulates the stability of loess slopes under hydraulic coupling conditions.

Description

Technical Field

[0001] This invention relates to the fields of mechanics and civil engineering, and in particular to a multi-field coupling analysis method for loess slopes that takes into account the influence of porosity. Background Technology

[0002] Loess slopes possess unique engineering characteristics, and the heterogeneity of their porosity is one of the key factors affecting slope stability. Loess is a special type of soil with a typical large-pore structure, formed during geological history through soil formation and weathering. Due to variations in sedimentary environment and organic matter content, the porosity distribution of loess often exhibits significant heterogeneity.

[0003] The distribution of heterogeneous porosity in loess slopes is mainly influenced by factors such as sedimentary environment, soil formation, weathering processes, and human activities. This heterogeneity results in significant spatial differences in the permeability, shear strength, and deformation characteristics of loess slopes. In high-porosity areas, the soil is more prone to water absorption, leading to soil softening and reduced strength, thus increasing the risk of landslides. Conversely, in low-porosity areas, the soil is relatively dense and has higher shear strength, but it may still fail under saturated conditions due to increased pore water pressure.

[0004] Extensive research has been conducted both domestically and internationally on the stability of loess slopes, employing methods such as laboratory experiments, field observations, analytical solutions, and numerical simulations to analyze the behavior of unsaturated soils or slope stability. Laboratory experiments are helpful in understanding the hydraulic coupling characteristics of soil, but they suffer from limitations such as high time consumption and cost. While analytical solutions provide a fast and effective method for analyzing slope stability, they still have limitations when dealing with more complex engineering projects. Numerical simulation technology also plays an important role in studying loess slope stability. Numerical simulation methods establish detailed numerical models to simulate the impact of external factors such as rainfall on slope stability and evaluate the effectiveness of different protective measures. Various constitutive models are used to analyze rainfall-induced landslides, such as the improved Drucker-Prager model, the Barcelona Basic Model (BBM), or the dual-porosity constitutive model. However, current constitutive models do not effectively couple porosity with shear strength parameters, especially for hydraulic coupling scenarios, and do not couple porosity with the inherent permeability coefficient. Currently, regarding the stability of heterogeneous slopes, some researchers have introduced porosity stochastic fields to characterize the stability of heterogeneous slopes, but they have not considered the influence of porosity on shear strength. Other researchers have proposed cohesion stochastic fields to characterize the heterogeneity of slopes, studying the influence of anisotropic depositional directions on slope stability and failure mechanisms; however, these stochastic fields do not consider the impact of porosity on cohesion. Furthermore, although some researchers have considered the porosity of loess slopes, improved the constitutive equations of loess, and studied the creep effect of homogeneous loess slopes, they have not considered the heterogeneity of porosity in loess slopes. Regarding the determination of the safety factor, some studies use abrupt displacement changes or computational non-convergence as the criterion for strength reduction methods, but the rate of displacement change sometimes does not exhibit an increasing pattern, making it difficult to use as a basis.

[0005] In conclusion, studying the distribution patterns of heterogeneous porosity on loess slopes is of great significance for accurately predicting slope stability and designing effective protective measures. This not only helps prevent landslide disasters but also provides a scientific basis for engineering construction in loess regions.

[0006] Current constitutive models do not effectively couple porosity with shear strength parameters, especially in the case of hydraulic coupling. They do not couple porosity with the inherent permeability coefficient, nor do they couple porosity with both shear strength parameters and inherent permeability coefficient simultaneously.

[0007] Current strength reduction methods for calculating safety factors rely on sudden displacement changes as the criterion. However, the rate of displacement change sometimes does not exhibit an increasing pattern, making it difficult to use as a reliable basis. This is an area that this application needs to focus on improving. Summary of the Invention

[0008] The technical problem to be solved by this invention is to provide a multi-field coupling analysis method for loess slopes that takes into account the influence of porosity. The improved strength reduction method can be better used for slope safety factor calculation and better simulate the stability of loess slopes under hydraulic coupling conditions.

[0009] To address the above technical problems, this invention provides a multi-field coupling analysis method for loess slopes that considers the influence of porosity, comprising the following steps:

[0010] Step S1: Geometric modeling using finite element software;

[0011] Step S2: Establish a constitutive model related to porosity;

[0012] 1) Solid phase;

[0013] For the mechanical analysis of heterogeneous loess slopes, the total strain rate considered is divided into two parts: elastic (dε) e / dt) and plasticity (dε) vp / dt);

[0014]

[0015] For the elastic part, D e It is the elastic stiffness matrix, E is the elastic modulus, and v is Poisson's ratio;

[0016] For the plastic part, F m Γ is the yield surface, G is the plastic potential, Γ is the fluidity, and m is the stress power.

[0017] The elastic constitutive model employs a standard linear elastic model with Young's modulus and Poisson's ratio, and the elastic modulus E varies with porosity:

[0018]

[0019] In the formula: E0 is the reference elastic modulus, φ is the porosity, and φ 0E Reference porosity;

[0020] The plastic constitutive model uses a modified Mohr-Coulomb model:

[0021]

[0022] In the formula: ψ is the friction angle, ω is the shear dilatation angle, p′, q, and θ are the effective mean stress, deviatoric stress, and Rhodes angle, respectively.

[0023] Meanwhile, the effect of porosity on cohesion is also considered:

[0024] c = (a c +b cs)g(φ);

[0025] g(φ) = (f + |f|) / 2;

[0026] f(φ) = 1 - (φ / φ) 0c ) n ;

[0027] In the formula: c is the cohesive force, a c 、b c φ is the cohesive force parameter. 0c Reference porosity;

[0028] 2) Liquid phase and gas phase;

[0029] The Van-Genuchten model is used as the water-holding curve to describe the suction variation of loess, and its specific expression is as follows:

[0030]

[0031] In the formula: S e This represents the actual moisture content; S rl Residual moisture content, i.e., the minimum moisture content in the soil; S ls The maximum water content is equal to the water content when the soil is saturated; S l The current moisture content; P g P is the pore pressure. l Let be the pore water pressure, and the difference between λ and ρ be the matrix suction; λ is the shape control parameter of the water holding curve, which is only related to the corresponding soil; P is the air inlet value of the soil, expressed as follows:

[0032]

[0033] P0(φ)=P0Exp{a(φ-φ 0S )};

[0034] λ(φ)=λExp{b(φ 0s )};

[0035] In the formula: P0 is the intake air value corresponding to temperature T0; σ is the surface tension; φ is the porosity. 0s The reference porosity is shown; a and b are relevant parameters.

[0036] Darcy's law is used to characterize the liquid and gas phases, and its specific expression is:

[0037]

[0038] Where: k is the inherent permeability; k rα It is relative permeability; μ α It's viscosity; It is the pore water pressure gradient between two points; ρ α ρ is the density of the liquid; g is the gravitational acceleration vector.

[0039] Using the exponential law of inherent permeability:

[0040] k=k0Exp{b k (φ-φ 0k )};

[0041] Where: φ 0k The reference porosity is k0, which is the corresponding intrinsic permeability.

[0042] The relative permeability of the liquid and gas phases is defined as follows:

[0043]

[0044] k rg =1-k rl ;

[0045] The constitutive model used takes into account the effect of porosity and couples porosity with elastic modulus, cohesion, and intrinsic permeability coefficient.

[0046] Step S3: Introduce a porosity random field and atmospheric boundary conditions;

[0047] Step S31: Introduce a porosity random field to characterize the heterogeneous properties of loess slopes;

[0048] The generation of the spatially correlated heterogeneous porosity field is based on Davis's research and follows the semivariogram relationship of geostatistical theory:

[0049]

[0050] Where C0 is the nugget value; C0+C1 is the base value; a is the associated range; |h| is the distance between two mutually influencing elements;

[0051] This is further transformed into a spatial problem. Under the condition of second-order stationarity, we obtain a spatially constant mean and variance. The relationship between the covariance function and the semivariogram function is C(h) = C(0) - γ(h), which gives:

[0052]

[0053] Generate a porosity random field with average porosity μ and covariance C0+C1. The formulas for its mean and covariance are as follows:

[0054] E(φ)=μ

[0055] Cov(φ)=C0+C1

[0056] The vector h is modified as follows to effectively characterize the anisotropy of space:

[0057]

[0058] Where, a x a y and a z These are the corresponding ranges of 'a' on the x, y, and z axes, respectively; x ani y ani z ani These are the original coordinates; α, β, and γ are the angles by which the original coordinate axes x, y, and z are rotated to the current coordinate axes, respectively.

[0059] Step S32: Introduce atmospheric boundary conditions;

[0060] Atmospheric boundary conditions are introduced, and boundary conditions are applied to evaporation and precipitation to simulate complex soil-atmosphere interactions; flux boundary conditions are used to represent water, air and energy, while the mass change functions of the three components are determined by three variables: liquid pressure, gas pressure and soil temperature;

[0061] The equation for calculating gas flux at the atmospheric interface is:

[0062] q g -γ g (P g -P at );

[0063] In the formula: γ g P is the gas infiltration coefficient. at At atmospheric pressure, this flux consists of two parts: dry air and water vapor, denoted as... Let the water vapor mass fraction be the dry air flux, then:

[0064]

[0065] Calculating water flux requires considering multiple factors, and its overall expression is as follows:

[0066]

[0067] In the formula: the water flux between the Earth's surface and the external environment is the rainfall intensity I, and the surface evaporation is E. v Steam flow rate Surface runoff j sr The sum of four parts; where the coefficient k rain and K evap This is a coefficient for rainfall and evaporation, used to adjust the mass ratio of externally input water; evaporation is expressed by the following formula:

[0068]

[0069] Where veg is defined as the vegetation fraction of green area per unit ground area, r a It is aerodynamic drag. It is the stability coefficient, k is the von Kármán constant, and v a It is wind speed, z0 is ground roughness related to canopy height, z a It measures wind speed and atmospheric absolute humidity (ρ). va The grid height, ρ v It refers to the absolute moisture content of the soil;

[0070] Where, ρ va and ρ v These represent the absolute humidity of the atmosphere and the atmosphere at the boundary condition nodes, respectively.

[0071] The horizontal flow rate of water is shown in the following formula:

[0072]

[0073] In the formula: ρ at Atmospheric gas density;

[0074] When actual rainfall is applied, considering that the infiltrated portion of rainfall is the difference between actual rainfall and the portion of rainfall intercepted outside the system, this portion is considered runoff:

[0075]

[0076] In the formula: γ w This is the water infiltration coefficient. Adjusting this value simulates infiltration conditions in different scenarios or with different types of soil and rock. The formula shows that the amount of runoff depends on the difference between pore water pressure and atmospheric pressure at the surface. The closer the soil is to saturation, the more difficult it is for water to infiltrate, which is consistent with reality. Finally, if the soil is saturated (P... l >P at If all the non-permeable rainfall becomes runoff, it means that this rainfall does not enter the system.

[0077] Step S4: Improved strength reduction method;

[0078] The slope safety factor is determined using the strength reduction method, with the reduction factor F. r The definition is the current actual strength parameter value. With reduced strength The ratio and the relationship between relevant parameters:

[0079]

[0080] In the formula: c and These are the cohesion and friction angle of the soil and rock materials, respectively.

[0081] Step S5: Obtain the safety factor;

[0082] Using the penetration of the plastic shear zone as a criterion for slope instability, the safety factor of loess slopes is calculated, including the following steps:

[0083] Step S51: Run precipitation model A according to atmospheric boundary conditions to obtain water pressure, and fix it according to the grid onto another set of precipitation models B with no atmospheric boundary conditions but otherwise the same conditions.

[0084] Step S52: Precipitation model B generates two time periods. The first time period is not reduced, and the material parameters are the initial cohesion c and friction angle. The second time period inputs the reduction factor and reduces it according to the strength reduction method formula, where c=(a+bs)g(φ), which reduces the cohesive parameters a and b for cohesive force reduction;

[0085] Step S53: For the calculation results, determine whether the plastic deviatoric strain at the top edge is greater than the set ε. If "yes", proceed to step S54; if "no", return to step S52, and simultaneously input the reduction coefficient of the second time period +0.01.

[0086] Step S54: The reduction factor for the second time period is -0.01;

[0087] Step S55: Determine whether the plastic deviatoric strain at the top edge is less than the set ε. If yes, proceed to step S56; if no, return to step S54.

[0088] Step S56: Safety factor = reduction factor at this time + 0.01, output safety factor.

[0089] When the plastic shear zone penetrates, a landslide occurs on the slope, which is consistent with the actual situation.

[0090] The superior effects of this invention are as follows:

[0091] 1) Effectively consider the influence of the non-homogeneity of porosity on loess slopes, and introduce a porosity random field to characterize the non-homogeneity of loess slopes;

[0092] 2) The improved strength reduction method can be better used to calculate the slope safety factor and better simulate the stability of loess slopes under hydraulic coupling conditions. Attached Figure Description

[0093] The accompanying drawings, which form part of this application, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0094] Figure 1This is a flowchart of a specific embodiment of the present invention;

[0095] Figure 2 This is a flowchart illustrating the calculation of the safety factor for loess slopes in a specific embodiment of the present invention;

[0096] Figure 3 This is a schematic diagram showing the dimensions and fixed boundaries of models A and B in a specific embodiment of the present invention;

[0097] Figure 4 This is a schematic diagram of porosity distribution in a specific embodiment of the present invention;

[0098] Figure 5 This is a schematic diagram of the penetration rate distribution in a specific embodiment of the present invention;

[0099] Figure 6 This is a schematic diagram of pore water pressure distribution in a specific embodiment of the present invention;

[0100] Figure 7 This is a schematic diagram of plastic deviatoric strain (reduction factor of 3.03) in a specific embodiment of the present invention. Detailed Implementation

[0101] The embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0102] In this specific embodiment of the invention, the slope angle of the loess slope is assumed to be θ (tangent of slope angle tanθ = 1 / 2), and the slope height is 10 meters. The loess porosity exhibits a heterogeneous distribution, with a mean porosity of 0.5 and a variance of 0.002, forming a random strip-like distribution. The rainfall intensity is 10 mm / day, and the slope stability is analyzed after 5 days of rainfall.

[0103] like Figure 1 As shown, this invention provides a multi-field coupling analysis method for loess slopes that considers the influence of porosity, including the following steps:

[0104] Step S1: Geometric modeling using finite element software;

[0105] Geometric modeling was performed using the finite element software CODE_BRIGHT as required, such as... Figure 3 As shown, the left and right boundaries x restrict its x-axis direction, and the lower boundary restricts its x and y directions. Two identical models A and B are established.

[0106] Step S2: Establish a constitutive model related to porosity;

[0107] Based on loess material, constitutive parameters were designed. In the solid phase: elastic component: elastic modulus E is 100, reference porosity φ 0E =0.5; Plastic part: cohesion a c It is 0.25, b c φ is 0.015.0c The internal friction angle is 0.7, the exponent n is 5, and the internal friction angle is 0.7. The angle is 15°. In the liquid phase: water holding curve, P0 is 0.02, λ is 0.18, σ is 0.072, a is 2, b is 2, φ 0S The intrinsic permeability is 0.7; k0 is 1.023 x 10⁻⁶. -12 , φ 0k It is 0.5;

[0108] 1) Solid phase;

[0109] For the mechanical analysis of heterogeneous loess slopes, the total strain rate considered is divided into two parts: elastic (dε) e / dt) and plasticity (dε) vp / dt);

[0110]

[0111] For the elastic part, D e It is the elastic stiffness matrix, E is the elastic modulus, and v is Poisson's ratio;

[0112] For the plastic part, F m Γ is the yield surface, G is the plastic potential, Γ is the fluidity, and m is the stress power.

[0113] The elastic constitutive model employs a standard linear elastic model with Young's modulus and Poisson's ratio, and the elastic modulus E varies with porosity:

[0114]

[0115] In the formula: E0 is the reference elastic modulus, φ is the porosity, and φ 0E Reference porosity;

[0116] The plastic constitutive model uses a modified Mohr-Coulomb model:

[0117]

[0118] In the formula: ψ is the friction angle, ω is the shear dilatation angle, p′, q, and θ are the effective mean stress, deviatoric stress, and Rhodes angle, respectively.

[0119] Meanwhile, the effect of porosity on cohesion is also considered:

[0120] c = (a c +b c s)g(φ);

[0121] g(φ) = (f + |f|) / 2;

[0122] f(φ) = 1 - (φ / φ)0c ) n ;

[0123] In the formula: c is the cohesive force, a c 、b c φ is the cohesive force parameter. 0c Reference porosity;

[0124] 2) Liquid phase and gas phase;

[0125] The Van-Genuchten model is used as the water-holding curve to describe the suction variation of loess, and its specific expression is as follows:

[0126]

[0127] In the formula: S e This represents the actual moisture content; S rl Residual moisture content, i.e., the minimum moisture content in the soil; S ls The maximum water content is equal to the water content when the soil is saturated; S l The current moisture content; P g P is the pore pressure. l Let be the pore water pressure, and the difference between λ and ρ be the matrix suction; λ is the shape control parameter of the water holding curve, which is only related to the corresponding soil; P is the air inlet value of the soil, expressed as follows:

[0128]

[0129] P0(φ)=P0Exp{a(φ-φ 0S )};

[0130] λ(φ)=λExp{b(φ 0s )};

[0131] In the formula: P0 is the intake air value corresponding to temperature T0; σ is the surface tension; φ is the porosity. 0S The reference porosity is shown; a and b are relevant parameters.

[0132] Darcy's law is used to characterize the liquid and gas phases, and its specific expression is:

[0133]

[0134] Where: k is the inherent permeability; k rα It is relative permeability; μ α It's viscosity; It is the pore water pressure gradient between two points; ρ α ρ is the density of the liquid; g is the gravitational acceleration vector.

[0135] Using the exponential law of inherent permeability:

[0136] k=k0Exp{b k (φ-φ 0k )};

[0137] Where: φ 0k The reference porosity is k0, which is the corresponding intrinsic permeability.

[0138] The relative permeability of the liquid and gas phases is defined as follows:

[0139]

[0140] k rg =1-k rl ;

[0141] The constitutive model used takes into account the effect of porosity and couples porosity with elastic modulus, cohesion, and intrinsic permeability coefficient.

[0142] Step S3: Introduce a porosity random field and atmospheric boundary conditions;

[0143] For the porosity random field, the generation has a mean of 0.5, a variance of 0.002, and is horizontally oriented, randomly distributed in stripes, fixed on a grid, such as... Figure 4 As shown. Due to the influence of the constitutive model, its permeability has changed with porosity, as... Figure 5 As shown;

[0144] For model A, atmospheric boundary conditions are applied, with a rainfall intensity of 10 mm / day, 5 days of rainfall, and a temperature of 20℃. The resulting pore water pressure distribution map is shown below. Figure 6 As shown;

[0145] Step S31: Introduce a porosity random field to characterize the heterogeneous properties of loess slopes;

[0146] The generation of the spatially correlated heterogeneous porosity field is based on Davis's research and follows the semivariogram relationship of geostatistical theory:

[0147]

[0148] Where C0 is the nugget value; C0+C1 is the base value; a is the associated range; |h| is the distance between two mutually influencing elements;

[0149] This is further transformed into a spatial problem. Under the condition of second-order stationarity, we obtain a spatially constant mean and variance. The relationship between the covariance function and the semivariogram function is C(h) = C(0) - γ(h), which gives:

[0150]

[0151] Generate a porosity random field with average porosity μ and covariance C0+C1. The formulas for its mean and covariance are as follows:

[0152] E(φ)=μ

[0153] Cov(φ)=C0+C1

[0154] The vector h is modified as follows to effectively characterize the anisotropy of space:

[0155]

[0156] Where, a x a y and a z These are the corresponding ranges of 'a' on the x, y, and z axes, respectively; x ani y ani z ani These are the original coordinates; α, β, and γ are the angles by which the original coordinate axes x, y, and z are rotated to the current coordinate axes, respectively.

[0157] Step S32: Introduce atmospheric boundary conditions;

[0158] Atmospheric boundary conditions are introduced, and boundary conditions are applied to evaporation and precipitation to simulate complex soil-atmosphere interactions; flux boundary conditions are used to represent water, air and energy, while the mass change functions of the three components are determined by three variables: liquid pressure, gas pressure and soil temperature;

[0159] The equation for calculating gas flux at the atmospheric interface is:

[0160] q g =γ g (P g -P at );

[0161] In the formula: γ g P is the gas infiltration coefficient. at At atmospheric pressure, this flux consists of two parts: dry air and water vapor, denoted as... Let the water vapor mass fraction be the dry air flux, then:

[0162]

[0163] Calculating water flux requires considering multiple factors, and its overall expression is as follows:

[0164]

[0165] In the formula: the water flux between the Earth's surface and the external environment is the rainfall intensity I, and the surface evaporation is E. v Steam flow rate Surface runoff j sr The sum of four parts; where the coefficient k rain and k evap This is a coefficient for rainfall and evaporation, used to adjust the mass ratio of externally input water; evaporation is expressed by the following formula:

[0166]

[0167] Where veg is defined as the vegetation fraction of green area per unit ground area, r a It is aerodynamic drag. It is the stability coefficient, k is the von Kármán constant, and v a It is wind speed, z0 is ground roughness related to canopy height, z a It measures wind speed and atmospheric absolute humidity (ρ). va The grid height, ρ v It refers to the absolute moisture content of the soil;

[0168] Where, ρ va and ρ v These represent the absolute humidity of the atmosphere and the atmosphere at the boundary condition nodes, respectively.

[0169] The horizontal flow rate of water is shown in the following formula:

[0170]

[0171] In the formula: ρ at Atmospheric gas density;

[0172] When actual rainfall is applied, considering that the infiltrated portion of rainfall is the difference between actual rainfall and the portion of rainfall intercepted outside the system, this portion is considered runoff:

[0173]

[0174] In the formula: γ w This is the water infiltration coefficient. Adjusting this value simulates infiltration conditions in different scenarios or with different types of soil and rock. The formula shows that the amount of runoff depends on the difference between pore water pressure and atmospheric pressure at the surface. The closer the soil is to saturation, the more difficult it is for water to infiltrate, which is consistent with reality. Finally, if the soil is saturated (P... l >P at If all the non-permeable rainfall becomes runoff, it means that this rainfall does not enter the system.

[0175] Step S4: Improved strength reduction method;

[0176] Fix the water pressure of model A onto model B according to... Figure 2Strength reduction is applied, with ε set to 0.02. When the reduction factor is 3.03, the maximum value of the upper boundary plastic deviatoric strain is 0.0239, and when the reduction factor is 3.02, the maximum value of the upper boundary plastic deviatoric strain is 0.015.

[0177] The slope safety factor is determined using the strength reduction method, with the reduction factor F. r The definition is the current actual strength parameter value. With reduced strength The ratio and the relationship between relevant parameters:

[0178]

[0179] In the formula: c and These are the cohesion and friction angle of the soil and rock materials, respectively.

[0180] Step S5: Obtain the safety factor.

[0181] like Figure 7 As shown, the safety factor is 3.03, which means that for loess slopes with this porosity distribution, the safety factor is 3.03 under the condition of 10 mm / day rainfall for 5 days.

[0182] When the plastic shear band penetrates, a landslide occurs on the slope, which is consistent with actual conditions. Considering that finite element numerical simulation software has certain errors due to the influence of mesh accuracy, a small amount of plastic deviatoric strain value is used as the judgment criterion, rather than a zero value.

[0183] 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 multi-field coupling analysis method for loess slopes considering the influence of porosity, comprising the following steps: Step S1: Geometric modeling using finite element software; Step S2: Establish a constitutive model related to porosity; Establish solid, liquid, and gas phases, and couple porosity with elastic modulus, cohesion, and intrinsic permeability coefficient; including the following: 1) Solid phase; Mechanical analysis of heterogeneous loess slopes, considering the total strain rate. It is divided into two parts: elasticity ( ) and plasticity ( ); ； For the elastic part, It is the elastic stiffness matrix. It is the elastic modulus. It is Poisson's ratio; For the plastic part, It is the yield surface. It is a plastic potential. It's about liquidity. It is stress power. It is the effective stress; The elastic constitutive model employs a standard linear elastic model with Young's modulus and Poisson's ratio, and the elastic modulus... As porosity changes: ； In the formula: For reference elastic modulus, Porosity The elastic modulus is referenced to the porosity. The plastic constitutive model uses a modified Mohr-Coulomb model: ； In the formula: It is the angle of friction. It is the shear expansion angle. It is a plastic potential parameter. , , These are the effective mean stress, deviatoric stress, and Rhodes angle, respectively. Meanwhile, the effect of porosity on cohesion is also considered: ； ； ； In the formula: For cohesion, , For cohesion parameters, The porosity is used as a reference for cohesion. 2) Liquid phase and gas phase; The Van-Genuchten model is used as the water-holding curve to describe the suction variation of loess, and its specific expression is as follows: ； In the formula: This represents the actual moisture content. Residual moisture content, which is the minimum moisture content in the soil. The maximum moisture content is equal to the moisture content when the soil is saturated. This represents the current moisture content. Pore ​​pressure, The difference between the pore water pressure and the matrix suction is the matrix suction. The shape control parameter of the water-holding curve is only related to the corresponding soil. This is the air intake value of the soil, expressed as follows: ； ； ； In the formula: for The corresponding intake air value at the temperature; Surface tension; Porosity Reference porosity; , For relevant parameters; Darcy's law is used to characterize liquid flux. Gas flux The specific expression is: ； in: It is the inherent penetration rate; It is the relative penetration rate; It's viscosity; It is the pore water pressure gradient between two points; The density of the liquid; It is the vector of gravitational acceleration; At this time, all parameters represent liquid phase parameters. At that time, all parameters represent gas phase parameters; Using the exponential law of inherent permeability: ; in: Porosity is used as a reference for permeability. This corresponds to the inherent penetration rate; Relative permeability of the liquid phase Relative permeability of the gas phase The definition is as follows: ； ； The constitutive model used takes into account the effect of porosity and couples porosity with elastic modulus, cohesion, and intrinsic permeability coefficient. Step S3: Introduce a porosity random field and atmospheric boundary conditions; Step S31: Introduce a porosity random field to characterize the heterogeneous properties of loess slopes; The spatially correlated heterogeneous porosity field is generated, and follows the semivariogram function relationship of geostatistical theory: ； in, Value of a nugget; This is the base value; For the related range; The distance between two mutually influencing elements; This is further transformed into a spatial problem. Under the condition of second-order stationarity, we obtain a spatially constant mean and variance. The relationship between its covariance function and semivariogram function is as follows: ,get: ； Generate an average porosity of covariance is The porosity random field has the following formulas for its mean and covariance: ； ； vector The following changes are made to effectively characterize spatial anisotropy: ； ； ； ； ； in, , and They are Corresponding to the respective ranges on the x, y, and z axes; , , It is a vector The coordinates; , , These are the original coordinates; , , These are the angles by which the original coordinate axes x, y, and z are rotated to the current coordinate axes; Step S32: Introduce atmospheric boundary conditions to simulate complex soil-atmosphere interactions; Step S4: Calculate the slope safety factor using the improved strength reduction method; The slope safety factor is determined using the strength reduction method, and the reduction factor is... The definition is the current actual strength parameter value ( , ) and strength after reduction ( , The ratio of ) and the relationship between relevant parameters: ； In the formula: and These are the cohesion and internal friction angle of the soil and rock materials, respectively. Using the penetration of plastic shear bands as a criterion for slope instability, the safety factor of loess slopes is calculated. Step S5: Obtain the safety factor.

2. A computer system comprising a memory, a processor, and a computer program stored in the memory, characterized in that: The processor executes the computer program to implement the steps of the multi-field coupling analysis method of claim 1.

Images