An analytical method for seepage problems with free surfaces based on the virtual element method

Abstract

This invention discloses an analysis method for seepage problems with free surfaces based on the virtual element method, comprising: Step 1: collecting observable data and using it as basic information to establish a seepage model of an earth-rock dam with free surfaces; Step 2: discretizing the seepage region to generate a fixed mesh; Step 3: calculating the element permeability matrix and aggregating it into a global permeability matrix, while simultaneously calculating the global flow vector; Step 4: applying boundary conditions and solving for the seepage free surface; Step 5: defining an error threshold; if the difference between the calculated result and the actual result is less than the threshold, proceed to Step 7; if the error check fails, proceed to Step 6; Step 6: using the generated free surface as a reference, locally modify the mesh, update the global permeability matrix and flow vector, and repeat Step 4; Step 7: outputting parameters within the seepage domain; This invention achieves automatic densification of nodes on the free surface during the iteration process, which is highly efficient and ensures good numerical accuracy.

Description

Technical Field

[0001] This invention relates to the field of seepage calculation and analysis technology in geotechnical engineering and hydraulic engineering, specifically to an analysis method for seepage problems with free surfaces based on the virtual element method. Background Technology

[0002] Free-face seepage is an important research topic frequently encountered in geotechnical engineering, hydraulic engineering, and other fields. The difficulty in its analysis lies in the fact that the location of the free surface is unknown beforehand, generally requiring iterative solutions. Finite element methods for solving this problem can be divided into two categories: the adjusted mesh method and the fixed mesh method. The adjusted mesh method treats the free surface as a movable boundary of the analysis domain, continuously adjusting the mesh shape and modifying the free surface position during iteration until the free surface position stabilizes. The adjusted mesh method suffers from significant drawbacks, such as high computational cost and the potential for mesh distortion or overlapping when the assumed free surface position differs greatly from the actual free surface position. Therefore, it is trending towards being completely replaced by the fixed mesh method. Nevertheless, due to its simplicity and intuitiveness, many scholars are still working to improve upon the adjusted mesh method.

[0003] Scholars both domestically and internationally have developed various fixed-mesh methods, which Zheng Hong et al. categorized into two types: intuitionistic methods and variational inequality methods. Intuitionistic methods include three approaches: adjusting flow rate, adjusting element permeation matrix, and locally modifying the mesh. Adjusting flow rate is similar to the initial stress method in stress analysis, requiring the calculation or consideration of the free surface. Adjusting the element permeation matrix is ​​similar to the tangent stiffness method, using the Heaviside function to extend Darcy's law, which is only applicable to wet regions, to the entire domain. For elements traversed by the free surface, artificial parameters are often introduced to fine-tune the Heaviside function to avoid oscillations in the solution. Locally modifying the mesh, such as the virtual element method, sub-element method, and element-dropping method, maintains a fixed mesh in form but modifies the mesh near the free surface during iteration; the area below the free surface is the effective region for analysis. However, all three approaches to intuitionistic methods cannot circumvent the problem of determining the free surface during the solution process.

[0004] The variational inequality method with a fixed mesh typically constructs a new problem within a fixed region, ensuring that the free surface and its conditions are not explicitly included in the problem. However, once this problem is solved, the free surface can be determined using certain functional conditions, such as the cutoff negative pressure method and the Signorini variational inequality formulation. It is worth noting that Zheng Hong et al. proposed Signorini-type conditions for the boundary conditions on the potential seepage surface, theoretically eliminating the singularity of the seepage point. The variational inequality method does not require attention to the location of the free surface during iterative solution, resulting in a more rigorous mathematical foundation, but its theory is relatively complex.

[0005] In addition, Li et al. attempted to solve the problem using a meshless method, but the biggest obstacle was the accurate application of the essential boundary conditions. Zheng et al. overcame this bottleneck by combining the meshless method with the manifold element method, and successfully solved the seepage problem with free surfaces, which is of great significance for the application of meshless methods in this problem.

[0006] From an engineering application perspective, while the traditional adjusted mesh method still has wide applications, it has many shortcomings; the variational inequality method of the fixed mesh method has not been widely used by engineers due to its theoretical complexity; and research on meshless methods for solving problems still needs to be enriched and deepened. Among the three intuitive fixed mesh methods, local mesh modification does not require the introduction of additional physical concepts or mathematical functions, making the approach more intuitive. However, due to the limitations imposed on element shapes by conventional finite element methods, careful handling of elements cut from free surfaces is often necessary when locally modifying the mesh.

