A non-invasive stochastic seepage analysis method based on two-dimensional finite element

Abstract

The present invention provides a non-intrusive stochastic seepage analysis method based on two-dimensional finite elements, belonging to the technical field of random field simulation in geotechnical engineering. The characteristics are as follows: 1) Establish a finite element model for deterministic unsaturated seepage analysis; 2) Based on the midpoint method, the statistical characteristics of the hydraulic parameters of the materials in the finite element model, and the spatial coordinate information of the grids, perform random field discretization on the hydraulic parameters of different materials; 3) Perform batch non-intrusive calculations on the discretely obtained stochastic seepage model; 4) Perform batch post-processing to obtain the seepage control indexes in the calculation results; 5) Obtain the statistical characteristics of the stochastic seepage in geotechnical engineering and complete the safety status assessment. The present invention can solve the problem of stochastic seepage analysis of the spatial variability of hydraulic parameters of different materials in complex geotechnical engineering. By performing random field discretization on the hydraulic parameters of multiple materials in complex geotechnical engineering and completing the coupled analysis of multiple random parameters, the seepage control indexes in complex geotechnical engineering can be obtained more reasonably.

Description

Technical Field

[0001] The present invention belongs to the field of stochastic seepage analysis in geotechnical engineering, and relates to a non-intrusive stochastic seepage analysis method based on two-dimensional finite elements. Technical Background

[0002] Seepage prevention and seepage control are important guarantees for the safe and stable operation of complex geotechnical engineering. However, due to the influence of geological actions such as material differences, sedimentation, and other uncontrollable artificial factors, the hydraulic parameters of the seepage prevention structure have strong variability in spatial distribution, and this spatial variability poses challenges to the seepage safety of complex geotechnical engineering.

[0003] The research on the spatial variability of soil parameters mainly focuses on aspects such as the statistical distribution law of material parameters, the spatial discretization method of random fields, and stochastic numerical simulation. First, from a statistical perspective, the spatial distribution law of soil parameters can be described by Gaussian and non-Gaussian distributions. Second, in order to perform stochastic finite element calculations, it is necessary to discretize the random field into a set of random variables that can be coupled with the finite element mesh. Finally, combining the random field of material parameters with deterministic analysis is an important means to study the influence of the spatial variability of material parameters on structural safety, which mainly includes two modes: intrusive and non-intrusive. The intrusive stochastic analysis method is limited due to the difficulties in establishing and solving equations in complex non-linear stochastic systems, while the non-intrusive stochastic analysis decouples the random field analysis from the deterministic numerical analysis, and has fewer limitations and is more acceptable, thus being widely applied. Currently, there is a lack of research on stochastic seepage considering the coupling of multiple materials and multiple hydraulic parameters in complex geotechnical engineering. Therefore, it is necessary to establish a non-intrusive stochastic seepage analysis method suitable for complex geotechnical engineering. Summary of the Invention

[0004] To solve the problems existing in the above-mentioned prior art, the present invention provides a non-intrusive stochastic seepage analysis method based on two-dimensional finite elements. On the basis of discretizing the hydraulic parameters of different materials into random fields by the midpoint method, a non-intrusive stochastic seepage analysis method is established, providing a feasible solution for obtaining the stochastic seepage characteristics of multiple materials and multiple hydraulic parameter couplings in complex geotechnical engineering, effectively obtaining the stochastic seepage characteristics and seepage safety status of complex geotechnical engineering, and providing guidance for the seepage safety control of complex geotechnical engineering.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is as follows:

[0006] A non-intrusive stochastic seepage analysis method based on two-dimensional finite elements, comprising the following steps:

[0007] The first step is to establish a deterministic unsaturated seepage analysis model

[0008] S1: Establish a finite element model according to geological data;

[0009] S2: Subdivide the finite element model into saturated and unsaturated zones;

[0010] S3: Determine the unsaturated parameters;

[0011] (1) Unsaturated seepage materials have strong fine-wet characteristics, which are described by the soil-water characteristic curve. The commonly used soil-water characteristic curve is the Van Genuchten model;

[0012] (2) For the permeability coefficient K, the material in the unsaturated seepage state is also a function of the saturation degree, as shown in Equation (1);

[0013] K = K SAT ·S e μ (1)

[0014] where K SAT is the permeability coefficient at saturation, μ is the empirical coefficient, and S e is the effective saturation degree.

[0015] S4: Determine the materials that require random seepage analysis;

[0016] The second step is the discretization of the non-Gaussian random field based on the midpoint method;

[0017] S5: Determine the distribution type of the parameters;

[0018] For the hydraulic parameters (permeability coefficient K, Van Genuchten model parameters) in different materials, the random field discretization needs to be carried out according to their statistical characteristics, which generally includes distribution types such as normal distribution, lognormal distribution, and Weibull distribution. In this patent, the lognormal distribution is mainly used.

[0019] S6: Select the autocorrelation function;

[0020] S7: Carry out the midpoint method non-Gaussian random field discretization process according to the autocorrelation function in S6;

[0021] S8: Random field discretization based on the deterministic analysis model;

[0022] The third step is the non-intrusive random seepage analysis;

[0023] S9: Carry out non-intrusive random seepage calculation to obtain the calculation results of flow velocity, pore pressure, and water head;

[0024] The fourth step is post-processing to obtain the seepage safety control index and safety status assessment;

[0025] S10: Extract the seepage safety control index;

[0026] S11: Seepage safety status assessment.

[0027] Preferably, the specific method for establishing a finite element model according to geological data in step S1 is as follows:

