Loess slope multi-field coupling analysis method considering porosity influence

Abstract

The invention discloses a loess slope multi-field coupling analysis method considering the porosity influence. The loess slope multi-field coupling analysis method comprises the steps that finite element software geometric modeling is carried out; establishing a constitutive model related to the porosity; introducing a porosity random field and atmospheric boundary conditions; solving a slope safety coefficient through an improved strength reduction method; and the safety coefficient is obtained. The method has the advantages that the improved strength reduction method can be better used for slope safety coefficient calculation, and the stability of the loess slope under the hydraulic coupling condition can be better simulated.

Description

Technical Field

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

[0002] Loess slopes have unique engineering properties, 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, which is formed by soil formation and weathering during geological history. Due to changes in the sedimentary environment and the different organic matter content in the soil, the porosity distribution of loess often shows significant heterogeneity.

[0003] The distribution of heterogeneous porosity in loess slopes is mainly affected by factors such as sedimentary environment, soil formation, weathering process and human activities. This heterogeneity makes the permeability, shear strength and deformation characteristics of loess slopes show large spatial differences. In high-porosity areas, the soil is more likely to absorb water, resulting in soil softening and reduced strength, thereby increasing the risk of landslides. Conversely, in low-porosity areas, the soil is denser and has higher shear strength, but it may still be damaged due to increased pore water pressure in a saturated state.

[0004] A large number of studies have been conducted at home and abroad on the stability of loess slopes, using indoor experiments, field observations, analytical solutions or numerical simulations to analyze the behavior of unsaturated soil or slope stability. Indoor experiments help to understand the hydraulic coupling characteristics of soil, but indoor experiments have limitations such as long time consumption and high cost. Although analytical solutions provide a fast and effective method for analyzing slope stability, they still have limitations when it comes to more complex projects. Numerical simulation technology also plays an important role in studying the stability of loess slopes. Numerical simulation methods simulate the influence of external factors such as rainfall on slope stability by establishing a sophisticated numerical model, and evaluate the effects of different protective measures. Various constitutive models are used by this method to analyze rainfall-induced landslides, such as the improved Drucker-Prager model, the Barcelona Basic Model (BBM) or the dual-porosity constitutive model. However, the existing constitutive models do not effectively couple porosity with shear strength parameters, especially for hydraulic coupling, and do not couple porosity with intrinsic permeability. At present, in order to solve the problem of stability of heterogeneous slopes, some researchers have introduced porosity random fields to characterize the stability of heterogeneous slopes, but have not considered the influence of porosity on shear strength. Some researchers have also proposed random fields of cohesion to characterize the heterogeneity of slopes, and studied the influence of anisotropic deposition direction on slope stability and failure mechanism, but this random field does not take into account the influence of porosity on cohesion. In addition, although some researchers have considered the porosity problem of loess slopes, improved the constitutive equation of loess, and studied the creep effect of homogeneous loess slopes, they have not considered the heterogeneity of porosity of loess slopes. For the determination of safety factor, some studies have used displacement mutation or non-convergence of calculation as the basis for judging the strength reduction method, but the displacement change rate sometimes does not show an increasing law, which makes it difficult to use it as a basis.

[0005] In conclusion, studying the distribution law of heterogeneous porosity in 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 areas.

[0006] The existing constitutive models do not effectively couple the porosity with the shear strength parameters, especially for the hydraulic coupling case, the porosity is not coupled with the intrinsic permeability, let alone the porosity is coupled with the shear strength parameters and the intrinsic permeability at the same time.

[0007] The current strength reduction method for solving the safety factor uses displacement mutation as the basis for judgment, but the displacement change rate sometimes does not show an increasing law, which makes it difficult to use it as a basis. This is the area that needs to be improved in this application. Summary of the invention

[0008] The technical problem to be solved by the present invention is to provide a multi-field coupling analysis method for loess slopes taking 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] In order to solve the above technical problems, the present invention provides a multi-field coupling analysis method for loess slope considering the influence of porosity, comprising the following steps:

[0010] Step S1: Finite element software geometric modeling;

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

[0012] 1) Solid phase;