[0007] The virtual element method (VEM), applicable to general polygonal meshes, constructs the matrices and vectors in the element equations through projection. Compared to the polygonal finite element method, it is computationally more convenient and faster. Its greatest advantage lies in its inclusiveness regarding element shapes. Polygonal elements formed by free surface cutting when performing fixed mesh methods with local mesh modification can be calculated and analyzed by VEM without additional processing. Therefore, combining VEM with the local mesh modification method to solve seepage problems has significant practical engineering implications. Summary of the Invention

[0008] The purpose of this invention is to overcome the above-mentioned shortcomings and provide an analysis method for seepage problems with free surfaces based on the virtual element method. This method overcomes the limitations of traditional local mesh modification and achieves automatic densification of nodes on the free surface during the iteration process, which is highly efficient and can ensure good numerical accuracy.

[0009] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: an analysis method for seepage problems with free surfaces based on the virtual element method, which includes the following steps:

[0010] Step 1: Collect observable data and use it as basic information to establish a seepage model of an earth-rock dam with free surfaces;

[0011] Step 2: Discretize the seepage region to generate a fixed grid;

[0012] Step 3: Calculate the unit permeability matrix and combine them into a global permeability matrix, and simultaneously calculate the global flow vector;

[0013] Step 4: Apply boundary conditions and solve for the free surface of seepage;

[0014] Step 5: Define an error threshold. If the difference between the calculated result and the actual result is less than the threshold, the calculated result is considered to be the true free surface position, and proceed to Step 7; if the error check fails, proceed to Step 6.

[0015] Step 6: Using the generated free surface as a reference, locally modify the mesh, update the overall permeation matrix and flow vector, and then execute Step 4 again;

[0016] Step 7: Output parameters within the permeation domain.

[0017] Preferably, step one specifically includes:

[0018] S1.1: Observe the object to be studied and determine its external dimensions;

[0019] S1.2: Determine the upstream and downstream water levels and the dam width, and define them as H1, H2, and H3;

[0020] S1.3: Determine the rock strata distribution and permeability coefficient. If the problem under study is homogeneous, use the same permeability coefficient; if the problem is heterogeneous, use different permeability coefficients.

[0021] S1.4: Based on the above information, establish a seepage model for an earth-rock dam with a free surface.

[0022] Preferably, step two specifically includes:

[0023] S2.1: Perform fixed mesh generation on the object to be analyzed, and store the fixed mesh information and node information;

[0024] S2.2: Based on experience, assume the location of a free seepage surface and set it as the initial free surface.

[0025] Preferably, step three specifically includes:

[0026] S3.1: Calculate the permeability matrix K of the unit e The specific form of the unit penetration matrix is ​​as follows:

[0027]

[0028] S3.2: Matrix element k ij The calculation format is as follows:

[0029] In the formula, D is the permeation tensor. Let i and j represent the basis functions in the i and j directions, respectively, and neither of them is explicitly expressed. For Hamiltonian operators;

[0030] S3.3: Fill the permeability matrix of each unit into the overall permeability matrix K according to the one-to-one correspondence principle;

[0031] S3.4: Calculate the overall flow vector F e Specifically:

[0032] In the formula, For the normalized unit basis function vector:

[0033] q n For traffic; therefore projection Specifically, it can be expressed as follows:

[0034]

[0035] The l here i and l i+1 Let n be the lengths of the two adjacent edges associated with node i. i and n i+1 Let n be the unit outward normal vectors of two adjacent edges. i =(n ix ,n iy ) T A is the polygonal area, x i y i These are the x and y coordinates of the node; n equals the number of sides of the polygon.

[0036] S3.5: Establish basis functions In the polynomial function space P1(Ω) e After projecting ), the basis functions Decomposed into:

[0037]

[0038] because and

[0039] In the formula, P1(Ω) e V1(Ω) represents the polynomial function space. e ) represents the test function space;

[0040] Due to projection residuals Then we have:

[0041]

[0042] Based on the orthogonality condition used in constructing the projection, we know that the second and third terms on the right side of the above equation are zero, thus further simplifying to:

[0043]

[0044] S3.6: Determine the overall permeation matrix form. The first term on the right-hand side of the above equation is accurately calculable; the virtual unit method refers to it as k. ij The consistent term; the second term on the right is called k. ij The stabilizing term is constructed using the values ​​of the projected residuals at the element nodes:

[0045]

