Form function KL expansion random field discretization method considering irregular domain of slope soil body

Abstract

The invention provides a shape function KL expansion random field discretization method considering a slope soil body irregular domain, and belongs to the technical field of soil body parameter random field analysis. Comprising the following steps: establishing a slope numerical model and dividing a discrete grid and a form function grid; establishing a probability distribution model of a soil parameter random field; calculating an element stiffness submatrix ke and an element mass matrix Ee; calculating an autocorrelation matrix rho and a global stiffness matrix B; assembling a global matrix; calculating a characteristic value lambda k, a characteristic function phi k and an expansion term number M; and random field discretization of the slope numerical model is realized. The random field discrete technology provided by the invention is used for efficiently and accurately analyzing the random distribution characteristics of the soil parameters in the rock-soil structure. The discrete random field model can be efficiently constructed in a complex geological structure, the calculation cost is remarkably reduced, the method can be conveniently combined with numerical simulation tools such as a finite element method, and the requirement of actual engineering for efficiency is met while the precision is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of geotechnical parameter random field analysis, and specifically to a KL expansion random field discretization method for shape functions considering irregular domains of slope soils. Background Technique

[0002] In the field of geotechnical engineering, the spatial variability of soil parameters is crucial for the safety of engineering design and construction. These parameters include, but are not limited to, aspects such as foundation settlement, slope stability, and tunnel excavation. To accurately describe and analyze this spatial variability, the Karhunen-Loève (KL) expansion method is widely used to discretize the random field into a linear combination of a set of orthogonal random variables. This method significantly reduces the computational complexity by simplifying the mathematical representation and analysis of the random field, and preserves the key statistical properties of the random field, thereby improving the efficiency and accuracy of simulating the spatial variability of soil parameters.

[0003] Reference 1 "Tan Xiaohui, Dong Xiaole, Fei Suozhu, Gong Wenping, Xiu Lintian, Hou Xiaoliang, Ma Haichun, Reliability Analysis Method Based on KL Expansion and Its Application, Chinese Journal of Geotechnical Engineering, 2020, 42(05), 808-816" demonstrates the application advantages of KL expansion in the discretization of random fields and the calculation of geotechnical structure reliability. Although the KL expansion method has significant advantages in theoretical analysis, it still faces many challenges in the practical application of geotechnical random fields. First, the traditional KL expansion has difficulties in solving the second-kind Fredholm integral equation, which usually involves high computational costs. Especially when dealing with multi-dimensional data, this computational cost will further increase, limiting the practicality of the KL expansion method. To solve this problem, researchers have proposed various improvement methods. For example, in Reference 2 "Gu Xin, Zhang Wengang, Ou Qiang, Reliability Analysis of Soil Slope Stability Based on Chebyshev-Galerkin-KL Expansion [J], Chinese Journal of Geotechnical Engineering, 2023, 45(12), 2472-2480", the Galerkin method and the second-kind Chebyshev polynomial are used to solve the second-kind Fredholm integral equation. This method simplifies the solution process to a certain extent, but still requires the calculation of multiple integrals, resulting in low computational efficiency, especially when dealing with complex geotechnical engineering problems. In Reference 3 "Lin X, Tan X, Yao Y, Dong X, Fei S, Ma L, Realization of multi-dimensional random field based on Jacobi-Lagrange-Galerkin method in geotechnical engineering, Computers and Geotechnics, 2022, 144, 104533 (Lin X, Tan X H, Yao Y C, Dong X L, Fei S Z, Ma L, Realization of multi-dimensional random field based on Jacobi-Lagrange-Galerkin method in geotechnical engineering, Computers and Geotechnics, 2022, 144, 104533)", an autocorrelation matrix is established based on the interpolation method, and the multi-dimensional integral is transformed into the form of multiplying multiple one-dimensional integrals through multi-dimensional Jacobi polynomials, optimizing the computational efficiency, but still cannot completely avoid the calculation of multiple integrals, which may be a limiting factor in practical applications.

[0004] When dealing with the discretization problem of irregular random fields, the existing Karhunen-Loève (KL) expansion methods also have some limitations. The current KL expansion methods usually discretize by complementing the irregular region into a regular region, which will introduce additional discretization errors and affect the accuracy of random field discretization. For example, in Document 4 "Pranesh S, Ghosh D, Addressing the curse of dimensionality in S SFEM using the dependence of eigenvalues in KL expansion on domain size, Computer Methods in Applied Mechanics and Engineering, 2016, 311, 457-475 (Pranesh S, Ghosh D, Using the dependence of eigenvalues in KL expansion on domain size to address the curse of dimensionality in SSFEM, Computer Methods in Applied Mechanics and Engineering, 2016, 311, 457-475)", it is mathematically proven that the size of the discretized domain will significantly affect the calculation of eigenvalues in the KL expansion, indicating that the size and shape of the domain have an important impact on the accuracy of the KL expansion. Although Document 4 proposes using domain decomposition techniques to reduce the number of expansion terms in the calculation process to improve computational efficiency, this method may still not be efficient enough when dealing with large-scale or complex geotechnical engineering problems. Document 5 "Basmaji A A, Dannert M M, Nackenhorst U, Implementation of karhunen-loève expansion using discontinuous legendre polynomial based galerkin approach, Probabilistic Engineering Mechanics, 2022, 67, 103176 (Basmaji AA, Dannert M M, Nackenhorst U, Implementation of KL expansion using the Galerkin method based on discontinuous Legendre polynomials, Probabilistic Engineering Mechanics, 2022, 67, 103176)" proposes a Galerkin method based on discontinuous Legendre polynomials, which constructs Legendre basis functions on each local element domain to handle irregular random field domains. However, the applicability and efficiency of this method in multi-dimensional complex cases still need to be improved.

[0005] In summary, when dealing with irregular random fields, the existing KL expansion methods can effectively express the statistical characteristics of random fields, but there are significant limitations, including:

[0006] 1. Since the KL expansion depends on the choice of the discrete domain, when the random field has a complex geometry, traditional methods often fail to accurately capture these irregular features. This can lead to large errors during the discretization process, affecting the accuracy of the final simulation results.

[0007] 2. The computational efficiency and accuracy of the KL expansion also face challenges. Traditional KL expansion methods often improve the accuracy by integrating high-order basis functions, but this often results in a decrease in computational efficiency. At the same time, in the case of a limited number of expansion terms, it may still be unable to fully retain the main characteristics of the random field, affecting the reliability of the final calculation results. Summary of the Invention

[0008] The technical problem to be solved by the present invention is the problem existing in the above-mentioned prior art. Specifically, the present invention provides a shape function KL expansion random field discretization method considering the irregular domain of slope soil mass. This method not only considers the influence of the discrete domain on the KL expansion method, can efficiently and accurately calculate the eigenvalues and eigenfunctions, but also can discretize the irregular random field with fewer random variables.

[0009] To achieve the above object, the technical solution adopted by the present invention is:

[0010] A shape function KL expansion random field discretization method considering the irregular domain of slope soil mass, comprising the following steps:

[0011] Step 1, establish a slope numerical model and divide the discrete grid and the shape function grid

[0012] Step 1.1, establish a slope numerical model in the finite element software ABAQUS

[0013] The slope is composed of the upper slope body and the lower soil mass. The cross-section of the slope body part is trapezoidal, and the cross-section of the soil mass part is rectangular. Establish a slope model in the finite element software ABAQUS. The specific process is as follows: Set the slope model in a cross-section perpendicular to the earth's surface. The slope model includes the upper slope body part and the lower soil mass part, and set the bottom side length of the soil mass to be horizontal and the left side length to be vertical. Specifically, denote the left-bottom corner endpoint of the soil mass as endpoint A S , the width of the soil mass in the horizontal direction is K S , the height in the vertical direction is H S , the upper width of the slope body part is K F , the lower width is K M , the height in the vertical direction is H F ; Taking endpoint A S as the origin, the horizontal direction as the x direction, and the vertical direction as the y direction, establish an xy plane coordinate system;

[0014] Step 1.2, Separating the Discrete Mesh and the Shape Function Mesh

[0015] The slope numerical model is discretized to generate a discrete mesh and a shape function mesh respectively. Both meshes use quadrilateral elements. In ABAQUS, the inp files corresponding to the two meshes are exported, and the mesh node coordinate information and element node information are extracted from them. The mesh node coordinate information records the geometric positions of all nodes in the spatial coordinate system, and the element node information records the node numbers associated with each mesh element and their connection relationships.

[0016] Step 2, Calculating the Center Points of Discrete Mesh Elements and the Effective Element Indexes of the Shape Function Mesh