[0013] For the mechanical analysis of heterogeneous loess slope, 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 is the elastic stiffness matrix, E is the elastic modulus, and v is the 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 uses a standard linear elastic model with Young's modulus and Poisson's ratio, and the elastic modulus E varies with the porosity:

[0018]

[0019] Where: E0 is the reference elastic modulus, φ is the porosity, φ 0E is the reference porosity;

[0020] The plastic constitutive model adopts the improved Mohr-Coulomb model:

[0021]

[0022] Where: is the friction angle, ψ is the dilatancy angle, ω is the plastic potential parameter, p′, q, θ are the effective mean stress, deviatoric stress and Rhodes angle, respectively;

[0023] At the same time, consider the effect of porosity on cohesion:

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

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

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

[0027] Where: c is the cohesive force, a c , b c is the cohesion parameter, φ 0c is the 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 change of loess. Its specific expression is:

[0030]

[0031] Where: S e is the actual moisture content; S rl is the residual moisture content, i.e. the minimum value of the moisture content in the soil; S ls is the maximum moisture content, which is equal to the moisture content when the soil is saturated; S l is the current moisture content; P g is the pore pressure, P l is the pore water pressure, and the difference between the two is the matrix suction; λ is the shape control parameter of the water retention curve, which is only related to the corresponding soil; P is the air intake value of the soil, which is expressed as follows:

[0032]

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

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

[0035] Where: P0 is the corresponding air intake value at T0 temperature; σ is the surface tension; φ is the porosity, φ 0s is the reference porosity; a and b are related parameters;

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

[0037]

[0038] Where: k is the intrinsic permeability; k rα is the relative permeability; μ α is viscosity; 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 intrinsic permeability:

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

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

[0042] The relative permeabilities of the liquid and gas phases are defined as follows:

[0043]

[0044] k rg =1-k rl ;

[0045] The constitutive model used considers the effect of porosity and couples porosity with elastic modulus, cohesion, and intrinsic permeability;

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

[0047] Step S31: introducing porosity random field to characterize the heterogeneous characteristics of loess slope;

[0048] The generation of spatially correlated heterogeneous porosity fields is based on the work of Davis and follows the semivariogram relationship of geostatistical theory:

[0049]

[0050] Among them, C0 is the nugget value; C0+C1 is the base value; a is the associated range; |h| is the distance between two elements that affect each other;

[0051] Further transformed into a spatial problem, under the condition of second-order stationarity, the spatially constant mean and variance are obtained. The relationship between its covariance function and semivariogram is C(h)=C(0)-γ(h), and we get:

[0052]

[0053] Generate a porosity random field with a mean porosity of μ and a covariance of C0+C1. The mean and covariance formulas are:

[0054] E(φ)=μ

[0055] Cov(φ)=C0+C1

[0056] The vector h undergoes the following changes to effectively characterize the spatial anisotropy:

[0057]

[0058] Among them, a x 、a y and a z are the corresponding ranges of a on the x, y and z axes respectively; ani ,y ani 、z ani is the original coordinate; α, β, γ are the angles of rotation of the original coordinate axes x, y and z to the current coordinate axes;

[0059] Step S32: introducing atmospheric boundary conditions;

[0060] Atmospheric boundary conditions are introduced to impose boundary conditions on evaporation and rainfall to simulate complex soil-atmosphere interactions; they are represented by flux boundary conditions representing 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 calculation equation for the gas flux at the atmospheric interface is:

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

[0063] Where: γ g is the gas infiltration coefficient, P at is atmospheric pressure. This flux consists of dry air and water vapor. is the water vapor mass fraction, the dry air flux is expressed as:

[0064]

[0065] Calculating the water flux requires considering several effects, and the overall expression is as follows:

[0066]

[0067] Where: the water flux between the surface and the outside world is the rainfall intensity I, the ground evaporation E v , steam advection flux Surface runoff sr The sum of the four parts; the coefficient k rain and K evap It is the coefficient of rainfall and evaporation, which adjusts the mass ratio of external input water; evaporation is expressed as follows:

[0068]