[0046] In the above formula, λ is a positive approximation multiplier used to ensure that the order of magnitude of the stable term is comparable to that of the consistency term. It is directly taken as 1. To reasonably consider the influence of the penetration tensor D, given the element k ij The calculation consists of a uniform term and a stable term, and the unit penetration matrix is ​​also divided into two corresponding parts:

[0047] K e =K c +λK s (10)

[0048] S3.7: Derive the multiplier representation, where K c The element is k ij Consistent term, K s The element is k calculated according to equation (9) without using the multiplier λ. ij The stable term; the calculation method for the multiplier λ in equation (10) is as follows: K c λK represents the contribution of the projection of the basis functions to the cell penetration matrix. s This represents the contribution of the projection residuals of the basis functions to the cell penetration matrix; since equation (9) uses the projection residuals at the cell nodes to estimate their contribution, the contribution of the projection can be estimated in the same way, as follows:

[0049]

[0050] S3.8: Determine the multipliers and introduce the reference matrix K refer Its elements are k calculated according to equation (11) without using the multiplier λ. ij Given the reference term, the multiplier λ can be determined as follows: trace represents finding the trace of the matrix.

[0051]

[0052] S3.9: Store the overall penetration matrix K and flow vector information F of the fixed grid. e .

[0053] Preferably, step four specifically includes:

[0054] S4.1: Apply boundary conditions based on the current mesh distribution;

[0055] S4.2: Solve the system of equations: K e H e =F e (13) Determine the position of the free surface;

[0056] S4.3: The location of the seepage point is determined by the quadratic curve intersection method.

[0057] Preferably, step five specifically includes: defining an error threshold ε based on the actual engineering situation and the expected result accuracy, and verifying the position of the free surface to be obtained, wherein the free surface should satisfy hy < ε;

[0058] In the formula, h is the head value at any point on the free surface, and y is the ordinate of that point.

[0059] Preferably, step six specifically includes:

[0060] S6.1: Read the fixed mesh information and cut it using a new free surface;

[0061] S6.2: Divide the elements into three categories: deleted elements (all nodes are above the free surface), native elements (all nodes are below the free surface), and through elements (some nodes are below the free surface, and some nodes are above the free surface);

[0062] S6.3: Define the intersection of the free surface and the penetrating element as the new node. The new node and the nodes of the penetrating element located below the free surface constitute the new element.

[0063] S6.4: Eliminate deleted units and through units, add new units, and retain original units;

[0064] S6.5: The nodes are processed in the same way as the element form to finally form the element node information of the new seepage zone mesh;

[0065] S6.6: Update the overall permeability matrix and flow vector based on the cell node information of the newly generated permeability domain grid.

[0066] Preferably, step seven specifically includes:

[0067] S7.1: Output the location of the free surface of seepage;

[0068] S7.2: Output the head value of each node;

[0069] S7.3: Calculate the location of the seepage point and the iteration error.

[0070] The beneficial effects of this invention are as follows: This invention overcomes the limitations of traditional local mesh modification, achieving automatic densification of nodes on the free surface during the iteration process, resulting in high efficiency and excellent numerical accuracy. This invention is a numerical calculation method for determining the position of the free surface based on the virtual element method, avoiding the impact of mesh distortion on calculation accuracy. Simultaneously, the adoption of local mesh modification during the iteration process reduces storage requirements, significantly improving solution efficiency. It is mainly applied in seepage analysis calculations in geotechnical engineering, hydraulic engineering, and other fields. When designing earth dams, seepage calculations are needed to determine leakage losses and reasonable seepage prevention and drainage measures. Dam types and cross-sectional dimensions often require comparison and selection using seepage calculations. Therefore, accurately calculating the position of the free surface using this method is of great significance to practical engineering. Attached Figure Description

[0071] Figure 1 This is a flowchart of the algorithm of the present invention;

[0072] Figure 2 A schematic diagram comparing the steady-state seepage flow and analytical solution of a rectangular earth-rock dam;

[0073] Figure 3 This is a schematic diagram of local mesh modification;

[0074] Figure 4 This is a schematic diagram of steady-state seepage in a trapezoidal earth-rock dam.

[0075] Figure 5 This is a schematic diagram of the target seepage surface and the calculated seepage surface. Detailed Implementation

[0076] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0077] like Figure 1 As shown, an analytical method for seepage problems with free surfaces based on the virtual element method includes the following steps:

[0078] Step 1: Collect observable data and use it as basic information to establish a seepage model of an earth-rock dam with free surfaces;