[0028] (1) Using Abaqus CAE, establish a two-dimensional solid model according to the geological data of the earth-rock dam and complete the finite element mesh division;

[0029] (2) In the load module of Abaqus CAE, set the initial water head boundary conditions for the upstream and downstream and the drainage boundary conditions on the downstream dam surface according to the data.

[0030] Preferably, the specific method for subdividing multiple groups of materials in the finite element model of complex geotechnical engineering into saturated and unsaturated zones and determining the unsaturated seepage control equation in step S2 is as follows;

[0031] (1) The seepage of complex geotechnical engineering is mainly unsaturated seepage. The established finite element model needs to divide the materials in the project into saturated and unsaturated zones (for example, the core wall, transition layer, and rockfill body in the core wall rockfill dam are usually unsaturated zones, and the bedrock and curtain grouting in the dam foundation are saturated zones);

[0032] (2) When giving attributes to unsaturated materials, the relationship between pore pressure and saturation needs to be given in the property module of Abaqus CAE. This relationship can be obtained through the soil-water characteristic curve (SWCC curve). The frequently used Van Genuchten model is adopted in the present invention, and the expression is as shown in (2)

[0033]

[0034] Among them, is the volumetric water content function of soil, is the matric suction, S e is the saturation degree, a represents the suction corresponding to the inflection point of the soil-water characteristic curve, n represents the slope at the inflection point of the soil-water characteristic curve, m = 1 - 1 / n, and the matric suction and pore pressure can be regarded as opposite numbers;

[0035] (3) For the unsaturated materials in the dam body, use equation (3) to represent the seepage control equation for unsaturated seepage;

[0036]

[0037] Among them, is the water head, K X is the horizontal permeability coefficient, K y is the vertical permeability coefficient, and C is the water capacity.

[0038] (4) The finite element model obtained after partitioning the above materials, assigning material parameters, applying boundary conditions, and discretizing the finite element mesh is exported as an INP format file that conforms to the Abaqus CAE calculation format and contains the information of each element and node. The deterministic seepage analysis model is within this file, which is an important basis for the following stochastic seepage analysis.

[0039] Preferably, the specific method for determining the materials and parameters that require stochastic seepage analysis in step S4 is as follows:

[0040] (1) The permeability coefficient K is one of the parameters that need to be considered key in stochastic seepage analysis. However, the stochastic characteristics of the permeability coefficient K of each material in geotechnical engineering do not all need to be considered. The specific screening criteria are as follows: For materials with a large difference in permeability coefficient K (more than two orders of magnitude difference, such as the permeability coefficient of the core wall of an earth-rock dam being 10 -8 m / s, which is much smaller than that of the rockfill of 10 -3 m / s), the influence of the non-seepage prevention material with a large permeability coefficient on stochastic seepage can be ignored. Therefore, only the stochastic seepage characteristics of the material with a small permeability coefficient need to be considered (such as all the anti-seepage structures in an earth-rock dam);

[0041] (2) The two parameters a and m in the Van Genuchten model are parameters for unsaturated stochastic seepage analysis, which determine the seepage characteristics of the material. For example, the core wall in a core wall dam, as an important anti-seepage material, mainly operates in an unsaturated seepage mode. The soil-water characteristic curve of the core wall is an important influencing factor for unsaturated seepage characteristic analysis and safety state evaluation. Therefore, it is necessary to consider the influence of the parameters a and m of the core wall on stochastic seepage.

[0042] Preferably, the specific method for selecting the autocorrelation function in step S6 is as follows:

[0043] When performing random field simulation calculations, the material parameters at different positions have autocorrelation, which can be described by the autocorrelation function. The autocorrelation function should be considered and selected according to factors such as the smoothness, continuity, and discrete error of the parameter random field distribution. It includes exponential type, Gaussian type, triangular type, exponential cosine type, etc. The present invention only shows the commonly used exponential type autocorrelation function;

[0044]

[0045] Among them, ρ is the autocorrelation coefficient, Z is the position coordinate of the discrete space point of the random field, τ x is the horizontal spatial distance between any two points, τ y is the vertical spatial distance between any two points, l h is the autocorrelation distance in the horizontal direction of the random field, l vis the autocorrelation distance in the vertical direction of the random field; the autocorrelation distance reflects the degree of spatial autocorrelation of material parameters. The larger the autocorrelation distance, the stronger the autocorrelation degree of the parameters.

[0046] Preferably, the specific method of the midpoint method non-Gaussian random field discretization process in step S7 is as follows:

[0047] The generation of the non-Gaussian random field is generated from the standard Gaussian random field through a series of transformations. The specific steps of a single random field are as follows:

[0048] (1) Select the material in the INP file of the finite element model in step S2 that needs to be discretized by the random field. According to the number of nodes in the finite element mesh of the selected material, generate independent standard normal random samples by means of Latin hypercube sampling to form an independent standard normal random sample matrix δ;

[0049] (2) According to the coordinate information of the finite element mesh of the selected material in the INP file of the finite element model in step S2, combined with the autocorrelation function selected in step S6, form the autocorrelation coefficient matrix R of the material parameters at any two points in the random field region;

[0050] (3) Perform Cholesky decomposition on the matrix R according to equation (5) to obtain the lower triangular matrix L, and then obtain the standard Gaussian random field F according to equation (6);

[0051] L L T =R (5)

[0052] F=Lδ (6)

[0053] (4) Convert the standard Gaussian random field F into a non-Gaussian random field P through equation (7)

[0054]

[0055] where G i -1 (·) is the inverse function of the non-Gaussian distribution marginal cumulative distribution, is the cumulative distribution function of the standard normal distribution;