[0017] Let the coordinate set of the center points of the discrete mesh elements be denoted as set u * , Set u * belongs to E * ×1 space, denoted as where E * is the number of discrete mesh elements, I is the discrete mesh element number, represents the coordinate of the center point of the I-th discrete mesh element, and the superscript "T" represents the transpose of the matrix; let the coordinate set of the center points of the discrete mesh elements in the local coordinate space be denoted as set Set belongs to E * ×1 space, denoted as where, represents the local coordinate of the center point of the I-th discrete mesh element; let the set of the areas of the discrete mesh elements be denoted as set S * , Set S * belongs to E * ×1 space, denoted as where represents the area of the I-th discrete mesh element;

[0018] Let the set of the effective element indexes of the shape function mesh be denoted as set e * , Set e * belongs to E * ×1 space, denoted as where represents the number of the shape function mesh element where the center point of the I-th discrete mesh element is located;

[0019] Step 3, Establishing the Probability Distribution Model of the Random Field of Soil Parameters

[0020] Based on the slope numerical model established in Step 1, a random field probability distribution model of soil parameters is further constructed. The soil parameters involved include, but are not limited to, cohesion, internal friction angle, elastic modulus, Poisson's ratio, and density. The components of the random field probability distribution model cover the type, mean, standard deviation, and autocorrelation function of the random field probability distribution.

[0021] Let (x, y) and (x′, y′) be the coordinates of the center points of any two discrete grid cells in the slope numerical model. The random field of soil parameters is the random field H(x, y), and its probability distribution type is lognormal distribution.

[0022] Let the mean of the random field H(x, y) be μ(x, y), the standard deviation be σ(x, y), and the logarithmic mean be μ LN (x, y), μ LN (x, y) = lnμ(x, y) - [ln(1 + δ 2 (x, y))] / 2, where δ(x, y) is the coefficient of variation. Let the logarithmic standard deviation of the random field H(x, y) be σ LN (x, y), σ LN (x, y) = [ln(1 + δ 2 (x,y))] 1 / 2 ;

[0023] Let the autocorrelation function of the random field H(x, y) be the autocorrelation function ρ((x, y), (x′, y′)). This autocorrelation function ρ((x, y), (x′, y′)) is a single-exponential function, and its expression is:

[0024]

[0025] In the formula, L h is the horizontal autocorrelation distance, and L v is the vertical autocorrelation distance;

[0026] Step 4, calculate the element stiffness submatrix k e and the element mass matrix E e

[0027] The element stiffness submatrix k e and the element mass matrix E e correspond to each shape function grid cell and are both calculated in the local coordinate space. A consistent element Gaussian integration point and shape function order P are selected for each shape function grid cell, and their expressions are respectively:

[0028]

[0029] In the formula, is the weight of the element Gaussian integration point, and Je is the value of the Jacobi determinant at the unit Gaussian integration point, and diag(·) is the symbol for converting a vector into a diagonal matrix; N e is the element shape function matrix, and is the matrix composed of the values of the element shape function N at the unit Gaussian integration point ; the element shape function N belongs to the N * ×1 space, denoted as N * is the number of element shape functions;

[0030] The element stiffness submatrix k e belongs to the m * ×N * space, denoted as where m * is the number of unit Gaussian integration points; the element mass matrix E e belongs to the N * ×N * space, denoted as

[0031] Step 5, calculate the autocorrelation matrix ρ and the global stiffness matrix B

[0032] First, generate the shape function index. The shape function index is the number of all shape functions in the shape function grid, denoted as IEN, which belongs to the e s ×N * space, where e s is the number of elements in the shape function grid; the generation rule of IEN depends on the shape function order; when the shape function order P = 1, the value of IEN corresponds to the node number of the shape function grid element; when P = 2, 3, IEN is based on the previous order number and adds the number of each side of the shape function grid element; when P = 4, IEN is based on the previous order number and adds the side number and face number of the shape function grid element;

[0033] The autocorrelation matrix ρ is the matrix calculated by the autocorrelation function for all unit Gaussian integration points. Specifically, map the coordinates of all unit Gaussian integration points from the local coordinate system to the global coordinate system, and denote the set of coordinates of all unit Gaussian integration points as z * , where is the coordinate of the t-th all-unit Gaussian integration point, and M * is the number of all unit Gaussian integration points; finally, calculate the autocorrelation matrix ρ formed by all unit Gaussian integration points through the autocorrelation function;

[0034] The global stiffness matrix B belongs to the N s ×N s space, denoted as is defined by the global stiffness submatrix K and the autocorrelation matrix ρ, B = KT ρK, where N s is the total number of shape functions for the shape function grid; the global stiffness sub - matrix K is assembled from the element stiffness sub - matrices k e according to the shape function index IEN;

[0035] Step 6, calculate the eigenvalue λ k , eigenfunction φ k and the number of expansion terms M

[0036] Step 6.1, introduce the eigenvalue matrix A, where E g is the global mass matrix, assembled from the element mass matrices E e according to the shape function index IEN, and the superscript "-1" represents the inverse operation of the matrix; the eigenvalue matrix A belongs to N s ×N s space, denoted as

[0037] The eigenvalue λ k is the k - th eigenvalue in the eigenvalue vector λ obtained by sorting the eigenvalues of the eigenvalue matrix A in descending order, λ belongs to N s ×1 space, denoted as

[0038] The eigenfunction corresponding to the eigenvalue λ k is denoted as φ k , which belongs to E * ×1 space, denoted as where D u is the eigenfunction matrix, and D uk is the column vector of D u ;

[0039] Step 6.2, calculate the number of expansion terms M

[0040] Step 6.2.1, initialize the number of expansion terms M of the KL method to 1;

[0041] Step 6.2.2, denote the discretization error of the number of expansion terms M as the discretization error ε M , and its expression is:

[0042]