[0069] where veg is the vegetation fraction defined as the green area per unit ground area, r a is the aerodynamic drag, is the stability coefficient, k is the von Karman constant, v a is the wind speed, z0 is the ground roughness associated with the canopy height, and z a It is used to measure wind speed and absolute humidity (ρ va ) grid height, ρ v is the absolute moisture of the soil;

[0070] Among them, ρ va and ρ v are the absolute humidity of the atmosphere and the atmosphere at the boundary condition nodes, respectively;

[0071] The advective flux of water is given by the following formula:

[0072]

[0073] Where: at is the atmospheric gas density;

[0074] In the case of actual rainfall impact, the part of rainfall that is considered to be infiltrated is the difference between the actual rainfall and the part of rainfall that is intercepted outside the system, which is the runoff:

[0075]

[0076] Where: γ w is the water infiltration coefficient. Adjust this value to simulate the infiltration of different scenarios or different types of rock and soil. From the formula, it can be concluded that the size of the runoff depends on the difference between the pore water pressure at the surface and the atmospheric pressure. The closer the soil is to saturation, the more difficult it is for water to infiltrate. This is also in line with the actual situation. Finally, if the soil is saturated (P l >P at ), then all rainfall that cannot infiltrate will become runoff, which means that this rainfall does not enter the system;

[0077] Step S4: improved strength reduction method;

[0078] The strength reduction method is used to calculate the slope safety factor, and the reduction factor F r is defined as the current actual strength parameter value Reduced strength The ratio of related parameters:

[0079]

[0080] Where: c and are the cohesion and friction angle of geomaterials, respectively;

[0081] Step S5: Obtaining the safety factor;

[0082] Taking the penetration of plastic shear band as the criterion for slope instability, the safety factor of loess slope is calculated, including the following steps:

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

[0084] Step S52: The precipitation model B generates two time periods. The first time period is not reduced, and the material parameters are the initial cohesion c, friction angle In the second time period, input the reduction factor and perform the reduction according to the strength reduction method, where c = (a + bs) g (φ), and the cohesion parameters a and b are reduced for the cohesion.

[0085] Step S53: for the calculation result, determine whether the plastic deviatoric strain of the top edge is greater than the set ε, "yes" to go to step S54; "no" to return to step S52, and at the same time input the reduction coefficient of the second time period +0.01;

[0086] Step S54: the reduction coefficient of the second time period is -0.01;

[0087] Step S55: determine whether the plastic deviatoric strain of 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 = current reduction factor + 0.01, output the safety factor.

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

[0090] The superior effects of the present invention are:

[0091] 1) Effectively consider the influence of porosity heterogeneity of loess slopes and introduce porosity random field to characterize the heterogeneous characteristics of loess slopes;

[0092] 2) 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. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] The drawings constituting a part of the present application are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0094] Figure 1is a flow chart of a specific embodiment of the present invention;

[0095] Figure 2 It is a flow chart of calculating the safety factor of loess slope in a specific embodiment of the present invention;

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

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

[0098] Figure 5 is a schematic diagram of permeability distribution of a specific embodiment of the present invention;

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

[0100] Figure 7 It is a schematic diagram of the plastic deviatoric strain (reduction factor is 3.03) of a specific embodiment of the present invention. DETAILED DESCRIPTION

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

[0102] The specific embodiment of the present invention assumes that the slope angle of the loess slope is θ (slope angle tangent value tanθ=1 / 2), and the slope height is 10 meters. The porosity of the loess is non-homogeneous, with a mean porosity of 0.5 and a variance of 0.002, and is randomly distributed in strips. The rainfall intensity is 10 mm / day, and the slope stability after 5 days of rainfall is analyzed.

[0103] like Figure 1 As shown, the present invention provides a multi-field coupling analysis method for loess slope considering the influence of porosity, comprising the following steps:

[0104] Step S1: Finite element software geometric modeling;

[0105] Finite element software CODE_BRIGHT is used to perform geometric modeling as required, such as Figure 3 As shown, the left and right boundaries x limit its x-axis direction, and the lower boundary limits its x, y directions, and two identical models A and B are established.

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