[0056] (5) The distribution types of the parameter random fields K, a, and m to be analyzed are all lognormal distributions. By performing the transformation of equation (8) on the standard Gaussian field F, the lognormal random field P of the parameters can be obtained i ;

[0057] P i =exp(μ 1 +σ 1 ·F) (8)

[0058] In the formula, μ 1 =lnμ 2 -σ2 2 / 2; μ 1 is the mean value of the parameters (K, a, m) of the lognormal variable, and σ 1 is the standard deviation of the parameters of the lognormal variable, and μ 2 is the mean value of the normal variables (lnK, lna, lnm), and σ 2 is the standard deviation of the normal variable.

[0059] (6) Repeat the above (1) to (5) until the random field discretization of all materials is completed.

[0060] Preferably, step S8 combines the deterministic analysis model obtained in S1 to S3 with the random field obtained in S7 to obtain a non-intrusive random seepage calculation file, and the specific settings are as follows:

[0061] (1) Generally, in order to fully consider the influence of material spatial variability, several random fields need to be discretized for each working condition to obtain stable index distribution characteristics. And to avoid excessive reduction in calculation efficiency due to the number of random fields, in the present invention, the process of S7 is repeated 100 times as a standard, 100 random fields are discretized, and a batch calculation file is generated for subsequent calculation and safety evaluation;

[0062] (2) According to the random field discretization results of different materials, on the basis of the deterministic analysis model, the material units in the deterministic analysis model are replaced with random material units and saved in a file format that meets the Abaqus calculation. Finally, 100 non-intrusive random seepage calculation files are obtained for each working condition;

[0063] Preferably, in step S9, the specific method for performing random seepage calculation based on the non-intrusive random seepage calculation file obtained in S8 to obtain the calculation results of flow velocity, pore pressure, and water head is as follows:

[0064] (1) By calling the calculation function of Abaqus, batch calculate the non-intrusive random seepage calculation files obtained in step S8 to generate a calculation result file;

[0065] (2) According to the calculation result file, batch obtain the number information of the units and nodes, the spatial coordinate information of each node, the pore pressure, the pressure water head information, and the flow velocity information of each unit in the random field calculation results;

[0066] (3) Among them, the water head needs to be further calculated from the pore pressure and the node coordinates, and can be obtained by formula (9);

[0067]

[0068] Where Z is the water head, P is the pore pressure, ρ is the density of water, and Y is the coordinate value of the node in the Y direction.

[0069] Preferably, the specific method for extracting the seepage safety control index in step S10 is as follows:

[0070] The seepage safety control index includes the phreatic line, hydraulic gradient, and sectional flow rate;

[0071] For the extraction of the phreatic line, the distribution of the phreatic line within the element is determined based on the pore pressure of the element nodes. The specific steps are as follows;

[0072] (1) According to the calculation result information of the random field pore pressure obtained in step S9, first screen out the elements that intersect the phreatic line, which are called interface elements in the present invention. Assume that the number of nodes with pore pressure P < 0 among the four nodes on the quadrilateral element in the random field is D, and 0 > D > 4 is the interface element, and the phreatic line is distributed within these elements;

[0073] (2) Find the intersection lines of the four sides of the interface element and the phreatic line. According to the pore pressure of the nodes within the interface element, when the pore pressures P of two nodes on a certain side are not both positive or both negative at the same time, there is an intersection point between this side and the phreatic line. Each interface element has two intersection points, and thus two sides intersecting the phreatic line can be found;

[0074] (3) After finding the two sides with intersection points of the interface element in (2), use the linear proportion method to find the intersection points. The intersection point coordinates can be obtained from equations (10) and (11);

[0075]

[0076]

[0077] Where X 0 , Y 0 are the horizontal and vertical coordinates of the intersection point of the interface element and the phreatic line, X 1 , X 2 , Y 1 , Y 2 are the horizontal and vertical coordinates of the two nodes on the side where the intersection point of the interface element and the phreatic line is located, and P 1 , P 2 are the pore pressures of the two nodes on the side where the intersection point of the interface element and the phreatic line is located;

[0078] (4) Repeat (1) to (3) for all elements of the finite element model to find all interface elements and the intersection points within the interface elements, and connect the two intersection points within the interface element in a straight line to obtain the distribution of the phreatic line for the entire project.

[0079] For the extraction of the gradient, the gradient is obtained according to Darcy's law of seepage, and the hydraulic gradient is obtained according to equation (12);

[0080]

[0081] Among them, J is the hydraulic gradient, V is the total flow velocity, and K is the permeability coefficient

[0082] Flow rate acquisition: The cross-sectional flow rate is obtained based on the flow velocity of the unit on the cross-section, and the cross-sectional flow rate of the required cross-section is obtained according to Equation (13);

[0083] Q = V·A (13)

[0084] Among them, Q is the cross-sectional flow rate, V is the total flow velocity, A is the cross-sectional area, and in the two-dimensional model, A = Δy, where Δy is the ordinate difference of the required interface;

[0085] Preferably, the specific steps for evaluating the seepage safety state in step S11 are as follows:

[0086] (1) Based on the phreatic line obtained in step S10, the distribution range of the phreatic line can be drawn, and a reasonable distribution range of the interface between the saturated area and the unsaturated area and the overflow point can be given, providing a target range for seepage control and reinforcement;

[0087] (2) Analyze the influence of the variation coefficient values of the hydraulic parameters in different materials on the seepage control index, and obtain the materials that are vulnerable to spatial variability, providing a reference for subsequent anti-seepage;