[0043] where, is the autocorrelation function value corresponding to the center points of the I - th and I′ - th discrete grid cells; are the eigenfunction values corresponding to the center points of the I - th and I′ - th discrete grid cells in the k - th eigenfunction φ k respectively; S is the total area of the discrete grid region; are the areas of the \(I\)-th and \(I'\)-th discrete grid cells respectively;

[0044] Step 6.2.3, compare the discrete error \(\varepsilon\) M with the given allowable discrete error \(\varepsilon\) a : If \(\varepsilon\) M \(>\varepsilon\) a , then increase the expansion term number \(M\) by 1 and return to Step 6.2.2; if \(\varepsilon\) M \(\leq\varepsilon\) a , then the iteration ends and the expansion term number \(M\) is output;

[0045] Step 7, realize the discretization of the random field of the slope numerical model

[0046] The discretization of the random field is the product of the independent standard normal random variable \(\zeta\) and the eigenvalue \(\lambda\) k , the eigenfunction \(\varphi\) k to obtain the estimated value of the random field at each point; the independent standard normal random variable \(\zeta = [\zeta_1,\zeta_2,\cdots,\zeta\) k ,\(\cdots,\zeta\) M T , representing the randomness of the geotechnical parameters, belonging to the \(M\times1\) space, denoted as \(\zeta\in R\) M×1 , where \(\zeta\) k is the random variable corresponding to the eigenvalue \(\lambda\) k , with a mean of 0 and a standard deviation of 1;

[0047] For the geotechnical parameters that conform to the normal random distribution, its mean is denoted as the normal mean \(\mu(u\) * ), the variance is denoted as the normal variance \(\sigma(u\) * ), and the discretized value of the random field is denoted as the normal discretized value \(H(u\) * ). The discretization expression of the random field of this geotechnical parameter is:

[0048]

[0049] In the formula, \(\zeta\) kn is the \(n\)-th sampling result of the \(k\)-th random variable \(\zeta\) k ;

[0050] For the geotechnical parameters that conform to the lognormal distribution, its mean is denoted as the lognormal mean \(\mu\) L (u * ), the variance is denoted as the lognormal variance \(\sigma\) L (u * ), the discretized value of the random field is denoted as the lognormal discretized value \(H\) L (u * ). The discretization expression of the random field of this geotechnical parameter is: ​

[0051]

[0052] Preferably, the implementation process of step 2 is as follows:

[0053] Calculate the central point coordinates of the I-th discrete grid cell based on the node coordinate information of the discrete grid described in step 1 And the area of the discrete grid cell And obtain the set u * And the set The central point coordinates of the discrete grid cell are the average of the node coordinates of each discrete grid cell; the area of the discrete grid cell is the area of each quadrilateral cell;

[0054] The effective element index refers to the shape function grid cell number that contains at least one central point of the discrete grid cell inside. The specific determination method is as follows:

[0055] Use the geometric determination method to evaluate the spatial relationship between the central point of each discrete grid cell and each shape function grid cell; check each central point of the discrete grid cell to determine whether it is inside the shape function grid cell; when the central point of the discrete grid cell is on the left boundary or the lower boundary of the shape function grid cell, it is determined to be inside the shape function grid cell; for each central point of the discrete grid cell inside the shape function grid cell, record its corresponding shape function grid cell number; according to the above rules, determine the shape function grid cell number where the central point of the I-th discrete grid cell is located And obtain the set e accordingly * ;

[0056] The local coordinates of the central point of the discrete grid cell are the coordinates in the local coordinate space Ω std = [-1, 1]. Map the effective elements in the shape function grid into standard shape elements in the local coordinate space, that is, map the quadrilateral element into a square element through the first-order shape function, and the coordinate range is [-1, 1]; inside the standard shape element, calculate the local coordinates of the central point of the I-th discrete grid cell through the first-order shape function weighted And obtain the set

[0057] Preferably, the element Gauss integration points described in step 4 Are the integration nodes of the Gauss-Legendre integration, Belong to the m * ×1 space, denoted as Where m * Is the number of element Gauss integration points, Represents the r-th element Gauss integration point in the e-th element; the weight of the element Gauss integration point Belong to the m *The ×1 space is denoted as is the weight corresponding to the Gaussian integration point of the r-th element in the e-th element; the Jacobian determinant value J of the element Gaussian integration point e belongs to m * The ×1 space is denoted as where J er is the Jacobian determinant value of the Gaussian integration point of the r-th element in the e-th element.

[0058] Preferably, the number N of element shape functions in step 4 * is related to the shape function order P, specifically as follows:

[0059] When P = 1, the element shape function N is a first-order shape function, N * = 4, which is consistent with the number of nodes of the shape function grid element, defined in the local coordinate space Ω std = [-1, 1]. The expressions of the 4 first-order shape functions are respectively:

[0060]

[0061] In the formula, (ξ, η) are the coordinates in the local coordinate space Ω std = [-1, 1], corresponding to (x, y) in the global coordinate space. N1(ξ, η), N2(ξ, η), N3(ξ, η), N4(ξ, η) are respectively the 1st, 2nd, 3rd, and 4th shape functions in the element shape function N;

[0062] When P = 2, the element shape function N adds 4 second-order shape functions on the basis of the previous order, that is, N * = 8, where the second-order shape functions are constructed by one-dimensional shape functions and Legendre polynomials. The expressions of the 4 second-order shape functions are respectively:

[0063]

[0064] In the formula, L0(ξ), L2(ξ) are the values of the 0th-order and 2nd-order Legendre polynomials at ξ respectively; L0(-ξ), L2(-ξ) are the values of the 0th-order and 2nd-order Legendre polynomials at -ξ respectively; L0(η), L2(η) are the values of the 0th-order and 2nd-order Legendre polynomials at η respectively; L0(-η), L2(-η) are the values of the 0th-order and 2nd-order Legendre polynomials at -η respectively; N5(ξ, η), N6(ξ, η), N7(ξ, η), N8(ξ, η) are respectively the 5th, 6th, 7th, and 8th shape functions in the element shape function N; the Legendre polynomials are realized by calling the legendreP function in the MATLAB environment;

[0065] When P = 3, the element shape function N increases by 4 third-order shape functions on the basis of the previous order, that is, N * = 12, where the expressions of the 4 third-order shape functions are respectively:

[0066]

[0067] In the formula, L1(ξ) and L3(ξ) are the values of the first-order and third-order Legendre polynomials at ξ respectively; L1(-ξ) and L3(-ξ) are the values of the first-order and third-order Legendre polynomials at -ξ respectively; L1(η) and L3(η) are the values of the first-order and third-order Legendre polynomials at η respectively; L1(-η) and L3(-η) are the values of the first-order and third-order Legendre polynomials at -η respectively; N9(ξ, η), N 10 (ξ, η), N 11 (ξ, η), N 12 (ξ, η) are the 9th, 10th, 11th, and 12th shape functions in the element shape function N respectively;

[0068] When P = 4, the element shape function N increases by 5 fourth-order shape functions on the basis of the previous order, that is, N * = 17, where the expressions of the 5 fourth-order shape functions are respectively:

[0069]

[0070] In the formula, L2(ξ) and L4(ξ) are the values of the second-order and fourth-order Legendre polynomials at ξ respectively, L2(-ξ) and L4(-ξ) are the values of the second-order and fourth-order Legendre polynomials at -ξ respectively, L2(η) and L4(η) are the values of the second-order and fourth-order Legendre polynomials at η respectively, L2(-η) and L4(-η) are the values of the second-order and fourth-order Legendre polynomials at -η respectively; N 13 (ξ, η), N 14 (ξ, η), N 15 (ξ, η), N 16 (ξ, η), N 17 (ξ, η) are the 13th, 14th, 15th, 16th, and 17th shape functions in the element shape function N respectively.

[0071] Compared with the prior art, the beneficial effects of the present invention include:

[0072] 1. Enhance the accuracy and applicability of the random field discretization of slopes: By combining with the slope numerical model, the present invention can flexibly adapt to different random field characteristics and complex boundaries, and is applicable to complex geological and structural random models. Compared with the traditional KL expansion method, the present invention effectively reduces the number of expansion terms, can reduce the number of subsequent random samplings, significantly improves the efficiency of subsequent random analysis, and is easy to be applied to fields such as numerical simulation and structural reliability analysis.

[0073] 2. Improve the efficiency and accuracy of random field discretization: The present invention innovatively combines the shape function method, the Gaussian quadrature method and the interpolation method, and proposes a new method to transform the complex multi-dimensional random field discretization and solution problem into a finite simple matrix calculation and assembly problem, significantly reducing the computational complexity and improving the computational efficiency. In addition, through the matrix processing ability of MATLAB, the present invention can greatly accelerate the calculation process while maintaining high accuracy, meeting the dual requirements of efficiency and accuracy in engineering practice.

[0074] 3. Improve the flexibility and versatility of the discretization process: The calculation accuracy can be improved by increasing the number of shape function grids or the order of shape function, and the optimal solution can be freely selected according to specific requirements. In addition, through algorithm optimization, the large matrix calculation is transformed into a block matrix calculation, effectively reducing the consumption of computing resources, so as to better adapt to the needs of large-scale data processing and complex model analysis common in engineering. Based on this, the present invention has high feasibility and computational efficiency in practical engineering applications, and is applicable to geotechnical engineering, structural engineering and other fields that require accurate simulation and analysis of random field characteristics. Description of the Drawings

[0075] Figure 1 It is a flow chart of a KL expansion random field discretization method based on the shape function method of the present invention;

[0076] Figure 2 It is a schematic diagram of the discretization grid in the embodiment of the present invention;

[0077] Figure 3 It is a schematic diagram of the shape function grid in the embodiment of the present invention;

[0078] Figure 4 It is a schematic diagram of the position relationship between the center point of the discretization grid and the shape function grid in the embodiment of the present invention;

[0079] Figure 5 It is a schematic diagram of the global coordinate transformation to the local coordinate in the embodiment of the present invention;

[0080] Figure 6 It is a schematic diagram of the element shape function in the embodiment of the present invention;

[0081] Figure 7 It is the calculation result of the eigenvalues of four shape function orders under the shape function grid h-1;

[0082] Figure 8 For the calculation results of the eigenvalues of four shape function meshes under the shape function order p - 1;

[0083] Figure 9 For the schematic diagram of the eigenfunction results in the embodiment of the present invention;

[0084] Figure 10 For the schematic diagram of the comparison results of the discrete errors in the embodiment of the present invention;

[0085] Figure 11 For the schematic diagram of the discrete results of the cohesion c in the embodiment of the present invention;

[0086] Figure 12 For the internal friction angle in the embodiment of the present invention;

[0087] Figure 13 For the statistical distribution characteristic curve of the cohesion c obtained by performing 1000 random samplings;

[0088] Figure 14 For the internal friction angle obtained by performing 1000 random samplings; Specific implementation manner

[0089] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings.

[0090] Figure 1 For the flowchart of the shape function KL expansion random field discretization method considering the irregular domain of slope soil mass in the present invention, as Figure 1 can be seen, it includes the following steps:

[0091] Step 1, establish a slope numerical model and divide the discrete grid and the shape function grid

[0092] Step 1.1, establish a slope numerical model in the finite element software ABAQUS

[0093] The slope is composed of the upper slope body and the lower soil mass. The cross-section of the slope body part is trapezoidal, and the cross-section of the soil mass part is rectangular; establish a slope model in the finite element software ABAQUS. The specific process is as follows: set the slope model in a section perpendicular to the earth's surface. The slope model includes the upper slope body part and the lower soil mass part, and set the bottom side length of the soil mass to be horizontal and the left side length to be vertical. Specifically, denote the lower left corner endpoint of the soil mass as endpoint A S , the width of the soil mass in the horizontal direction is K S , the height in the vertical direction is H S , the upper width of the slope body part is KF , the lower width is K M , and the height in the vertical direction is H F ; Taking the endpoint A S as the origin, with the horizontal direction as the x-axis and the vertical direction as the y-axis, an xy-plane coordinate system is established.

[0094] Step 1.2, dividing the discrete grid and the shape function grid

[0095] The slope numerical model is discretized to generate a discrete grid and a shape function grid respectively. Both grids use quadrilateral elements; the inp files corresponding to the two grids are exported in ABAQUS, and the coordinate information of the grid nodes and the element node information are extracted from them; the grid node coordinate information records the geometric positions of all nodes in the space coordinate system, and the element node information records the node numbers associated with each grid element and their connection relationships.

[0096] The slope model established in this embodiment is shown in Figure 2 . In the embodiment of the present invention, K S = 10m, H s = 2m, K F = 3m, K M = 8m, H F = 4m.

[0097] In this embodiment, 181 grid elements and 212 grid nodes are set for the discrete grid in ABAQUS. To illustrate the influence of the shape function grid on the discretization accuracy, the slope numerical model is divided into four shape function grids with different element sizes, as shown in detail in Figure 3 . Among them, Figure 3 a. The approximate size of the grid elements is 1, with a total of 34 grid elements and 48 grid nodes, denoted as h-1; Figure 3 b. The approximate size of the grid elements is 0.75, with a total of 89 grid elements and 111 grid nodes, denoted as h-0.75; Figure 3 c. The approximate size of the grid elements is 0.5, with a total of 181 grid elements and 212 grid nodes, denoted as h-0.5; Figure 3 d. The approximate size of the grid elements is 0.25, with a total of 754 grid elements and 816 grid nodes, denoted as h-0.25. The inp files of the grids in ABAQUS are exported, and the corresponding node coordinate information and the element node composition information are extracted.

[0098] Step 2, calculating the center points of the discrete grid elements and the effective element indices of the shape function grid

[0099] Denote the set of the center point coordinates of the discrete grid elements as set u * , Set u *Belonging to E * ×1 space, denoted as where E * is the number of discrete grid cells, I is the discrete grid cell number, represents the coordinates of the center point of the I-th discrete grid cell, and the superscript "T" represents the transpose of the matrix; the set of coordinates of the center points of the discrete grid cells in the local coordinate space is denoted as the set The set belongs to E * ×1 space, denoted as where represents the local coordinates of the center point of the I-th discrete grid cell; the set of areas of the discrete grid cells is denoted as the set S * , The set S * belongs to E * ×1 space, denoted as where represents the area of the I-th discrete grid cell.

[0100] The set of the effective element indices of the shape function grid is denoted as the set e * , The set e * belongs to E * ×1 space, denoted as where represents the number of the shape function grid cell where the center point of the I-th discrete grid cell is located.

[0101] In this embodiment, the implementation process of step 2 is as follows:

[0102] Calculate the coordinates of the center point of the I-th discrete grid cell according to the node coordinate information of the discrete grid described in step 1 and the area of the discrete grid cell and obtain the set u * and the set The coordinates of the center point of the discrete grid cell are the average of the node coordinates of each discrete grid cell; the area of the discrete grid cell is the area of each quadrilateral cell;

[0103] The effective element index refers to the number of the shape function grid cell that contains at least one center point of the discrete grid cell inside. The specific determination method is as follows:

[0104] Using a geometric determination method, evaluate the spatial relationship between the center points of each discrete grid cell and each shape function grid cell; check the center point of each discrete grid cell to determine whether it is located inside the shape function grid cell; when the center point of the discrete grid cell is on the left boundary or the lower boundary of the shape function grid cell, determine that it is located inside the shape function grid cell; for each center point of the discrete grid cell located inside the shape function grid cell, record the corresponding shape function grid cell number; according to the above rules, determine the shape function grid cell number where the center point of the I-th discrete grid cell is located And obtain set e accordingly * ;

[0105] The local coordinates of the center point of the discrete grid cell are the coordinates within the local coordinate space Ω std = [-1, 1]. Map the effective cells in the shape function grid into standard shape cells in the local coordinate space, that is, map the quadrilateral cells into square cells through the first-order shape function, and the coordinate range is [-1, 1]; within the standard shape cell, calculate the local coordinates of the center point of the I-th discrete grid cell through the weighted first-order shape function And obtain the set

[0106] Figure 4 Shows the spatial relationship between the center points of discrete grid cells and four shape function grid cells. The black dots in the figure represent the centers of the discrete grids. Among them, Figure 4 a shows the spatial relationship between the center point of the discrete grid and the shape function grid h-1; Figure 4 b shows the spatial relationship between the center point of the discrete grid and the shape function grid h-0.75; Figure 4 c shows the spatial relationship between the center point of the discrete grid and the shape function grid h-0.5; Figure 4 d shows the spatial relationship between the center point of the discrete grid and the shape function grid h-0.25. Figure 5 Shows the coordinate transformation process of converting the center point of the discrete grid cell in the shape function grid cell to the local coordinates.

[0107] Step 3, establish the probability distribution model of the soil parameter random field

[0108] Based on the slope numerical model established in Step 1, further construct the probability distribution model of the soil parameter random field. Among them, the soil parameters involved include but are not limited to cohesion, internal friction angle, elastic modulus, Poisson's ratio, and density; the components of this probability distribution model of the random field cover the type, mean, standard deviation, and autocorrelation function of the probability distribution of the random field.

[0109] Let (x, y) and (x′, y′) be the coordinates of the centers of any two discrete grid cells in the slope numerical model, and the soil parameter random field be the random field H(x, y), whose probability distribution type is lognormal distribution.

[0110] Let the mean of the random field H(x, y) be μ(x, y), the standard deviation be σ(x, y), and the logarithmic mean be μ LN (x, y), μ LN (x, y) = lnμ(x, y) - [ln(1 + δ 2 (x, y))] / 2, where δ(x, y) is the coefficient of variation; let the logarithmic standard deviation of the random field H(x, y) be σ LN (x, y), σ LN (x, y) = [ln(1 + δ 2 (x,y))] 1 / 2 .

[0111] Let the autocorrelation function of the random field H(x, y) be the autocorrelation function ρ((x, y), (x′, y′)), and this autocorrelation function ρ((x, y), (x′, y′)) is a single-exponential function, and its expression is:

[0112]

[0113] In the formula, L h is the horizontal autocorrelation distance, L v is the vertical autocorrelation distance;

[0114] In this embodiment, the random field parameters of the slope soil are the cohesion c and the internal friction angle The mean of the cohesion c is μ(x, y) = 34 kPa, and the coefficient of variation δ(x, y) = 0.37; the internal friction angle has a mean of μ(x, y) = 23°, and the coefficient of variation δ(x, y) = 0.19. L h = 20 m, L v = 2 m.

[0115] Step 4, calculate the element stiffness submatrix k e and the element mass matrix E e

[0116] The element stiffness submatrix k e and the element mass matrix E e correspond to each shape function grid cell and are both calculated in the local coordinate space. Select consistent element Gaussian integration points and shape function order P for each shape function grid cell, and their expressions are respectively:

[0117]

[0118] In the formula, is the weight of the unit Gauss integration point, and J e is the value of the Jacobi determinant at the unit Gauss integration point, and diag(·) is the symbol for converting a vector into a diagonal matrix; N e is the unit shape function matrix, and is the matrix composed of the values of the unit shape function N at the unit Gauss integration point ; the unit shape function N belongs to the N * ×1 space, denoted as N * is the number of unit shape functions;

[0119] The unit stiffness sub-matrix k e belongs to the m * ×N * space, denoted as where m * is the number of unit Gauss integration points; the unit mass matrix E e belongs to the N * ×N * space, denoted as

[0120] In this embodiment, the unit Gauss integration point described in step 4 is the integration node of the Gauss-Legendre integration, belongs to the m * ×1 space, denoted as where m * is the number of unit Gauss integration points, represents the r-th unit Gauss integration point in the e-th unit; the weight of the unit Gauss integration point belongs to the m * ×1 space, denoted as is the weight corresponding to the r-th unit Gauss integration point in the e-th unit; the value of the Jacobi determinant J of the unit Gauss integration point e belongs to the m * ×1 space, denoted as where J er is the value of the Jacobi determinant of the r-th unit Gauss integration point in the e-th unit.

[0121] In this embodiment, the number of unit shape functions N * is related to the shape function order P, specifically as follows:

[0122] When P = 1, the unit shape function N is a first-order shape function, N * = 4, which is consistent with the number of nodes of the shape function grid element, and is defined in the local coordinate space Ω Std = [-1, 1], and the expressions of the 4 first-order shape functions are respectively:

[0123]

[0124] wherein, (ξ, η) are the coordinates in the local coordinate space Ω std = [-1, 1], corresponding to (x, y) in the global coordinate space. N1(ξ, η), N2(ξ, η), N3(ξ, η), N4(ξ, η) are the 1st, 2nd, 3rd, and 4th shape functions in the element shape function N respectively;

[0125] When P = 2, the element shape function N adds 4 second-order shape functions on the basis of the previous order, i.e., N * = 8, where the second-order shape functions are constructed by one-dimensional shape functions and Legendre polynomials. The expressions of the 4 second-order shape functions are respectively:

[0126]

[0127] wherein, L0(ξ), L2(ξ) are the values of the 0th-order and 2nd-order Legendre polynomials at ξ respectively; L0(-ξ), L2(-ξ) are the values of the 0th-order and 2nd-order Legendre polynomials at -ξ respectively; L0(η), L2(η) are the values of the 0th-order and 2nd-order Legendre polynomials at η respectively; L0(-η), L2(-η) are the values of the 0th-order and 2nd-order Legendre polynomials at -η respectively; N5(ξ, η), N6(ξ, η), N7(ξ, η), N8(ξ, η) are the 5th, 6th, 7th, and 8th shape functions in the element shape function N respectively; The Legendre polynomials are implemented by calling the legendreP function in the MATLAB environment;

[0128] When P = 3, the element shape function N adds 4 third-order shape functions on the basis of the previous order, i.e., N * = 12, wherein, the expressions of the 4 third-order shape functions are respectively:

[0129]

[0130] wherein, L1(ξ), L3(ξ) are the values of the 1st-order and 3rd-order Legendre polynomials at ξ respectively; L1(-ξ), L3(-ξ) are the values of the 1st-order and 3rd-order Legendre polynomials at -ξ respectively; L1(η), L3(η) are the values of the 1st-order and 3rd-order Legendre polynomials at η respectively; L1(-η), L3(-η) are the values of the 1st-order and 3rd-order Legendre polynomials at -η respectively; N9(ξ, η), N 10 (ξ, η), N 11 (ξ, η), N 12($\xi$, $\eta$) are the 9th, 10th, 11th, and 12th shape functions within the element shape function $N$ respectively;

[0131] When $P = 4$, the element shape function $N$ adds 5 fourth-order shape functions on the basis of the previous order, i.e., $N$ * $= 17$, where the expressions of the 5 fourth-order shape functions are respectively:

[0132]

[0133] In the formula, $L_2(\xi)$ and $L_4(\xi)$ are the values of the 2nd-order and 4th-order Legendre polynomials at $\xi$ respectively, $L_2(-\xi)$ and $L_4(-\xi)$ are the values of the 2nd-order and 4th-order Legendre polynomials at $-\xi$ respectively, $L_2(\eta)$ and $L_4(\eta)$ are the values of the 2nd-order and 4th-order Legendre polynomials at $\eta$ respectively, and $L_2(-\eta)$ and $L_4(-\eta)$ are the values of the 2nd-order and 4th-order Legendre polynomials at $-\eta$ respectively; $N$ 13 ($\xi$, $\eta$), $N$ 14 ($\xi$, $\eta$), $N$ 15 ($\xi$, $\eta$), $N$ 16 ($\xi$, $\eta$), $N$ 17 ($\xi$, $\eta$) are the 13th, 14th, 15th, 16th, and 17th shape functions within the element shape function $N$ respectively.

[0134] In this embodiment, Figure 6 shows the schematic diagrams of shape functions under different orders. The first row corresponds to $P = 1$, and from left to right are $N_1(\xi, \eta)$, $N_2(\xi, \eta)$, $N_3(\xi, \eta)$, and $N_4(\xi, \eta)$. The second row corresponds to $P = 2$, and from left to right are the schematic diagrams of $N_5(\xi, \eta)$, $N_6(\xi, \eta)$, $N_7(\xi, \eta)$, and $N_8(\xi, \eta)$. The third row corresponds to $P = 3$, and from left to right are $N_9(\xi, \eta)$, $N$ 10 ($\xi$, $\eta$), $N$ 11 ($\xi$, $\eta$), and $N$ 12 ($\xi$, $\eta$). The fourth row corresponds to $P = 4$, and from left to right are $N$ 13 ($\xi$, $\eta$), $N$ 14 ($\xi$, $\eta$), $N$ 15 ($\xi$, $\eta$), $N$ 16 ($\xi$, $\eta$), and $N$ 17 ($\xi$, $\eta$) schematic diagrams. Among them, the first row are all nodal shape functions, related to the number of nodes of the shape function grid unit, $N$ 17 ($\xi$, $\eta$) is a surface function, related to the number of surfaces of the shape function grid unit, and the others are all edge functions, related to the number of edges of the shape function grid unit.

[0135] To illustrate the influence of high-order shape functions on the discrete accuracy, the shape functions of the lower elements of orders 1, 2, 3, and 4 are selected for the shape function mesh h-1 respectively, and the shape functions of the lower elements of only order 1 are selected for other shape function meshes.

[0136] In this embodiment, the number of Gaussian integration points in each dimension is set to be 1 more than the order of the shape function, and the numbers of integration points are 2, 3, 4, and 5 in sequence; correspondingly, the sizes of the two-dimensional element shape function matrices under the four orders are 4×4, 9×8, 16×12, and 25×17 in sequence; when the order of the shape function p = 1, N e The specific values are as follows:

[0137]

[0138] In this embodiment, corresponding to the numbers of Gaussian integration points under the four orders, the unit Jacobi determinant matrix J e has sizes of 4×1, 9×1, 16×1, and 25×1 in sequence. When the order of the shape function p = 1, the Jacobi determinant matrix J e = [0.0652 0.06260.0652 0.0626] T ; corresponding to the unit shape function matrices N e under the four orders, the unit stiffness sub-matrix k e has sizes of 4×4, 9×8, 16×12, and 25×17 in sequence; when the order of the shape function p = 1, k e The specific values are as follows:

[0139]

[0140] In this embodiment, corresponding to the unit shape function matrices N e under the four orders, the unit mass matrix E e has sizes of 4×4, 8×8, 12×12, and 17×17 in sequence; when the order of the shape function p = 1, E e The specific values are as follows:

[0141]

[0142] Step 5, calculate the autocorrelation matrix ρ and the global stiffness matrix B

[0143] First, generate a shape function index, which is the number of all shape functions of the shape function mesh, denoted as IEN, belonging to e s ×N * space, where e sLet \(n\) be the number of shape function grid cells; the IEN generation rule depends on the order of the shape function; when the order of the shape function \(P = 1\), the value of IEN corresponds to the node number of the shape function grid cell; when \(P = 2, 3\), on the basis of the previous order number, IEN increases by the number of each side of the shape function grid cell; when \(P = 4\), on the basis of the previous order number, IEN increases by the side number and face number of the shape function grid cell.

[0144] The autocorrelation matrix \(\rho\) is a matrix calculated by the autocorrelation function for all element Gaussian integration points. Specifically, the coordinates of all element Gaussian integration points are mapped from local coordinates to the global coordinate system, and the set of coordinates of all element Gaussian integration points is denoted as \(z\). * , where is the coordinate of the \(t\)-th all-element Gaussian integration point, and \(M\) * is the number of all-element Gaussian integration points; finally, the autocorrelation matrix \(\rho\) is calculated by the autocorrelation function for all element Gaussian integration points.

[0145] The global stiffness matrix \(B\) belongs to the \(N\) s × \(N\) s space, denoted as is defined by the global stiffness submatrix \(K\) and the autocorrelation matrix \(\rho\), \(B = K\) T \(\rho K\), where \(N\) s is the number of all shape functions of the shape function grid; the global stiffness submatrix \(K\) is assembled from the element stiffness submatrix \(k\) e according to the shape function index IEN.

[0146] In this embodiment, the global stiffness submatrix \(K\) is assembled from the elements in the element stiffness submatrix \(k\) e in the order of the shape function grid cells according to the shape function index IEN, and belongs to the \(M\) * × \(N\) s space, denoted as The correspondence rule between the elements in the global stiffness submatrix \(K\) and those in \(k\) e is as follows:

[0147]

[0148] \(i=(e - 1)m\) * +\(i\) *

[0149]

[0150] In the formula, are the elements in the \(i\)-th row and \(j\)-th column of the global stiffness submatrix \(K\) in the \(e\)-th and \((e - 1)\)-th iterations respectively, where \(i\) is determined by the \(i\)-th * element Gaussian integration point in the \(e\)-th shape function grid cell; \(j\) is determined by the \(j\)-th* The shape function index IEN corresponding to each shape function is determined; this process is repeated until all elements in the element stiffness matrix k e are assembled into the global stiffness submatrix K;

[0151] The expression of the autocorrelation matrix ρ is:

[0152]

[0153] In the formula, sym. represents a symmetric matrix; represents the t-th global element Gauss integration point and the t'-th global element Gauss integration point the autocorrelation function value obtained through the autocorrelation function.

[0154] In this embodiment, when the shape function mesh h - 1 and P are 1, 2, 3, 4 in sequence, the sizes of ρ are 136×136, 306×306, 544×544, 850×850 respectively; when the shape function mesh h - 0.75 and P = 1, the size of ρ is 356×356; when the shape function mesh h - 0.5 and p = 1, the size of ρ is 724×724; when the shape function mesh h - 0.25 and P = 1, the size of ρ is 3016×3016.

[0155] In this embodiment, when the shape function mesh h - 1 and P are 1, 2, 3, 4 in sequence, the sizes of K are 136×48, 306×129, 544×210, 850×325 respectively; when the shape function mesh h - 0.75 and P = 1, the size of K is 356×111; when the shape function mesh h - 0.5 and P = 1, the size of K is 724×212; when the shape function mesh h - 0.25 and P = 1, the size of K is 3016×816.

[0156] In this embodiment, when the shape function mesh h - 1 and P are 1, 2, 3, 4 in sequence, the sizes of B are 48×48, 129×129, 210×210, 325×325 respectively; when the shape function mesh h - 0.75 and P = 1, the size of B is 111×111; when the shape function mesh h - 0.5 and P = 1, the size of B is 212×212; when the shape function mesh h - 0.25 and P = 1, the size of B is 816×816.

[0157] Step 6, calculate the eigenvalue λ k eigenfunction φ k and the number of expansion terms M

[0158] Step 6.1, introduce the eigenvalue matrix A, where E g is the global mass matrix, composed of the element mass matrix E eAssembled according to the shape function index IEN, where the superscript "-1" represents the inverse operation of the matrix; the eigenvalue matrix A belongs to N s ×N s space, denoted as

[0159] The eigenvalue λ k is the k-th eigenvalue in the eigenvalue vector λ obtained by sorting the eigenvalues of the eigenvalue matrix A in descending order, λ belongs to N s ×1 space, denoted as

[0160] The eigenfunction corresponding to the eigenvalue λ k is denoted as φ k , belonging to E * ×1 space, denoted as where D u is the eigenfunction matrix, and D uk is the column vector of D u .

[0161] In this embodiment, the global mass matrix E g is assembled from the elements in the element mass matrix E e , belonging to N s ×N s space, denoted as The corresponding assembly rule can be written as the following formula:

[0162]

[0163] In the formula, are the elements in the ii-th row and jj-th column of the global mass matrix E g in the e-th and (e - 1)-th iterations respectively, where ii and jj are determined by the shape function index IEN corresponding to the ii * -th, and the jj * -th shape functions in the e-th shape function element; repeat this process until all elements in the element mass matrix E e are assembled into the global mass matrix E g .

[0164] The calculation process of the eigenfunction matrix D u is as follows:

[0165] Introduce the discrete point element eigenfunction matrix D ue , belonging to E * ×N * space, is the matrix composed of the values of the element shape function N at the local coordinates of the center point of the discrete grid element; through the shape function index IEN, Due Assembled into the characteristic function matrix D u , D u The column elements in it are the eigenvalues λ k The corresponding characteristic functions, denoted as D uk .

[0166] Step 6.2, calculate the number of expansion terms M

[0167] Step 6.2.1, initialize the number of expansion terms M of the KL method to 1;

[0168] Step 6.2.2, denote the discretization error of the number of expansion terms M as the discretization error ε M , and its expression is:

[0169]

[0170] Where, is the autocorrelation function value corresponding to the center points of the I-th and I'-th discrete grid cells; are the characteristic function values corresponding to the center points of the I-th and I'-th discrete grid cells in the k-th characteristic function φ k respectively; S is the total area of the discrete grid region; are the areas of the I-th and I'-th discrete grid cells respectively;

[0171] Step 6.2.3, compare the discretization error ε M with the given allowable discretization error ε a : If ε M > ε a , then increment the number of expansion terms M by 1 and return to Step 6.2.2; if ε M ≤ ε a , then the iteration ends and the number of expansion terms M is output.

[0172] In this embodiment, the calculation results of the first 20 eigenvalues λ are shown in Figure 7 and Figure 8 . Where, Figure 7 are the calculation results of the eigenvalues of the four shape function orders under the shape function grid h-1, Figure 8The calculation results of the eigenvalues of four shape function meshes under the shape function order p - 1 are shown. It can be seen from the figure that when h - 1, p - 3 and h - 0.5, p - 1, the calculation results of the eigenvalues have highly converged. The sizes of the eigenvalue matrices to be solved are 325×325 and 212×212 respectively. Through simple matrix assembly operations, the random field can be discretized, significantly improving the efficiency and accuracy of the random field discretization. The final convergence results of the eigenvalues are λ1 = 17.7754, λ2 = 6.9628, λ3 = 3.2871, λ4 = 1.8570, λ5 = 1.4211, λ6 = 1.1160, λ7 = 0.7626, λ8 = 0.5702, λ9 = 0.5388, λ 10 = 0.4076, λ 11 = 0.3237, λ 12 = 0.3033, λ 13 = 0.2667, λ 14 = 0.2533, λ 15 = 0.2091, λ 16 = 0.1927, λ 17 = 0.1735, λ 18 = 0.1485, λ 19 = 0.1451 and λ 20 = 0.1451.

[0173] In this embodiment, the calculation results of the normalized eigenfunctions φ k corresponding to the first 20 eigenvalues λ k are shown in Figure 9 . Figure 9 Among them, 1 - 20 are the calculation results of φ1 - φ 20 respectively.

[0174] In this embodiment, Figure 10 shows the comparison between the KL method and the method of the present invention, where KL - 1 is the KL method and KL - 2 is the discrete method proposed in this paper. It can be seen from Figure 10 that the discrete error ε M obtained by the method of the present invention is always smaller than that of the KL method, and can effectively reduce the number of expansion terms, thereby improving the accuracy and applicability of the random field discretization of the slope. Let the allowable discrete error ε a = 1%, and the calculation results of the discrete error ε M in the iterative process of selecting the number of expansion terms M are: ε1 = 17.01%, ε2 = 8.86%, ε3 = 5.86%, ε4 = 4.38%, ε5 = 3.63%, ε6 = 2.90%, ε7 = 2.40%, ε8 = 2.13%, ε9 = 1.84%, ε 10 = 1.60%, ε 11 = 1.34%, ε 12= 1.22%, ε 13 = 1.12%, ε 14 = 0.94%, ε 15 = 0.93%, ε 16 = 0.86%, ε 17 = 0.83%, ε 18 = 0.81%, ε 19 = 0.78% and ε 20 = 0.75%, where ε 14 ≤ 1%, and finally obtained M = 14.

[0175] Step 7, realizing the discretization of the random field of the slope numerical model

[0176] The random field is discretized into an independent standard normal random variable ζ and an eigenvalue λ k , and a characteristic function φ k multiplied to obtain the estimated value of the random field at each point; the independent standard normal random variable ζ = [ζ1, ζ2,..., ζ k ,..., ζ M T , representing the randomness of the geotechnical parameters, belonging to the M×1 space, denoted as ζ ∈ R M×1 , where ζ k is the random variable corresponding to the eigenvalue λ k , with a mean of 0 and a standard deviation of 1;

[0177] For the geotechnical parameters that conform to the normal random distribution, their mean is denoted as the normal mean μ(u * ), the variance is denoted as the normal variance σ(u * ), and the discretized value of the random field is denoted as the normal discretized value H(u * ). The discretization expression of the random field of this geotechnical parameter is:

[0178]

[0179] In the formula, ζ kn is the nth sampling result of the kth random variable ζ k ;

[0180] For the geotechnical parameters that conform to the lognormal distribution, their mean is denoted as the lognormal mean μ L (u * ), the variance is denoted as the lognormal variance σ L (u * ), the discretized value of the random field is denoted as the lognormal discretized value H L (u * ). The discretization expression of the random field of this geotechnical parameter is: ​

[0181]

[0182] In this embodiment, for the random field parameters of the slope numerical model, the cohesion c and the internal friction angle The discrete results of random sampling are shown in Figure 11 and Figure 12 , where Figure 11 are the nine random results of the cohesion c, Figure 12 are the nine random results of the internal friction angle .

[0183] Figure 13 is the statistical distribution characteristic curve of the cohesion c statistically obtained from 1000 random samplings. The statistical mean μ(x, y) is 33.81, and the coefficient of variation δ(x, y) is 0.3714, which is highly consistent with the preset mean μ(x, y) = 34 and the preset coefficient of variation δ(x, y) = 0.37. Figure 14 is the statistical distribution characteristic curve of the internal friction angle statistically obtained from 1000 random samplings. The statistical mean μ(x, y) is 23.0216, and the coefficient of variation δ(x, y) is 0.1915, which is highly consistent with the preset mean μ(x, y) = 23 and the preset coefficient of variation δ(x, y) = 0.19, indicating the correctness of this method. By adopting the method of the present invention in combination with the numerical model of the irregular random field, it can effectively adapt to different random field characteristics and boundary conditions, and is widely applicable to complex geological conditions and structural random models, providing an effective solution for the discretization of the irregular slope random field.

Claims

1. A discrete method for KL expansion random field of shape functions considering irregular domains of slope soil masses, characterized in that, It includes the following steps: Step 1: Establish a slope numerical model and divide discrete grids and shape function grids Step 1.1: Establish a slope numerical model in the finite element software ABAQUS The slope is composed of the upper slope body and the lower soil body. The cross-section of the slope body part is trapezoidal, and the cross-section of the soil body part is rectangular. The slope model is established in the finite element software ABAQUS. The specific process is as follows: The slope model is set in a section perpendicular to the earth's surface. The slope model includes the upper slope body part and the lower soil body part. The bottom side length of the soil body is set in the horizontal direction, and the left side length is set in the vertical direction. Specifically, the lower left corner endpoint of the soil body is denoted as endpoint A s , the width of the soil body in the horizontal direction is K s , and the height in the vertical direction is H s , the upper width of the slope body part is K F , the lower width is K M , and the height in the vertical direction is H F ; Taking endpoint A s as the origin, the horizontal direction as the x direction, and the vertical direction as the y direction, an xy plane coordinate system is established; Step 1.2: Divide discrete grids and shape function grids Perform discretization processing on the slope numerical model to generate discrete grids and shape function grids respectively. Both types of grids use quadrilateral elements; Export the inp files corresponding to the two types of grids in ABAQUS, and extract the grid node coordinate information and element node information from them; the grid node coordinate information records the geometric positions of all nodes in the spatial coordinate system, and the element node information records the node numbers associated with each grid element and their connection relationships; Step 2: Calculate the center points of discrete grid elements and the effective element indices of shape function grids Denote the set of the central point coordinates of discrete grid cells as set \(u\). * , Set \(u\) * belongs to the \(E\times1\) space, denoted as * where \(E\) is the number of discrete grid cells, \(I\) is the discrete grid cell number, * , and \((x_{I})\) represents the central point coordinates of the \(I\)-th discrete grid cell, and the superscript "\(T\)" represents the transpose of the matrix; Denote the set of the coordinates of the central points of discrete grid cells in the local coordinate space as set where \((\xi_{I})\) represents the local coordinates of the central point of the \(I\)-th discrete grid cell; Denote the set of the areas of discrete grid cells as set \(S\). Set belongs to the \(E\times1\) space, denoted as * where and \((\xi_{I})\) represents the local coordinates of the central point of the \(I\)-th discrete grid cell; Denote the set of the areas of discrete grid cells as set \(S\). where \((S_{I})\) represents the area of the \(I\)-th discrete grid cell; * , Set \(S\) * belongs to the \(E\times1\) space, denoted as * where and \((S_{I})\) represents the area of the \(I\)-th discrete grid cell; ​ Denote the set of valid element indices of the shape function grid as set e * , Set e * belongs to E * ×1 space, denoted as where represents the number of the shape function grid cell where the center point of the I-th discrete grid cell is located; Step 3: Establish a probability distribution model of the random field of soil parameters Based on the slope numerical model established in Step 1, further construct a probability distribution model of the random field of soil parameters. Among them, the soil parameters involved include but are not limited to cohesion, internal friction angle, elastic modulus, Poisson's ratio, and density; the components of this probability distribution model of the random field cover the type, mean, standard deviation, and autocorrelation function of the probability distribution of the random field; Denote (x, y) and (x′, y′) as the coordinates of the center points of any two discrete grid elements in the slope numerical model. The random field of soil parameters is the random field H(x, y), and its probability distribution type is lognormal distribution; Let the mean of the random field H(x, y) be μ(x, y), the standard deviation be σ(x, y), and the logarithmic mean be μ LN (x, y), μ LN (x, y) = lnμ(x, y) - [ln(1 + δ 2 (x, y))] / 2, where δ(x, y) is the coefficient of variation; let the logarithmic standard deviation of the random field H(x, y) be σ LN (x, y), σ LN (x, y) = [ln(1 + δ 2 (x, y))] 1 / 2 ; Denote the autocorrelation function of the random field H(x, y) as the autocorrelation function ρ((x, y), (x′, y′)). This autocorrelation function ρ((x, y), (x′, y′)) is a single-exponential function, and its expression is: where L h is the horizontal autocorrelation distance, and L v is the vertical autocorrelation distance; Step 4, calculate the element stiffness sub-matrix k e and the element mass matrix E e The unit stiffness sub-matrix k e and the unit mass matrix E e correspond to each shape function mesh element and are both calculated in the local coordinate space. Consistent unit Gauss integration points are selected for each shape function mesh element and the shape function order P, and their expressions are respectively: In the formula, is the weight of the unit Gauss integration point, and J e is the value of the Jacobi determinant at the unit Gauss integration point, and diag(·) is the symbol for converting a vector into a diagonal matrix; N e is the unit shape function matrix, which is the matrix composed of the values of the unit shape function N at the unit Gauss integration point . The unit shape function N belongs to the N * ×1 space, denoted as N * is the number of unit shape functions; The element stiffness sub-matrix k e belongs to m * ×N * space, denoted as where m * is the number of element Gauss integration points; the element mass matrix E e belongs to N * ×N * space, denoted as Step 5: Calculate the autocorrelation matrix ρ and the global stiffness matrix B First, generate a shape function index, which is the numbering of all shape functions in the shape function grid, denoted as IEN, belonging to the e s ×N * space, where e s is the number of shape function grid cells; the generation rule of IEN depends on the shape function order; and when the shape function order P = 1, the value of IEN corresponds to the node number of the shape function grid cell; when P = 2, 3, IEN is incremented by the number of each side of the shape function grid cell based on the previous order numbering; when P = 4, IEN is incremented by the side number and face number of the shape function grid cell based on the previous order numbering. The autocorrelation matrix ρ is a matrix obtained by calculating the autocorrelation function for all element Gaussian integration points. Specifically, the coordinates of all element Gaussian integration points are mapped from local coordinates to the global coordinate system, and the set of coordinates of all element Gaussian integration points is denoted as z*, where is the coordinate of the t-th all-element Gaussian integration point, M* is the number of all-element Gaussian integration points; finally, the autocorrelation matrix ρ is calculated by the autocorrelation function for all element Gaussian integration points. The global stiffness matrix B belongs to N s ×N s space, denoted as defined by the global stiffness submatrix K and the autocorrelation matrix ρ, B = K T ρK, where N s is the total number of shape functions of the shape function mesh; the global stiffness submatrix K is assembled from the element stiffness submatrix k e according to the shape function index IEN; Step 6, calculate the eigenvalue λ k , the eigenfunction φ k and the number of expansion terms M Step 6.1, introduce the eigenvalue matrix A, where E g is the global mass matrix, assembled from the element mass matrix E E according to the shape function index IEN, and the superscript "-1" represents the inverse operation of the matrix; the eigenvalue matrix A belongs to N s ×N s space, denoted as The eigenvalue λ k is the k-th eigenvalue in the eigenvalue vector λ obtained by sorting the eigenvalue matrix A in descending order after eigen-decomposition, λ belongs to N s ×1 space, denoted as Let the eigenvalue be λ k and denote the corresponding eigenfunction as φ k , which belongs to the E * ×1 space, denoted as where D u is the eigenfunction matrix, and D uk is the column vector of D u ; Step 6.2: Calculate the expansion term number M Step 6.2.1: Initialize the KL method expansion term number M to 1; Step 6.2.2, denote the discrete error of the number of expansion terms M as the discrete error ε M , and its expression is: Among them, is the autocorrelation function value corresponding to the center point of the I-th discrete grid cell and the center point of the I'-th discrete grid cell; are respectively the eigenfunction values corresponding to the center points of the I-th and I'-th discrete grid cells in the k-th eigenfunction φ k ; S is the total area of the discrete grid region; are respectively the areas of the I-th and I'-th discrete grid cells; Step 6.2.3, take the discrete error ε M and compare it with the given allowable discrete error ε a : If ε M > ε a , then increment the number of expansion terms M by 1 and return to Step 6.2.2; if ε M ≤ ε A , then end the iteration and output the number of expansion terms M; Step 7: Realize the discretization of the random field of the slope numerical model The random field is discretized into independent standard normal random variables ζ and eigenvalues λ k , and the characteristic function φ k are multiplied to obtain the estimated value of the random field at each point; the independent standard normal random variable ζ = [ζ1, ζ2,..., ζ k ,..., ζ M T , representing the randomness of geotechnical parameters, belonging to the M×1 space, denoted as ζ ∈ R M×1 , where ζ k is the random variable corresponding to the eigenvalue λ k , with a mean of 0 and a standard deviation of 1;​ For the geotechnical parameters that conform to the normal random distribution, their mean value is denoted as the normal mean μ (u * ), the variance is denoted as the normal variance σ (u * ), and the discrete value of the random field is denoted as the normal discrete value H (u * ). The discrete expression of the random field of this geotechnical parameter is as follows: where ζ kn is the n-th sampling result of the k-th random variable ζ k ; For the geotechnical parameters that conform to the lognormal distribution, the mean is denoted as the lognormal mean μ L (u * ), and the variance is denoted as the lognormal variance σ L (u * ). The discrete value of the random field is denoted as the lognormal discrete value H L (u * ), and the discrete expression of the random field of this geotechnical parameter is as follows:

2. A discrete method for KL expansion random field of shape function considering irregular domain of slope soil body according to claim 1, characterized in that The implementation process of Step 2 is as follows: Calculate the center point coordinates of the $I$-th discrete grid cell based on the node coordinate information of the discrete grid described in step 1 and the area of the discrete grid cell and obtain the set $u$ * and the set The center point coordinates of the discrete grid cell are the average of the node coordinates of each discrete grid cell; the area of the discrete grid cell is the area of each quadrilateral cell The effective element index refers to the shape function grid element number that contains at least one center point of a discrete grid element inside. The specific determination method is: Using a geometric determination method, evaluate the spatial relationship between the center point of each discrete grid cell and each shape function grid cell; check the center point of each discrete grid cell to determine whether it is located inside the shape function grid cell; when the center point of the discrete grid cell is on the left boundary or lower boundary of the shape function grid cell, determine that it is located inside the shape function grid cell; for each center point of the discrete grid cell located inside the shape function grid cell, record the corresponding shape function grid cell number; according to the above rules, determine the shape function grid cell number where the center point of the I-th discrete grid cell is located And thus obtain the set e*; The local coordinates of the center point of the discrete grid element are the coordinates within the local coordinate space Ω std = Coordinates within [-1, 1]. Map the effective elements in the shape function grid into standard shape elements in the local coordinate space, that is, map the quadrilateral elements into square elements through the first-order shape function, and the coordinate range is [-1, 1]; within the standard shape element, calculate the local coordinates of the center point of the I-th discrete grid element by weighting through the first-order shape function and obtain the set 3. A discrete method for KL expansion random field of shape function considering irregular domain of slope soil mass according to claim 1, characterized in that The unit Gauss integration points described in step 4 are the integration nodes of Gauss-Legendre integration, belong to the m * ×1 space, denoted as where m * is the number of unit Gauss integration points, represents the r-th unit Gauss integration point in the e-th unit; the weight of the unit Gauss integration point belongs to the m * ×1 space, denoted as is the weight corresponding to the r-th unit Gauss integration point in the e-th unit; the Jacobi determinant value J of the unit Gauss integration point e belongs to the m * ×1 space, denoted as where J er is the Jacobi determinant value of the r-th unit Gauss integration point in the e-th unit.

4. A discrete method for KL expansion random field of shape function considering irregular domain of slope soil mass according to claim 1, characterized in that The number of element shape functions N* in Step 4 is related to the shape function order P, specifically as follows: When P = 1, the element shape function N is a first-order shape function, N* = 4, which is consistent with the number of nodes of the shape function mesh element, and is defined in the local coordinate space Ω std = [-1, 1], and the expressions of the four first-order shape functions are respectively: where (ξ, η) are the local coordinate spaces Ω respectively std = coordinates within [-1, 1], corresponding to (x, y) in the global coordinate space. N1(ξ, η), N2(ξ, η), N3(ξ, η), and N4(ξ, η) are the 1st, 2nd, 3rd, and 4th shape functions within the element shape function N respectively; When P = 2, the element shape function N adds 4 second-order shape functions on the basis of the previous order, i.e., N * = 8, where the second-order shape functions are constructed by one-dimensional shape functions and Legendre polynomials. The expressions of the 4 second-order shape functions are respectively: In the formula, L0(ξ) and L2(ξ) are the values of the 0th-order and 2nd-order Legendre polynomials at ξ respectively; L0(-ξ) and L2(-ξ) are the values of the 0th-order and 2nd-order Legendre polynomials at -ξ respectively; L0(η) and L2(η) are the values of the 0th-order and 2nd-order Legendre polynomials at η respectively; L0(-η) and L2(-η) are the values of the 0th-order and 2nd-order Legendre polynomials at -η respectively; N5(ξ, η), N6(ξ, η), N7(ξ, η), and N8(ξ, η) are the 5th, 6th, 7th, and 8th shape functions in the element shape function N respectively; The Legendre polynomials are realized by calling the legendreP function in the MATLAB environment; When P = 3, the element shape function N adds 4 third-order shape functions on the basis of the previous order, that is, N* = 12. Among them, the expressions of the 4 third-order shape functions are: where \(L_1(\xi)\) and \(L_3(\xi)\) are the values of the first - order and third - order Legendre polynomials at \(\xi\) respectively; \(L_1(-\xi)\) and \(L_3(-\xi)\) are the values of the first - order and third - order Legendre polynomials at \(-\xi\) respectively; \(L_1(\eta)\) and \(L_3(\eta)\) are the values of the first - order and third - order Legendre polynomials at \(\eta\) respectively; \(L_1(-\eta)\) and \(L_3(-\eta)\) are the values of the first - order and third - order Legendre polynomials at \(-\eta\) respectively; \(N_9(\xi,\eta)\), \(N\) 10 (\xi,\eta)\), \(N\) 11 (\xi,\eta)\), \(N\) 12 (\xi,\eta)\) are the 9th, 10th, 11th, and 12th shape functions within the element shape function \(N\) respectively; When P = 4, the element shape function N adds 5 fourth-order shape functions on the basis of the previous order, that is, N * = 17, where the expressions of the 5 fourth-order shape functions are respectively: where \(L_2(\xi)\) and \(L_4(\xi)\) are the values of the second - order and fourth - order Legendre polynomials at \(\xi\) respectively, \(L_2(-\xi)\) and \(L_4(-\xi)\) are the values of the second - order and fourth - order Legendre polynomials at \(-\xi\) respectively, \(L_2(\eta)\) and \(L_4(\eta)\) are the values of the second - order and fourth - order Legendre polynomials at \(\eta\) respectively, \(L_2(-\eta)\) and \(L_4(-\eta)\) are the values of the second - order and fourth - order Legendre polynomials at \(-\eta\) respectively; \(N 13 (\xi,\eta)\), \(N 14 (\xi,\eta)\), \(N 15 (\xi,\eta)\), \(N 16 (\xi,\eta)\), \(N 17 (\xi,\eta)\) are the 13th, 14th, 15th, 16th, and 17th shape functions within the element shape function \(N\) respectively.