[0079] Step 2: Discretize the seepage region to generate a fixed grid;

[0080] Step 3: Calculate the unit permeability matrix and combine them into a global permeability matrix, and simultaneously calculate the global flow vector;

[0081] Step 4: Apply boundary conditions and solve for the free surface of seepage;

[0082] Step 5: Define an error threshold. If the difference between the calculated result and the actual result is less than the threshold, the calculated result is considered to be the true free surface position, and proceed to Step 7; if the error check fails, proceed to Step 6.

[0083] Step 6: Using the generated free surface as a reference, locally modify the mesh, update the overall permeation matrix and flow vector, and then execute Step 4 again;

[0084] Step 7: Output parameters within the permeation domain.

[0085] Preferably, step one specifically includes:

[0086] S1.1: Observe the object to be studied and determine its external dimensions;

[0087] S1.2: Determine the upstream and downstream water levels and the dam width, and define them as H1, H2, and H3;

[0088] S1.3: Determine the rock strata distribution and permeability coefficient. If the problem under study is homogeneous, use the same permeability coefficient; if the problem is heterogeneous, use different permeability coefficients.

[0089] S1.4: Based on the above information, establish a seepage model for an earth-rock dam with a free surface.

[0090] Preferably, step two specifically includes:

[0091] S2.1: Perform fixed mesh generation on the object to be analyzed, and store the fixed mesh information and node information;

[0092] S2.2: Based on experience, assume the location of a free seepage surface and set it as the initial free surface.

[0093] Preferably, step three specifically includes:

[0094] S3.1: Calculate the permeability matrix K of the unit e The specific form of the unit penetration matrix is ​​as follows:

[0095]

[0096] S3.2: Matrix element k ij The calculation format is as follows:

[0097] In the formula, D is the permeation tensor. Let i and j represent the basis functions in the i and j directions, respectively, and neither of them is explicitly expressed. For Hamiltonian operators;

[0098] S3.3: Fill the permeability matrix of each unit into the overall permeability matrix K according to the one-to-one correspondence principle;

[0099] S3.4: Calculate the overall flow vector F e Specifically:

[0100] In the formula, For the normalized unit basis function vector:

[0101] q n For traffic; therefore projection Specifically, it can be expressed as follows:

[0102]

[0103] The l here i and l i+1 Let n be the lengths of the two adjacent edges associated with node i. i and n i+1 Let n be the unit outward normal vectors of two adjacent edges. i =(n ix ,n iy ) T A is the polygonal area, x i y i These are the x and y coordinates of the node; n equals the number of sides of the polygon.

[0104] S3.5: Establish basis functions In the polynomial function space P1(Ω) e After projecting ), the basis functions Decomposed into:

[0105]

[0106] because and

[0107] In the formula, P1(Ω) e V1(Ω) represents the polynomial function space. e ) represents the test function space;

[0108] Due to projection residuals Then we have:

[0109]

[0110] Based on the orthogonality condition used in constructing the projection, we know that the second and third terms on the right side of the above equation are zero, thus further simplifying to:

[0111]

[0112] S3.6: Determine the overall permeation matrix form. The first term on the right-hand side of the above equation is accurately calculable; the virtual unit method refers to it as k. ij The consistent term; the second term on the right is called k. ijThe stabilizing term is constructed using the values ​​of the projected residuals at the element nodes: in this step, the second term on the right is called k. ij The stable term can only be handled in an approximate way. When applying the virtual element method to solve various problems, the appropriate stable term is a hot research topic for scholars at home and abroad.

[0113]

[0114] In the above equation, λ is a positive approximation multiplier used to ensure that the order of magnitude of the stable term is comparable to that of the consistency term. It is directly taken as 1 (in classical literature, it is directly taken as 1). To reasonably consider the influence of the penetration tensor D, given the element k ij The calculation consists of a uniform term and a stable term, and the unit penetration matrix is ​​also divided into two corresponding parts:

[0115] K e =K c +λK s (10)

[0116] S3.7: Derive the multiplier representation, where K c The element is k ij Consistent term, K s The element is k calculated according to equation (9) without using the multiplier λ. ij The stable term; the calculation method for the multiplier λ in equation (10) is as follows: K c λK represents the contribution of the projection of the basis functions to the cell penetration matrix. s This represents the contribution of the projection residuals of the basis functions to the cell penetration matrix; since equation (9) uses the projection residuals at the cell nodes to estimate their contribution, the contribution of the projection can be estimated in the same way, as follows:

[0117]

[0118] S3.8: Determine the multipliers and introduce the reference matrix K refer Its elements are k calculated according to equation (11) without using the multiplier λ. ij Given the reference term, the multiplier λ can be determined as follows: trace represents finding the trace of the matrix.

[0119]

[0120] S3.9: Store the overall penetration matrix K and flow vector information F of the fixed grid. e .

[0121] Preferably, step four specifically includes:

[0122] S4.1: Apply boundary conditions based on the current mesh distribution;

[0123] S4.2: Solve the system of equations: K e H e =F e (13) Determine the position of the free surface;

[0124] S4.3: The location of the seepage point is determined by the quadratic curve intersection method.

[0125] Preferably, step five specifically includes: defining an error threshold ε based on the actual engineering situation and the expected result accuracy, and verifying the position of the free surface to be obtained, wherein the free surface should satisfy hy < ε;

[0126] In the formula, h is the head value at any point on the free surface, and y is the ordinate of that point.

[0127] Preferably, step six specifically includes:

[0128] S6.1: Read the fixed mesh information and cut it using a new free surface;

[0129] S6.2: Divide the elements into three categories: deleted elements (all nodes are above the free surface), native elements (all nodes are below the free surface), and through elements (some nodes are below the free surface, and some nodes are above the free surface);

[0130] S6.3: Define the intersection of the free surface and the penetrating element as the new node. The new node and the nodes of the penetrating element located below the free surface constitute the new element.

[0131] S6.4: Eliminate deleted units and through units, add new units, and retain original units;

[0132] S6.5: The nodes are processed in the same way as the element form to finally form the element node information of the new seepage zone mesh;

[0133] S6.6: Update the overall permeability matrix and flow vector based on the cell node information of the newly generated permeability domain grid.

[0134] Preferably, step seven specifically includes:

[0135] S7.1: Output the location of the free surface of seepage;

[0136] S7.2: Output the head value of each node;

[0137] S7.3: Calculate the location of the seepage point and the iteration error.

[0138] It is worth noting that the virtual element method is applicable to general polygonal meshes, including triangular and quadrilateral elements. When using the first-order virtual element method, there is no need to introduce the calculation of stability terms for triangular elements, and the results are completely consistent with those of three-node triangular elements in conventional finite element methods. For rectangular elements, the stability term calculation method recommended in this paper is used, and the resulting element penetration matrix is ​​highly comparable to that of four-node rectangular bilinear elements in conventional finite element methods.

[0139] Example 1:

[0140] S1.1: Observe the object to be studied and determine its external dimensions;

[0141] S1.2: Determine the upstream and downstream water levels and define them as H1 = 1m and H2 = 2m, and the bottom width of the dam body as H3 = 0.5m;

[0142] S1.3: The rock strata are isotropic, and the research object is a homogeneous earth dam, so the permeability coefficient is k = 1 m / d;

[0143] S1.4: Based on the above information, establish a seepage model for an earth-rock dam with a free surface;

[0144] S2.1: Perform fixed mesh generation on the object to be analyzed, using four-node quadrilateral elements for mesh generation, and store fixed mesh information and node information;

[0145] S2.2: Based on experience, assume the location of a free seepage surface and set it as the initial free surface. In this example, the horizontal line that is level with the upstream water level is taken as the initial free surface.

[0146] Step 3: Calculate the unit penetration matrix according to equation (8) and combine them into the overall penetration matrix. Specifically, this includes the calculation of the consistency term and the calculation of the stability term. At the same time, calculate the overall flow vector according to equation (3).

[0147] Step 4: Apply boundary conditions. For programming convenience, the zero-to-one method is used to apply boundary conditions. That is, the diagonal elements corresponding to the known head values ​​in the total permeability matrix are set to 1, and the rest are set to 0. Substitute the known head values ​​into the corresponding head column vector, and then solve the seepage free surface according to equation (13).

[0148] Step 5: Define the error threshold. In this example, the error threshold ε = 0.05 is taken according to the actual working conditions. If the difference between the calculated result and the actual result is less than the threshold, the calculated result is considered to be the true free surface position, and step 7 is executed; if the error check fails, step 6 is executed.

[0149] Step 6: Using the generated free surface as a reference, locally modify the mesh, update the overall permeation matrix and flow vector, and then execute Step 4 again.

[0150] Taking the 5th iteration as an example, this section will explain in detail the process of mesh modification and the update of the penetration matrix and flow vector.