[0088] (3) Analyze the influence of the change in the variation coefficient on the seepage safety control index, and obtain the range of the variation coefficient that makes the structure in a relatively safe state.

[0089] (4) Based on the cross-sectional flow rate and hydraulic gradient obtained in step S10, the probability density function PDF of different variation coefficients can be obtained, and combined with the cumulative probability distribution function CDF of the hydraulic gradient, the probability of seepage failure of this material can be given (including but not limited to directly obtaining the failure probability by the Monte Carlo simulation method).

[0090] Advantages of the present invention:

[0091] The present invention provides a non-intrusive random seepage analysis method based on two-dimensional finite elements for analyzing the random seepage characteristics of complex geotechnical engineering. On the basis of discretizing the hydraulic parameters of different materials by the midpoint method in the random field, a non-intrusive random seepage analysis method is established, providing a feasible solution for obtaining the random seepage characteristics of multi-material and multi-hydraulic parameter coupling in complex geotechnical engineering, effectively obtaining the random seepage characteristics and seepage safety state of complex geotechnical engineering, and providing guidance for the seepage safety control of complex geotechnical engineering. Description of the drawings

[0092] Figure 1 It is a schematic diagram of the non-intrusive random seepage analysis process;

[0093] Figure 2 It is a schematic diagram of material zoning for an ultra-high core rockfill dam;

[0094] Figure 3 It is a schematic diagram of a two-dimensional finite element model for an ultra-high core rockfill dam;

[0095] Figure 4 It is a schematic diagram of the permeability coefficient and hydraulic gradient of a homogeneous field; Figure 4 (a) is the permeability coefficient diagram of the homogeneous field; Figure 4 (b) is the hydraulic gradient diagram of the homogeneous field;

[0096] Figure 5 It is the coefficient of variation COV K = 0.6 of the schematic diagram of the hydraulic gradient of the permeability coefficient random field; Figure 5 (a) is the permeability coefficient diagram of the random field; Figure 5 (b) is the hydraulic gradient diagram of the random field;

[0097] Figure 6 It is the distribution range of the core wall seepage line (COV K = 0.6) schematic diagram; Figure 6 (a) is the distribution diagram of the core wall seepage line; Figure 6 (b) is the distribution diagram of the seepage line in the middle area of the core wall; Figure 6 (c) is the distribution of the core wall near the overflow point;

[0098] Figure 7 It is a schematic diagram of the distribution range of the core wall seepage line for random fields with different hydraulic parameters. Figure 7 (a) is the coefficient of variation COV K = 0.6, COV a = 0.15, COV n = 0.1 of the distribution range of the core wall seepage line; Figure 7 (b) is the coefficient of variation COV K = 0.6, COV a = 0.2, COV n = 0.15 of the seepage line distribution range; Figure 7 (c) is the coefficient of variation COV a = 0.2, COV n = 0.15 of the distribution range of the core wall seepage line; Figure 7 (d) is the coefficient of variation COV a = 0.3, COV n = 0.15 of the distribution range of the core wall seepage line;

[0099] Figure 8 It is a schematic diagram of the results of random seepage analysis of different materials. Figure 8 (a1) Coefficient of variation of the core wall COVK = 0.6, COV K = 0.3, COV K = 0.1 when the maximum hydraulic gradient PDF diagram; Figure 8 (a2) is the coefficient of variation COV of the core wall K = 0.6, COV K = 0.3, COV K = 0.1 when the PDF diagram of the cross-sectional flow rate; Figure 8 (a3) is the diagram of the variation law of the mean value of the maximum hydraulic gradient and the mean value of the flow rate of the core wall with the coefficient of variation; Figure 8 (b1) is the coefficient of variation COV of the curtain grouting K = 1, COV K = 0.5, COV K = 0.1 when the maximum hydraulic gradient PDF diagram; Figure 8 (b2) The coefficient of variation COV of the curtain grouting K = 1, COV K = 0.5, COV K = 0.1 when the PDF diagram of the cross-sectional flow rate; Figure 8 (b3) Diagram of the variation law of the mean value of the maximum hydraulic gradient and the mean value of the flow rate of the curtain grouting with the coefficient of variation; Figure 8 (c1) The coefficient of variation COV of the slightly weathered rock mass K = 1, COV K = 0.5, COV K = 0.1 when the maximum hydraulic gradient PDF diagram; Figure 8 (c2) is the coefficient of variation COV of the slightly weathered rock mass K = 1, COV K = 0.5, COV K = 0.1 when the PDF diagram of the cross-sectional flow rate; Figure 8 (c3) is the diagram of the variation law of the mean value of the maximum hydraulic gradient and the mean value of the flow rate of the slightly weathered rock mass with the coefficient of variation;

[0100] Figure 9 Schematic diagram of the coupling analysis results of the random seepage of the bedrock and the curtain grouting. Figure 9 (a) is the coefficient of variation COV K = 1, 0.5, 0.1 when the probability density distribution of the maximum hydraulic gradient of the coupling of the curtain grouting and the bedrock-curtain grouting; Figure 9 (b) is the coefficient of variation COV K = 1, 0.5, 0.1 when the probability density distribution of the seepage flow rate drop of the coupling of the curtain grouting and the bedrock-curtain grouting;