5. A discrete method for KL expansion random field of shape function considering irregular domain of slope soil mass according to claim 1, characterized in that The global stiffness submatrix K described in Step 5 is assembled by arranging the element stiffness submatrices k e The elements inside are assembled according to the shape function index IEN and belong to M * ×N s space, denoted as The correspondence rule between the elements inside the global stiffness submatrix K and the elements inside k e is as follows: i = (e - 1)m* + i* In the formula, are the elements of the i-th row and j-th column in the global stiffness sub-matrix K for the e-th and (e - 1)-th iterations, respectively, where i is determined by the i-th * element Gaussian integration point within the e-th shape function mesh element; j is determined by the shape function index IEN corresponding to the j-th * shape function in the e-th shape function mesh element; repeat this process until all the elements within the element stiffness matrix k e are assembled into the global stiffness sub-matrix K; The expression of the autocorrelation matrix ρ is as follows: where sym. represents a symmetric matrix; represents the t-th all-element Gauss integration point and the t'-th all-element Gauss integration point The autocorrelation function value obtained through the autocorrelation function.

6. A discrete method for KL expansion random field of shape function considering irregular domain of slope soil mass according to claim 1, characterized in that The global mass matrix E described in step 6 g is the element mass matrix E e The elements inside are assembled according to the shape function index IEN and belong to N s ×N s space, denoted as The corresponding assembly rule can be written as the following formula: In the formula, are the global mass matrices E for the e-th and (e - 1)-th iterations respectively g The elements in the ii-th row and jj-th column within, where ii and jj are determined by the shape function indices IEN corresponding to the ii-th * th and jj-th * shape functions in the e-th shape function element; repeat this process until all elements within the element mass matrix E e are assembled into the global mass matrix E g ; The feature function matrix D u is calculated as follows: Introduce the discrete point element eigenfunction matrix D ue , belonging to E * ×N * space, is the matrix composed of the values of the element shape function N at the local coordinates of the center point of the discrete grid element ; Assemble D ue into the eigenfunction matrix D u , and the column elements in D u are the eigenfunctions corresponding to the eigenvalues λ k , denoted as D uk .

Images