[0107] According to the loess material, the constitutive parameters are designed. In the solid phase: elastic part: elastic modulus E is 100, reference porosity φ 0E 0.5; Plastic part: cohesion a c is 0.25, b c is 0.015,φ0c is 0.7, index n is 5, internal friction angle In the liquid phase: water retention curve, P0 is 0.02, λ is 0.18, σ is 0.072, a is 2, b is 2, φ 0S is 0.7; inherent permeability: k0 is 1.023x10 -12 ,φ 0k is 0.5;

[0108] 1) Solid phase;

[0109] For the mechanical analysis of heterogeneous loess slope, 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 is the elastic stiffness matrix, E is the elastic modulus, and v is the 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 uses a standard linear elastic model with Young's modulus and Poisson's ratio, and the elastic modulus E varies with the porosity:

[0114]

[0115] Where: E0 is the reference elastic modulus, φ is the porosity, φ 0E is the reference porosity;

[0116] The plastic constitutive model adopts the improved Mohr-Coulomb model:

[0117]

[0118] Where: is the friction angle, ψ is the dilatancy angle, ω is the plastic potential parameter, p′, q, θ are the effective mean stress, deviatoric stress and Rhodes angle, respectively;

[0119] At the same time, consider the effect of porosity on cohesion:

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

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

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

[0123] Where: c is the cohesive force, a c , b c is the cohesion parameter, φ 0c is the 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 change of loess. Its specific expression is:

[0126]

[0127] Where: S e is the actual moisture content; S rl is the residual moisture content, i.e. the minimum value of the moisture content in the soil; S ls is the maximum moisture content, which is equal to the moisture content when the soil is saturated; S l is the current moisture content; P g is the pore pressure, P l is the pore water pressure, and the difference between the two is the matrix suction; λ is the shape control parameter of the water retention curve, which is only related to the corresponding soil; P is the air intake value of the soil, which is expressed as follows:

[0128]

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

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

[0131] Where: P0 is the corresponding air intake value at T0 temperature; σ is the surface tension; φ is the porosity, φ 0S is the reference porosity; a and b are related parameters;

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

[0133]

[0134] Where: k is the intrinsic permeability; k rα is the relative permeability; μ α is viscosity; 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 intrinsic permeability:

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

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

[0138] The relative permeabilities of the liquid and gas phases are defined as follows:

[0139]

[0140] k rg =1-k rl ;

[0141] The constitutive model used considers the effect of porosity and couples porosity with elastic modulus, cohesion, and intrinsic permeability;

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

[0143] For the porosity random field, the generated mean is 0.5, the variance is 0.002, the direction is horizontal, and it is randomly distributed in strips and fixed on the grid, such as Figure 4 As shown in Figure 2, due to the influence of the constitutive model, its permeability has changed with the porosity, as shown in Figure 2. Figure 5 As shown;

[0144] For model A, atmospheric boundary conditions are applied, the rainfall intensity is 10 mm / day, the rainfall lasts for 5 days, and the temperature is 20°C. The pore water pressure distribution diagram is obtained, as shown in Figure 6 As shown;

[0145] Step S31: introducing porosity random field to characterize the heterogeneous characteristics of loess slope;

[0146] The generation of spatially correlated heterogeneous porosity fields is based on the work of Davis and follows the semivariogram relationship of geostatistical theory:

[0147]

[0148] Among them, C0 is the nugget value; C0+C1 is the base value; a is the associated range; |h| is the distance between two elements that affect each other;

[0149] Further transformed into a spatial problem, under the condition of second-order stationarity, the spatially constant mean and variance are obtained. The relationship between its covariance function and semivariogram is C(h)=C(0)-γ(h), and we get:

[0150]

[0151] Generate a porosity random field with a mean porosity of μ and a covariance of C0+C1. The mean and covariance formulas are:

[0152] E(φ)=μ

[0153] Cov(φ)=C0+C1

[0154] The vector h undergoes the following changes to effectively characterize the spatial anisotropy:

[0155]

[0156] Among them, a x 、a y and a z are the corresponding ranges of a on the x, y and z axes respectively; ani ,y ani 、z ani is the original coordinate; α, β, γ are the angles of rotation of the original coordinate axes x, y and z to the current coordinate axes;