[0101] Figure 10It is a comparison diagram of a homogeneous field without considering spatial variability, single-material random seepage, and multi-material coupled random seepage. Figure 10(a) shows the multi-material coupling of curtain grouting, weakly weathered zone, and slightly weathered rock mass, and the coefficient of variation of single-material curtain grouting is COV K = 0.1, and the upstream side pressure distribution diagram of curtain grouting in a homogeneous field without considering spatial variability; Figure 10 (b) shows the multi-material coupling of curtain grouting, weakly weathered zone, and slightly weathered rock mass, and the coefficient of variation of single-material curtain grouting is COV K = 1.0, and the upstream side pressure distribution diagram of curtain grouting in a homogeneous field without considering spatial variability;

[0102] Figure 11 It is the probability density distribution function of the maximum hydraulic gradient when the permeability coefficient of the core wall is at different coefficients of variation;

[0103] Figure 12 It is a schematic diagram of the finite element model of a two-dimensional homogeneous dam;

[0104] Figure 13 It is the contour cloud diagram of pore pressure of a two-dimensional homogeneous dam;

[0105] Figure 14 It is a schematic diagram of the distribution range of the phreatic line of different hydraulic parameters of a two-dimensional homogeneous dam with a random field; Figure 14 (a) shows the distribution range of the phreatic line when the coefficient of variation of the dam body COV K = 0.6, COV a = 0.4, COV n = 0.1; Figure 14 (b) shows the distribution range of the phreatic line when the coefficient of variation of the dam body COV a = 0.4, COV n = 0.1;

[0106] Figure 15 It is a schematic diagram of the hydraulic gradient of a homogeneous field and a random field. Figure 15 (a) is a schematic diagram of the hydraulic gradient of a homogeneous field; Figure 15 (b) is a schematic diagram of the hydraulic gradient of a random field when the coefficient of variation of the dam body COV K = 0.6, COV a = 0.4, COV n = 0.1;

[0107] Figure 16 It is a schematic diagram of the permeability coefficient K, SWCC curve parameters a and n of a random field. Figure 16 (a) is a schematic diagram of the permeability coefficient K of the dam body when COV K = 0.6; Figure 16 (b) is a schematic diagram of the model parameter a when COV a = 0.4; Figure 16(c) is COV n Schematic diagram of model parameter n when COV = 0.1;

[0108] Figure 17 is COV K = 0.6, COV a = 0.4, COV n = 0.1, probability density function (PDF) of the maximum hydraulic gradient and discharge of the dam body. Figure 17 (a) is COV K = 0.6, COV a = 0.4, COV n = 0.1, PDF diagram of the maximum hydraulic gradient of the dam body; Figure 17 (b) is COV K = 0.6, COV a = 0.4, COV n = 0.1, PDF diagram of the hydraulic gradient;

[0109] In the figure, COV K , COV a , COV n respectively represent the coefficient of variation of the permeability coefficient K, the SWCC curve parameters a and n, and PDF is the probability density function. Specific implementation mode

[0110] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be described in detail and clearly below with reference to the accompanying drawings in the embodiments of the present invention. The specific embodiments described are only a part of the numerous embodiments of the present invention and do not represent all embodiments. Based on the specific embodiments described in the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the scope of protection of the present invention.

[0111] Embodiment 1

[0112] The present invention provides a non-invasive stochastic seepage analysis method based on two-dimensional finite elements: used to cooperate with software Abaqus to generate a hydraulic parameter random field and perform seepage calculation work, and established a non-invasive stochastic seepage analysis method for high core wall dams based on the interaction between Abaqus and Matlab. Specifically, the basic process can be as Figure 1 shown, including:

[0113] In step S1, a finite element model is established based on geological data. The two-dimensional core wall rockfill dam model adopted in the present invention has a dam height of 315 m. The dam building materials include rockfill, transition material, filter material, core wall, contact clay, etc. The dam foundation consists of curtain grouting, weakly unloaded, weakly weathered upper zone, weakly weathered lower zone and slightly weathered rock mass. Among them, the depth of curtain grouting is 126 m, which is grouted into the slightly weathered rock mass. The specific distribution pattern is shown in Figure 2 , and the two-dimensional finite element model diagram is shown in Figure 3 . The initial water head boundary of this model is that the upstream pore pressure boundary is 2798.3 kPa, the downstream pore pressure boundary is -2.9 kPa, and the downstream dam surface is a drainage boundary.

[0114] In step S2, the finite element model needs to be subdivided into saturated and unsaturated zones. Among them, the unsaturated zone can be divided into rockfill, filter material, transition material, and core wall, and the saturated zone can be divided into curtain grouting, weakly unloaded, upper and lower weakly weathered zones, and slightly weathered rock mass. In the unsaturated seepage calculation, the bedrock and curtain grouting materials in the saturated zone do not need to be considered, and the SWCC curve does not need to be considered; in the unsaturated materials, the core wall is used as the main anti-seepage structure and the SWCC curve needs to be considered. For the remaining fill materials, since the permeability coefficients are very different from those of the core wall, the SWCC curve does not need to be considered.

[0115] In step S3, the Van Genuchten model of the SWCC curve of the core wall is determined, and the following formula is obtained according to the test;

[0116]

[0117] In step 4, the materials for random seepage analysis are the core wall, curtain grouting, and foundation materials;

[0118] In step S5, the distribution type of the parameters is determined, and the statistical characteristics of the parameters of each material are shown in Table 1;

[0119] Table 1

[0120]

[0121]

[0122] In step S6, an autocorrelation function is selected. In the present invention, an exponential autocorrelation function is adopted, and then a non-Gaussian random field based on the midpoint method is generated according to step S7;

[0123] In step S8, the elements in the deterministic analysis model need to be replaced with random field elements, and a random seepage analysis INP file that meets the Abaqus calculation format is formed. In the present invention, 100 random fields are discretized to reduce the influence of the number of random fields;