[0151] S6.1: Read the fixed mesh information and cut it using a new free surface;

[0152] S6.2: As Figure 3 As shown, elements can be divided into three categories: deleted elements, native elements, and through elements. All nodes of a deleted element are above the free surface, and its element type number is set to 0. All nodes of a native element are below the free surface, and its element type number is set to 1. Some nodes of a through element are below the free surface, and some nodes are above the free surface.

[0153] S6.3: Further process the penetrating element. Calculate the intersection points of the free surface and the penetrating element. These intersection points, along with the nodes of the penetrating element located below the free surface, constitute new elements. Set the element type number of the new elements to 2, and simultaneously change the type of the penetrating element to deleted element, setting its element type number to 0.

[0154] S6.4: with Figure 3 Taking the processing within the dashed box as an example, the original unit 6-16-15-8 requires no modification, the three new units are B-6-8-C, C-8-C′ and C′-8-15-14-D, and the corresponding deleted units are 5-6-8-7, 7-8-10-9 and 8-15-14-10.

[0155] S6.5: Following the cell type classification method, cell nodes are also divided into three categories: deleted nodes, native nodes, and newly created nodes, with their type numbers set to 0, 1, and 2 respectively. Figure 3 For example, nodes 3, 5, 7, 9, 10, 11, and 12 are deleted nodes, nodes 2, 4, 6, 8, 13, 14, 15, and 16 are native nodes, and nodes A, B, C, C′, D, and E are newly created nodes. Based on the original fixed mesh information, the total number of elements and nodes increases by 1 when one new element and node appear. The node composition of the new element and the coordinates of the new node are recorded, ultimately forming the element node information of the new seepage region mesh. Although the new seepage region is actually composed only of native elements and newly created elements, its mesh contains all deleted elements, native elements, and through elements, and correspondingly, its nodes also contain all deleted nodes, native nodes, and newly created nodes.

[0156] Based on the cell node information of the newly generated seepage domain mesh, update the overall seepage matrix and flow vector as follows:

[0157] First, the overall penetration matrix K and the overall flow vector F obtained based on the fixed grid are stored;

[0158] Then, after the new seepage zone is cut and formed, considering that the number of nodes in the new grid is greater than that in the original fixed grid, K and F are expanded to K′ and F′ respectively based on the total number of nodes in the new grid.

[0159] K′ and F′ are updated according to the cell type. If the contribution of the original cell remains unchanged, there is no need to modify K′ and F′. ​​If the cell is deleted and no longer belongs to the new mesh, its cell permeability matrix and cell flow vector are calculated and removed from K′ and F′. ​​For newly created cells, its cell permeability matrix and cell flow vector are calculated and added to K′ and F′. ​​After all cells have been processed, the calculation of the overall permeability matrix and overall flow vector of the new mesh is completed.

[0160] Step 7: Output parameters within the permeation domain;

[0161] After seven iterations, the results stabilized. The iteration process and the final free surface position information are shown in Table 1. For ease of presentation, both the target information and the calculation results are plotted in Table 1. Figure 4 .

[0162] Table 1 Location of the free surface of seepage in rectangular earth-rock dams

[0163]

[0164] Using the method proposed in this invention, the calculation results obtained after 7 iterations are in high agreement with the analytical solution for the free surface. The final location of the seepage point is 0.6594m.

[0165] Example 2:

[0166] S1.1: Observe the object to be studied and determine its external dimensions;

[0167] S1.2: Determine the upstream and downstream water levels and define them as H1 = 5m and H2 = 1m, and the bottom width of the dam body as H3 = 7m;

[0168] S1.3: The rock strata are isotropic, and the research object is a homogeneous earth dam, so the permeability coefficient is k = 1 m / d;

[0169] S1.4: Based on the above information, establish a seepage model for an earth-rock dam with a free surface;

[0170] S2.1: A fixed mesh is generated for the object to be analyzed. Without loss of generality, a combination of three-node triangles and four-node quadrilaterals is used for mesh generation in this example. After processing the entire domain, the nodes and elements are numbered sequentially, and the fixed mesh information and node information are stored.

[0171] S2.2: Based on experience, assume the location of a free seepage surface and set it as the initial free surface. In this example, the horizontal line that is level with the upstream water level is taken as the initial free surface.

[0172] Step 3: Calculate the unit penetration matrix according to equation (8) and combine them into the overall penetration matrix. Specifically, this includes the calculation of the consistency term and the calculation of the stability term. At the same time, calculate the overall flow vector according to equation (3).