[0157] Step S32: introducing atmospheric boundary conditions;

[0158] Atmospheric boundary conditions are introduced to impose boundary conditions on evaporation and rainfall to simulate complex soil-atmosphere interactions; they are represented by flux boundary conditions representing 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 calculation equation for the gas flux at the atmospheric interface is:

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

[0161] Where: γ g is the gas infiltration coefficient, P at is atmospheric pressure. This flux consists of dry air and water vapor. is the water vapor mass fraction, the dry air flux is expressed as:

[0162]

[0163] Calculating the water flux requires considering several effects, and the overall expression is as follows:

[0164]

[0165] Where: the water flux between the surface and the outside world is the rainfall intensity I, the ground evaporation E v , steam advection flux Surface runoff sr The sum of the four parts; the coefficient k rain and k evap It is the coefficient of rainfall and evaporation, which adjusts the mass ratio of external input water; evaporation is expressed as follows:

[0166]

[0167] where veg is the vegetation fraction defined as the green area per unit ground area, r a is the aerodynamic drag, is the stability coefficient, k is the von Karman constant, v a is the wind speed, z0 is the ground roughness associated with the canopy height, and z a It is used to measure wind speed and absolute humidity (ρ va ) grid height, ρ v is the absolute moisture of the soil;

[0168] Among them, ρ va and ρ v are the absolute humidity of the atmosphere and the atmosphere at the boundary condition nodes, respectively;

[0169] The advective flux of water is given by the following formula:

[0170]

[0171] Where: at is the atmospheric gas density;

[0172] In the case of actual rainfall impact, the part of rainfall that is considered to be infiltrated is the difference between the actual rainfall and the part of rainfall that is intercepted outside the system, which is the runoff:

[0173]

[0174] Where: γ w is the water infiltration coefficient. Adjust this value to simulate the infiltration of different scenarios or different types of rock and soil. From the formula, it can be concluded that the size of the runoff depends on the difference between the pore water pressure at the surface and the atmospheric pressure. The closer the soil is to saturation, the more difficult it is for water to infiltrate. This is also in line with the actual situation. Finally, if the soil is saturated (P l >P at ), then all rainfall that cannot infiltrate will become runoff, which means that this rainfall does not enter the system;

[0175] Step S4: improved strength reduction method;

[0176] Fix the water pressure of model A on model B according to Figure 2Strength reduction is performed, where ε is taken as 0.02. When the reduction factor is 3.03, the maximum value of the upper boundary plastic deviator strain is 0.0239. When the reduction factor is 3.02, the maximum value of the upper boundary plastic deviator strain is 0.015.

[0177] The strength reduction method is used to calculate the slope safety factor, and the reduction factor F r is defined as the current actual strength parameter value Reduced strength The ratio of related parameters:

[0178]

[0179] Where: c and are the cohesion and friction angle of geomaterials, respectively;

[0180] Step S5: Obtain the safety factor.

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

[0182] When the plastic shear zone penetrates, the slope will slide, which is consistent with the actual situation. Considering that the finite element numerical simulation software has certain errors due to the influence of grid accuracy, the micro plastic deviator strain value is used as the judgment basis instead of 0 value.

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

Claims

1. A multi-field coupling analysis method for loess slope considering the influence of porosity, comprising the following steps: Step S1: Finite element software geometric modeling; Step S2: Establishing a constitutive model related to porosity; Establish the solid phase, liquid phase and gas phase, and couple the porosity with the elastic modulus, cohesion and intrinsic permeability; Step S3: introducing porosity random field and atmospheric boundary conditions; Step S31: introducing porosity random field to characterize the heterogeneous characteristics of loess slope; Step S32: introducing atmospheric boundary conditions to simulate complex soil-atmosphere interactions; Step S4: Calculate the slope safety factor using the improved strength reduction method; The strength reduction method is used to calculate the slope safety factor, and the reduction factor F t is defined as the current actual strength parameter value Reduced strength The ratio of related parameters: Where: C and are the cohesion and friction angle of geomaterials, respectively; The plastic shear band penetration is used as the criterion for slope instability and the safety factor of loess slope is calculated. Step S5: Obtain the safety factor.