[0124] Batch obtain the seepage safety control indicators according to steps S9 and S10, and finally conduct a seepage safety status assessment according to this indicator in step S11. As Figure 4 shown are the core wall permeability coefficient K and hydraulic gradient without considering spatial variability, Figure 5 shown is the coefficient of variation COV K = 0.6, the permeability coefficient and hydraulic gradient. By comparing the two figures, it is found that the spatial variability of the permeability coefficient has a greater impact on the hydraulic gradient. Figure 6 shown is COV K = 0.6, the distribution range of the phreatic line, Figure 7 is a schematic diagram of the distribution range of the phreatic line of the core wall with random fields of different hydraulic parameters. By comparing Figure 6 with Figure 7 it can be found that the spatial variability of the soil-water characteristic curve parameters is not much different from the permeability coefficient, and the phreatic line of the core wall is mainly distributed near the boundary on the downstream side of the core wall. Therefore, it is necessary to focus on this area. Figure 8 shown is a schematic diagram of the results considering the spatial variability of hydraulic parameters of different materials. When only considering the spatial variability of a single material, the PDF of the maximum hydraulic gradient in the core wall basically satisfies the lognormal distribution, and the value of the coefficient of variation of the permeability coefficient in the curtain grouting has little influence on the maximum hydraulic gradient and cross-sectional flow rate. Figure 9 is the analysis result of single-material curtain grouting and multi-material bedrock-curtain grouting coupled random seepage. The hydraulic gradient and cross-sectional flow rate of multi-material coupling are both greater than those of single-material. Figure 10 is the pressure head distribution diagram on the upstream side of the curtain grouting of single-material curtain grouting random seepage, multi-material curtain grouting-bedrock material coupled random seepage, and homogeneous field without random simulation. Figure 10 (a) is COV K = 0.1, Figure 10 (b) is COV K = 1.0, the analysis results. By comparison, it is found that the pressure head range of multi-material coupling is greater than that of single-material and homogeneous field. Therefore, for complex geotechnical engineering, it is necessary to consider its random seepage characteristics of multi-material-multi-hydraulic parameter random characteristic coupling. Figure 11 gives the cumulative distribution probability function of the maximum hydraulic gradient of the core wall material under different coefficients of variation. According to this figure, it can be found that when the coefficient of variation Cov takes 0.1 and 0.3, the maximum hydraulic gradient of the core wall is always less than the threshold of the critical hydraulic gradient. Therefore, it can be considered that the safety probability at this time is 100%, that is, the failure probability is infinitely close to 0%. However, when Cov takes 0.6, the probability that the maximum hydraulic gradient of the core wall is less than the critical hydraulic gradient 4 is 64%. It can be considered that the safety probability at this time is 64%, that is, the failure probability is infinitely close to 36%.

[0125] Example 2

[0126] Step S1 establishes a finite element model based on the earth-rock dam data. Figure 12 The schematic diagram of the finite element model of the two-dimensional homogeneous dam is shown. The dam height is 14m, the top width of the dam is 4m, the slope ratio is 1:3, the upstream water level is 13m, the thickness of the dam foundation is 10m, the permeability coefficient of the dam body is 1.6×10 -5 m / s, and the permeability coefficient of the dam foundation is 2.1×10 -4 m / s. The upstream pore pressure boundary is 130Kpa, the downstream pore pressure boundary is 0Kpa, and the downstream dam surface is the drainage boundary.

[0127] In step S2, it is necessary to subdivide the saturated and unsaturated zones of the finite element model, where the unsaturated zone is the dam body material and the saturated zone is the dam foundation material. In the unsaturated calculation, the dam foundation material in the saturated zone does not need to be considered, and the soil-water characteristic curve does not need to be considered; the soil-water characteristic curve needs to be considered for the unsaturated material dam body.

[0128] In step S3, the soil-water characteristic curve of the dam body needs to be determined, and the same curve as in Example 1 is adopted.

[0129] In step S4, the material for random seepage analysis is the dam body material.

[0130] In step S5, the distribution type of the dam body material parameters is determined, and the specific statistical characteristics are shown in Table 2.

[0131] Table 2

[0132]

[0133] In step S6, the autocorrelation function is selected. In the present invention, the exponential autocorrelation function is adopted, and then a non-Gaussian random field based on the midpoint method is generated according to step S7.

[0134] In step S8, the elements in the deterministic analysis model need to be replaced with random field elements, and a random seepage analysis INP file that meets the Abaqus calculation format is formed. In the present invention, 100 random fields are discretized to reduce the influence of the number of random fields.

[0135] According to steps S9 and S10, the seepage safety control indexes are obtained in batches, and finally, the seepage safety state assessment is carried out according to this index in step S11. Figure 14 (a) shows the distribution range of the phreatic line when the coefficient of variation COV K =0.6, COV a =0.4, COV n =0.1 of the dam body; Figure 14 (b) shows the distribution range of the phreatic line when the coefficient of variation COV a =0.4, COV n =0.1 of the dam body. It can be obtained that the influence of the spatial variability of the SWCC curve model parameters a and n on the phreatic line is much smaller than that of the permeability coefficient K; comparisonFigure 15 (a) and Figure 15 (b) It can be found that after considering the spatial variability of the permeability coefficient K, model parameters a and n, the hydraulic gradient is greatly affected; from Figure 16 it can be seen the spatial distribution after the random field discretization of different hydraulic parameters K, a and n; from Figure 17 it can be seen that after the non-invasive random seepage analysis of the permeability coefficient K, model parameters a and n of the dam body material, the calculated maximum hydraulic gradient and discharge both satisfy the lognormal distribution.