[0173] Step 4: Apply boundary conditions. For programming convenience, the zero-to-one method is used to apply boundary conditions. That is, the diagonal elements corresponding to the known head values ​​in the total permeability matrix are set to 1, and the rest are set to 0. Substitute the known head values ​​into the corresponding head column vector, and then solve the seepage free surface according to equation (13).

[0174] Step 5: Define the error threshold. In this example, the error threshold ε = 0.05 is taken according to the actual working conditions. If the difference between the calculated result and the actual result is less than the threshold, the calculated result is considered to be the true free surface position, and step 7 is executed; if the error check fails, step 6 is executed.

[0175] Step 6: Using the generated free surface as a reference, locally modify the mesh, update the overall permeation matrix and flow vector, and then execute Step 4 again.

[0176] Modifying the mesh follows the same method as the previous example, and will not be explained in detail here;

[0177] Step 7: Output parameters within the permeation domain;

[0178] After 17 iterations, the results tended to stabilize. The iteration process and the final free surface position information are shown in Table 2. For ease of presentation, both the target information and the calculation results are plotted on [Table 2]. Figure 5 .

[0179] Table 2 Location of the free surface of seepage in trapezoidal earth-rock dams

[0180]

[0181] Using the method proposed in this invention, the calculation results obtained after 17 iterations are consistent with the actual free surface height, and the final location of the seepage point is 3.579m. The following conclusions can be drawn:

[0182] Using the virtual element method (VEM) to solve seepage problems with free surfaces can effectively eliminate the impact of mesh distortion on computational accuracy.

[0183] During local mesh modification, the cut-based seepage region modification method can reduce the storage space requirements and improve computational efficiency to some extent.

[0184] The above embodiments are merely preferred technical solutions of the present invention and should not be considered as limitations on the present invention. The embodiments and features described in these embodiments can be arbitrarily combined without conflict. The scope of protection of the present invention should be limited to the technical solutions described in the claims, including equivalent substitutions of the technical features described in the claims. That is, equivalent substitutions and improvements within this scope are also within the scope of protection of the present invention.

Claims

1. A method for analyzing free-surface seepage problems based on the virtual element method, characterized in that: It includes the following steps: Step 1: Collect observable data and use it as the basic information to build a seepage model of an earth-rock dam with free surfaces; Step 2: Discretize the seepage region to generate a fixed grid; Step 3: Calculate the unit permeability matrix and combine them into a global permeability matrix, and simultaneously calculate the global flow vector; Step 4: Apply boundary conditions and solve for the free surface of seepage; Step 5: Define an error threshold. If the difference between the calculated result and the actual result is less than the threshold, the calculated result is considered to be the true free surface position, and proceed to Step 7. If the error check fails, proceed to step six; Step 6: Using the generated free surface as a reference, locally modify the mesh, update the overall permeation matrix and flow vector, and then execute Step 4 again; Step 7: Output parameters within the penetration domain; Step three specifically includes: S3.1: Calculate the permeability matrix of the unit K e The specific form of the unit penetration matrix is ​​as follows: （1）； S3.2: Matrix elements in the formula k ij The calculation format is as follows: (2); In the formula, D For the permeation tensor; φ i , φ j Let i and j represent the basis functions in the i and j directions, respectively, and neither of them is explicitly expressed. For Hamiltonian operators; S3.3: Fill the overall permeability matrix with the permeability matrix of each unit according to the one-to-one correspondence principle. K Inside; S3.4: Calculate the overall flow vector F e Specifically: (3); In the formula, φ e For the normalized unit basis function vector: (4); q n For traffic; therefore φ i projection Specifically, it can be expressed as follows: （5）； Here l i and l i+1 Let i be the length of the two adjacent edges associated with node i. n i and n i+1 Let be the unit outward normal vectors of the two adjacent edges. n i = ( n ix , n iy ) T , A For the polygonal area, x i y i Here are the x and y coordinates of the node, and n equals the number of sides of the polygon; S3.5: Establish basis functions φ i In the polynomial function space P 1 (Ω) e After projecting the basis functions, φ i Decomposed into: (6); because φ i ∈ V 1 (Ω) e ), Π φ i ∈ P 1 (Ω) e ),and P 1 (Ω) e ) V 1 (Ω) e ), In the formula, P 1 (Ω) e ) represents the polynomial function space, V 1 (Ω) e () represents the test function space; Due to the projection residual ( φ i -Π φ i )∈ V 1 (Ω) e Then we have: （7）； Based on the orthogonality condition used in constructing the projection, we know that the second and third terms on the right side of the above equation are zero, thus further simplifying to: （8）； S3.6: Determine the overall permeability matrix form. The first term on the right-hand side of the above equation is accurately calculable; the virtual unit method refers to it as... k ij The first term is the consistent term; the second term on the right is called the consistent term. k ij The stabilizing term is constructed using the values ​​of the projected residuals at the element nodes: （9）； In the above formula x m Let m be the x-coordinate of point m. λ To ensure that the order of magnitude of the stable term is comparable to that of the uniform term, the positive approximation multiplier is directly set to 1. This is to reasonably account for the penetration tensor. D The influence of elements k ij The calculation consists of a uniform term and a stable term, and the unit penetration matrix is ​​also divided into two corresponding parts: （10）； S3.7: Deriving the multiplier representation, here K c The elements are k ij Consistent terms, K s The element is one that does not use multipliers. λ Calculated according to formula (9) under the following circumstances k ij Stable term; multipliers in equation (10) λ The calculation approach is as follows: K c This represents the contribution of the projection of the basis functions to the cell penetration matrix. λK s This represents the contribution of the projection residuals of the basis functions to the cell penetration matrix; Since the contribution of the projection residual is estimated at the element node value in equation (9), the contribution of the projection can be estimated in the same way, as follows: （11）； S3.8: Determine the multipliers and introduce the reference matrix. K refer Its elements are those that do not use multipliers. λ Calculated according to formula (11) under the following circumstances k ij The reference term, then the multiplier λ It can be determined as follows: trace represents finding the trace of a matrix. （12）； S3.9: Store the overall penetration matrix of the fixed grid. K and flow vector information F e .