2. The multi-field coupling analysis method for loess slope considering the influence of porosity according to claim 1 is characterized in that: Step S31: introducing porosity random field to characterize the heterogeneous characteristics of loess slope; The generation of spatially correlated heterogeneous porosity fields is based on the work of Davis and follows the semivariogram relationship of geostatistical theory: Among them, C0 is the nugget value; C0+C1 is the base value; a is the associated range; |h| is the distance between two elements that affect each other; Further transformed into a spatial problem, under the condition of second-order stationarity, the spatially constant mean and variance are obtained. The relationship between its covariance function and semivariogram is C(h)=C(0)-γ(h), and we get: Generate a porosity random field with a mean porosity of μ and a covariance of C0+C1. The mean and covariance formulas are: E(φ)=μ Cov(φ)=C0+C1 The vector h undergoes the following changes to effectively characterize the spatial anisotropy: Among them, a x 、a y and a z are the corresponding ranges of a on the axes x, y and z; ani ,y ani , Z ani are the original coordinates; α, β, γ are the angles of rotation of the original coordinate axes x, y and z to the current coordinate axes.

3. The multi-field coupling analysis method for loess slope considering the influence of porosity according to claim 1 is characterized in that: The step S2 establishes a constitutive model related to porosity, including the following: 1) Solid phase; For the mechanical analysis of heterogeneous loess slope, the total strain rate considered is divided into two parts: elastic (dε e / dt) and plasticity (dε vp / dt); For the elastic part, D e is the elastic stiffness matrix, E is the elastic modulus, and v is the Poisson’s ratio; For the plastic part, F m is the yield surface, G is the plastic potential, Γ is the fluidity, and m is the stress power; The elastic constitutive model uses a standard linear elastic model with Young's modulus and Poisson's ratio, and the elastic modulus E varies with the porosity: Where: E0 is the reference elastic modulus, φ is the porosity, φ 0E is the reference porosity; The plastic constitutive model adopts the improved Mohr-Coulomb model: Where: is the friction angle, ψ is the dilatancy angle, ω is the plastic potential parameter, p′, q, θ are the effective mean stress, deviatoric stress and Rhodes angle, respectively; At the same time, consider the effect of porosity on cohesion: c=(a c +b c s)g(φ); g(φ)=(f+|f|) / 2; f(φ)=1-(φ / φ 0c ) n ; Where: c is the cohesive force, a c 、b c is the cohesion parameter, φ 0c is the reference porosity; 2) Liquid phase and gas phase; The Van-Genuchten model is used as the water holding curve to describe the suction change of loess. Its specific expression is: Where: S e is the actual moisture content; S rl is the residual moisture content, i.e. the minimum value of the moisture content in the soil; S ls is the maximum moisture content, which is equal to the moisture content when the soil is saturated; S1 is the current moisture content; P g is the pore air pressure, P1 is the pore water pressure, and the difference between the two is the matrix suction; λ is the shape control parameter of the water retention curve, which is only related to the corresponding soil; P is the air intake value of the soil, which is expressed as follows: P0(φ)=P0Exp{a(φ-φ0φ 0S )}; λ(φ)=λExp{b(φ-φ 0S )}; Where: P0 is the corresponding air intake value at T0 temperature; σ is the surface tension; φ is the porosity, φ 0S is the reference porosity; a and b are related parameters; Darcy's law is used to characterize the liquid phase and the gas phase, and the specific expression is: Where: k is the intrinsic permeability; k rα is the relative permeability; μ α is viscosity; is the pore water pressure gradient between two points; ρ α is the density of the liquid; g is the gravitational acceleration vector; Using the exponential law of intrinsic permeability: k=k0Exp{bk(φ-φ 0k )}; Where: φ 0k is the reference porosity, K0 is the corresponding intrinsic permeability; The relative permeabilities of the liquid and gas phases are defined as follows: k rg =1-k rl ; The constitutive model used takes into account the effect of porosity and couples porosity with elastic modulus, cohesion, and intrinsic permeability.

4. 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 according to any one of claims 1 to 3.

Images