[0136] The above-described embodiments only represent the implementation modes of the present invention, but should not be construed as limiting the scope of the patent for the present invention. It should be noted that for those skilled in the art, without departing from the concept of the present invention, several modifications and improvements can still be made, and these all belong to the protection scope of the present invention.

Claims

1. A non-intrusive stochastic seepage analysis method based on two-dimensional finite element, characterized in that, it includes the following steps: The first step is to establish a deterministic unsaturated seepage analysis model S1: Establish a finite element model according to geological data; S2: Subdivide the saturated and unsaturated zones of the finite element model; S3: Determine the unsaturated parameters; S4: Determine the materials that require stochastic seepage analysis; The second step is the non-Gaussian random field discretization based on the midpoint method; S5: Determine the distribution type of the parameters; S6: Select the autocorrelation function; S7: Perform the non-Gaussian random field discretization process based on the autocorrelation function in S6; S8: Random field discretization based on the deterministic analysis model; The third step is non-intrusive stochastic seepage analysis; S9: Perform non-intrusive stochastic seepage calculation to obtain the calculation results of flow velocity, pore pressure, and water head; The fourth step is post-processing to obtain seepage safety control indicators and safety status evaluation; S10: Extract seepage safety control indicators; S11: Evaluate the seepage safety status; In step S2, the saturated and unsaturated zones of multiple groups of materials in the finite element model of complex geotechnical engineering are subdivided, and the unsaturated seepage control equation is determined; the seepage of complex geotechnical engineering is mainly unsaturated seepage. The established finite element model needs to divide the materials in the project into saturated and unsaturated regions. The core wall, transition layer, and rockfill body in the core wall rockfill dam are unsaturated zones, while the bedrock and curtain grouting in the dam foundation are saturated zones; The specific method for extracting seepage safety control indicators in step S10 is: The seepage safety control indicators include the phreatic line, hydraulic gradient, and cross-sectional flow rate; For phreatic line extraction, determine the phreatic line distribution within the element according to the pore pressure at the element nodes. The specific steps are as follows; (1) According to the information of the pore pressure calculation results of the random field obtained in step S9, first screen out the elements that intersect the phreatic line, which are called interface elements. Assume that the number of nodes with pore pressure P < 0 among the four nodes of the quadrilateral element in the random field is D. When 0 > D > 4, it is an interface element, and the phreatic line is distributed within these elements; (2) Find the intersection lines of the four sides of the interface element and the phreatic line; according to the pore pressure of the nodes within the interface element, when the pore pressures P of two nodes on a side are not both positive or both negative at the same time, there is an intersection point between this side and the phreatic line. Each interface element has two intersection points, and two sides that intersect the phreatic line can be found; (3) After finding the two sides with intersection points of the interface element in (2), use the linear proportion method to find the intersection points. The intersection point coordinates can be obtained from equations (10) and (11); Among them, X 0 , Y 0 are the horizontal and vertical coordinates of the intersection point of the interface unit and the phreatic line. X 1 , X 2 , Y 1 , Y 2 are the horizontal and vertical coordinates of the two nodes on the side where the intersection point of the interface unit and the phreatic line is located. P 1 , P 2 are the pore pressures of the two nodes on the side where the intersection point of the interface unit and the phreatic line is located; (4) Repeat (1) to (3) for all elements of the finite element model to find all interface elements and the intersection points within the interface elements. Connect the two intersection points within the interface element in a straight line to obtain the phreatic line distribution of the entire project; For gradient extraction, obtain the gradient according to Darcy's seepage law and obtain the hydraulic gradient according to equation (12); where J is the hydraulic gradient, V is the total flow velocity, and K is the permeability coefficient For flow rate acquisition, obtain the cross-sectional flow rate according to the flow velocity of the elements on the cross-section and obtain the cross-sectional flow rate of the required cross-section according to equation (13); Q = V·A (13) Where Q is the sectional flow rate, V is the total flow velocity, A is the sectional area, and in the two-dimensional model, A = Δy, where Δy is the required ordinate difference of the interface.

2. A non-intrusive stochastic seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, in the step S1, a finite element model is established as follows: (1) Using Abaqus CAE, a two-dimensional solid model is established according to the geological data of the earth-rock dam, and finite element mesh division is completed; (2) In the load module of Abaqus CAE, the initial head boundary conditions of the upstream and downstream and the drainage boundary conditions of the downstream dam surface are set according to the data.

3. A non-intrusive stochastic seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, the step S4 to determine the materials that need stochastic seepage analysis is as follows: (1) For materials with a difference in permeability coefficient K exceeding two orders of magnitude, the influence of the non-seepage-proof material with a large permeability coefficient on stochastic seepage is negligible, and only the stochastic seepage characteristics of the material with a small permeability coefficient are considered; (2) The two parameters a and m in the Van Genuchten model are parameters for unsaturated stochastic seepage analysis, which determine the seepage characteristics of the material.

4. A non-intrusive stochastic seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, the non-saturated parameter determination method of S3 is as follows: (1) The fine-wet characteristics of unsaturated seepage materials are described by the soil-water characteristic curve, and the soil-water characteristic curve is the Van Genuchten model; (2) For the permeability coefficient K, the material in the unsaturated seepage state is also a function of the degree of saturation, as shown in Equation (1); K = K SAT ·S e μ (1) Among them, K SAT is the permeability coefficient at saturation, μ is the empirical coefficient, and S e is the effective saturation.