2. The analytical method for free surface seepage problems based on the virtual element method according to claim 1, characterized in that: Step one specifically includes: S1.1: Observe the object to be studied and determine its external dimensions; S1.2: Determine the upstream and downstream water levels and define them as H1 and H2, and define the dam width as H3; S1.3: Determine the rock strata distribution and permeability coefficient. If the problem under study is homogeneous, use the same permeability coefficient; if the problem is heterogeneous, use different permeability coefficients. S1.4: Based on the above information, establish a seepage model for an earth-rock dam with a free surface.

3. The analytical method for free surface seepage problems based on the virtual element method according to claim 1, characterized in that: Step two specifically includes: S2.1: Perform fixed mesh generation on the object to be analyzed, and store the fixed mesh information and node information; S2.2: Based on experience, assume the location of a free seepage surface and set it as the initial free surface.

4. The analytical method for free surface seepage problems based on the virtual element method according to claim 1, characterized in that: Step four specifically includes: S4.1: Apply boundary conditions based on the current mesh distribution; S4.2: Solve the system of equations: Determine the position of the free surface; S4.3: The location of the seepage point is determined by the quadratic curve intersection method.

5. The analytical method for free surface seepage problems based on the virtual element method according to claim 1, characterized in that: Step five specifically includes: defining an error threshold based on the actual engineering situation and the expected accuracy of the results. The position of the free surface is checked, and the free surface should satisfy the following conditions: hy < ; In the formula h Let the head value be at any point on the free surface. y Let be the ordinate of that point.

6. The analytical method for free surface seepage problems based on the virtual element method according to claim 1, characterized in that: Step six specifically includes: S6.1: Read the fixed mesh information and cut it using a new free surface; S6.2: The elements are divided into three categories: deleted elements, native elements, and through elements; where deleted elements are all nodes above the free surface, native elements are all nodes below the free surface, and through elements are some nodes below the free surface and some nodes above the free surface. S6.3: Define the intersection of the free surface and the penetrating element as the new node. The new node and the nodes of the penetrating element located below the free surface constitute the new element. S6.4: Eliminate deleted units and through units, add new units, and retain original units; S6.5: The nodes are processed in the same way as the element form to finally form the element node information of the new seepage zone mesh; S6.6: Update the overall permeability matrix and flow vector based on the cell node information of the newly generated permeability domain grid.

7. The analytical method for free surface seepage problems based on the virtual element method according to claim 1, characterized in that: Step seven specifically includes: S7.1: Output the location of the free surface of seepage; S7.2: Output the head value of each node; S7.3: Calculate the location of the seepage point and the iteration error.

Images