5. A non-intrusive stochastic seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, step S2 further includes: (1) When giving attributes to unsaturated materials, the relationship between pore pressure and degree of saturation needs to be given in the Abaqus CAE attribute module, and this relationship can be obtained through the soil-water characteristic curve (SWCC curve), and the Van Genuchten model with the expression as (2) is adopted; Among them, is the volumetric water content function of soil, is the matric suction, S e is the degree of saturation, a represents the suction corresponding to the inflection point of the soil-water characteristic curve, n represents the slope at the inflection point of the soil-water characteristic curve, m = 1 - 1 / n, and the matric suction and pore pressure can be regarded as opposite numbers; (2) For the unsaturated materials in the dam body, the seepage control equation for unsaturated seepage is expressed by Equation (3); where is the water head, K X is the horizontal permeability coefficient, K y is the vertical permeability coefficient, and C is the water capacity; (3) The finite element model obtained after zoning the above materials, assigning material parameters, applying boundary conditions, and discretizing the finite element mesh is exported as an INP format file that conforms to the Abaqus CAE calculation format and contains the information of each element and node. This file is the deterministic seepage analysis model.

6. A non-intrusive stochastic seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, the specific method of the midpoint method non-Gaussian random field discretization process in step S7 is: The generation of the non-Gaussian random field is generated from the standard Gaussian random field through a series of transformations. The specific steps of a single random field are as follows: (1) Select the materials in the INP file of the finite element model in step S2 that need to be discretized by the random field. According to the number of nodes of the finite element mesh of the selected materials, generate independent standard normal random samples by means of Latin hypercube sampling to form an independent standard normal random sample matrix δ; (2) According to the coordinate information of the finite element mesh of the selected materials in the INP file of the finite element model in step S2, combined with the autocorrelation function selected in step S6, form an autocorrelation coefficient matrix R of the material parameters at any two points in the random field region; (3) Perform Cholesky decomposition on the matrix R according to formula (5) to obtain a lower triangular matrix L, and then obtain a standard Gaussian random field F according to formula (6); LL T = R(5) F = Lδ (6) (4) Convert the standard Gaussian random field F into a non-Gaussian random field P through formula (7) Among them, G i -1 (·) is the inverse function of the marginal cumulative distribution of the non-Gaussian distribution, is the cumulative distribution function of the standard quasi-normal distribution; (5) The distribution types of the parameter random fields K, a, and m to be analyzed are all lognormal distributions. By performing the transformation in Equation (8) on the standard Gaussian field F, the lognormal random field P of the parameters can be obtained. i ; P i = exp(μ 1 + σ 1 · F)(8) where μ 1 = ln μ 2 - σ 2 2 / 2; μ 1 is the mean of the parameters (K, a, m) of the lognormal variable, σ 1 is the standard deviation of the parameters of the lognormal variable, μ 2 is the mean of the normal variables (ln K, ln a, ln m), σ 2 is the standard deviation of the normal variable; (6) Repeat the above (1) to (5) until the random field discretization of all materials is completed.

7. A non-intrusive random seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, In step S8, the deterministic analysis model obtained in S1 to S3 is combined with the random field obtained in S7 to obtain a non-intrusive random seepage calculation file, which is specifically set as follows: (1) In order to fully consider the influence of material spatial variability, several random fields need to be discretized for each working condition, and then stable index distribution characteristics are obtained. Taking repeating the process of S7 100 times as the standard, 100 random fields are discretized, and a batch calculation file is generated for subsequent calculation and safety evaluation; (2) According to the random field discretization results of different materials, on the basis of the deterministic analysis model, replace the material units in the deterministic analysis model with random material units and save them in a file format that meets the Abaqus calculation. Finally, 100 non-intrusive random seepage calculation files are obtained for each working condition.

8. A non-intrusive random seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, The specific method for obtaining the calculation results of flow velocity, pore pressure, and water head by performing random seepage calculation on the non-intrusive random seepage calculation file obtained in step S9 based on S8 is as follows: (1) By calling the calculation function of Abaqus, batch calculate the non-intrusive random seepage calculation files obtained in step S8 to generate a calculation result file; (2) According to the calculation result file, batch obtain the number information of the units and nodes, the spatial coordinate information of each node, the pore pressure, the pressure head information, and the flow velocity information of each unit in the random field calculation results; (3) Among them, the water head needs to be further calculated from the pore pressure and the node coordinates, and can be obtained by formula (9); In the formula, Z is the water head, P is the pore pressure, ρ is the density of water, and Y is the coordinate value of the node in the Y direction.

9. A non-intrusive random seepage analysis method based on two-dimensional finite elements according to claim 1, characterized in that, The specific steps of the seepage safety state evaluation in step S11 are: (1) The phreatic line obtained according to step S10 can be used to draw the distribution range of the phreatic line, and to give the reasonable distribution ranges of the interface between the saturated area and the unsaturated area and the overflow point, providing a target range for seepage control and reinforcement; (2) Analyze the influence degree of the variation coefficient values of the hydraulic parameters in different materials on the seepage control index, and obtain the materials vulnerable to spatial variability, providing a reference for subsequent anti-seepage; (3) Analyze the influence of the change of the variation coefficient on the seepage safety control index, and obtain the range of the variation coefficient that makes the structure in a relatively safe state; (4) According to the sectional flow rate and hydraulic gradient obtained in step S10, the probability density function PDF of different variation coefficients can be obtained, and combined with the cumulative probability distribution function CDF of the hydraulic gradient, the probability of seepage failure of this material can be given.

Images