Efficient numerical solution method for steady-state water head field of dam containing complex foundation

Abstract

The invention discloses an efficient numerical solution method for a steady-state water head field of a dam containing a complex foundation, and relates to the technical field of geotechnical engineering. The method comprises the following steps: establishing a seepage mathematical model in a Cartesian coordinate system by acquiring foundation parameters and boundary conditions, introducing similar lines parallel to a far-field boundary to construct an improved proportional boundary finite element coordinate system mapping relation, and discretizing and deducing a variable coefficient control equation by adopting a one-dimensional high-order spectrum unit; an equation is simplified into a Riccati equation by analyzing the characteristics of a coefficient matrix, an infinite seepage matrix is solved in combination with Schur decomposition, a seepage matrix at a boundary is obtained through iteration of a fourth-order Runge-Kutta method with the infinite seepage matrix as an initial value, and finally steady-state water head field distribution is obtained. According to the method, the continuity of the interface water head and the flowing boundary condition is guaranteed, the precision and efficiency of complex foundation seepage solving are greatly improved, and reliable technical support is provided for engineering seepage-proofing design and stability evaluation.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering technology, specifically to an efficient numerical solution method for the steady-state head field of dams with complex foundations. Background Technology

[0002] Steady-state seepage analysis of complex foundations is a core technical challenge in the fields of hydrogeology and civil engineering, and is widely applicable to practical engineering scenarios such as dam foundation seepage prevention design, levee seepage stability assessment, groundwater resource development, and pollutant migration prediction. In these engineering scenarios, the foundation often exhibits multi-layered, heterogeneous, and anisotropic characteristics, and its spatial range extends to a semi-infinite domain. Accurate prediction of the head field distribution and seepage flow directly affects the stability of the engineering structure, the rationality of the seepage prevention scheme, and the effectiveness of environmental safety protection, and is a key prerequisite for ensuring the long-term reliable operation of the project.

[0003] In the aforementioned engineering applications, traditional numerical methods face significant technical bottlenecks: the finite element method requires discretizing the entire computational domain into a mesh, and to ensure the accuracy of the far-field boundary, the model needs to be over-extended, leading to a surge in computational units and excessive resource consumption, making it difficult to meet the high-efficiency solution requirements of complex foundations; while the boundary element method has high solution efficiency, it has poor adaptability to the heterogeneity of the medium and complex boundary conditions, and cannot accurately capture the seepage laws inside the foundation; although the traditional scaled boundary finite element method (SBFEM) is naturally adapted to infinite domain modeling, a single scale center is difficult to describe the geometric non-uniformity of complex foundations, and the multi-scale center scheme is prone to causing discontinuities in interface head and flux, affecting computational accuracy. In addition, artificially truncated boundaries in half-space problems can easily introduce errors, further reducing the consistency between numerical solutions and engineering realities.

[0004] Engineering practice has placed urgent demands on the accuracy, efficiency, and adaptability of solutions for steady-state head fields in complex foundations: methods must accurately characterize core features such as heterogeneity, anisotropy, and semi-infinite domains, and possess efficient computational capabilities to support engineering design and decision-making, focusing on practical technical challenges in specific engineering scenarios.

[0005] Therefore, developing a numerical solution method that balances accuracy, efficiency, and engineering adaptability to address complex foundation seepage problems in dam foundations, embankments, and other engineering projects is crucial for advancing related engineering technologies and ensuring project safety. Summary of the Invention

[0006] Based on the above-mentioned technical problems, this application discloses an efficient numerical solution method for the steady-state head field of dams with complex foundations, specifically including:

[0007] The medium characteristic parameters, boundary conditions, and geometric parameters of complex foundations are obtained. The medium characteristic parameters include permeability coefficient, porosity, and compressibility coefficient. The boundary conditions include head boundary, flow boundary, and infinity boundary.

[0008] A mathematical model of steady-state seepage in a half-space is established in a Cartesian coordinate system. The governing equations are constructed based on Darcy's law and Laplace's equation. The Galerkin weighted residual method is used to weaken the function continuity requirement. The characteristics of the half-infinite domain are handled by Fourier transform.

[0009] By introducing similar lines parallel to the far-field boundary, a mapping relationship between the Cartesian coordinate system and the improved scaled boundary finite element coordinate system is established, and the differential operator transformation is realized through the Jacobian matrix.

[0010] In the improved proportional boundary finite element coordinate system, one-dimensional high-order spectral elements are used to discretize the boundary lines and similarity lines, and the variable coefficient steady-state seepage control equation is derived based on the weighted residual method and integral by parts.

[0011] By analyzing the asymptotic properties of the coefficient matrix at infinity, the governing equations are simplified to the Riccati equations. The Hamiltonian matrix is ​​introduced and the steady-state seepage matrix at infinity is solved by Schur decomposition.

[0012] Using the steady-state seepage matrix at infinity as the initial value, the fourth-order Runge-Kutta method is used to iterate backward along the radial coordinate to obtain the steady-state seepage matrix at the boundary, and then the steady-state head field distribution of the dam is obtained.

[0013] Preferably, the governing equation is: , ,in, It is the source per unit volume. It is the main pressure head. This is the derivative of the total pressure head with respect to time. It is a specific storage coefficient. It is the density of the fluid. It is gravitational acceleration. It is the porosity of the medium. and These are the compressibility coefficients of the solid and the fluid, respectively. It is the exhaust velocity vector. It is a penetration matrix. and respectively along the medium and Permeability coefficient in directional direction, in isotropic media In anisotropic media .

[0014] Preferably, the governing equations are transformed to the frequency domain using Fourier transform, resulting in: , ,exist Above, among which, , and They are , and Fourier transform;

[0015] The boundary conditions are specifically as follows: ,exist superior; ,exist Above, among which, For the entire boundary, The water head boundary, For the velocity boundary, To specify the head value, Specify the flow rate value.

[0016] Preferably, the equation for the mapping relationship is: , ,in, and The value of a point in the problem domain. and To scale the coordinates of the splicing line, and The coordinates of the boundary line, To improve the radial coordinates of the scaled boundary finite element coordinate system;

[0017] The coordinates of nodes in boundaries and similar lines are given by the shape function formula: , ,in, The shape function of the nodal line element. and In line unit Cartesian coordinates of the nodes.

[0018] Preferably, the derivative transformation between the Cartesian coordinate system and the improved scaled boundary finite element coordinate system is based on the Jacobian matrix, the determinant of which is: The derivative transformation relationship is: ,in, To improve the circumferential coordinates of the proportional boundary finite element coordinate system.

[0019] Preferably, the field variables at any point within the analysis domain are obtained through interpolation: ,in, The water head of the points on the ray connecting the nodes of similar line elements and the nodes of the boundary line elements is the radial coordinate. The function;

[0020] The variable-coefficient steady-state seepage control equation is: ,in, , , These depend on radial coordinates. The coefficient matrix, right The partial derivatives, The transpose of the matrix, right The partial derivatives, right The second-order partial derivative, right The first-order partial derivative, The nodal flux magnitude generated for the normal flux magnitude on the line connecting similar lines and artificial boundaries. The amplitude is the internal source amplitude.

[0021] Preferably, the elements of the coefficient matrix are expressed in rational fraction form as follows: , , ,in , , They are respectively , and The Middle Line number Column elements, For the row and column indices of the matrix elements, It is about polynomial, , , These are the polynomials corresponding to the coefficient matrices, with the first subscript indicating the polynomial's... The highest power, the second subscript corresponds to the coefficient matrix. It is the determinant of the Jacobian matrix, and is about... A first-order polynomial.

[0022] Preferably, the Riccati equation is: ,in, The steady-state seepage matrix at infinity; for The inverse matrix;

[0023] By introducing intermediate variables Transform the second-order differential equation into a first-order ordinary differential equation: ,in, It is a Hamiltonian matrix. This is the head-related vector. This is the node flux-related vector. right The first-order partial derivative, For Hamiltonian matrices; Schur decomposition yields ,in The eigenvector matrix, for The inverse matrix, For an eigenvalue matrix, select eigenvector submatrices corresponding to eigenvalues ​​with negative real parts. , ,get , for The inverse matrix.

[0024] Preferably, the iterative steps of the fourth-order Runge-Kutta method are as follows:

[0025] ;

[0026] ;

[0027] ;

[0028] ;

[0029] ;in, It is the iteration step size. , , , They are respectively , , , Gradient of position, For the first The radial coordinate values ​​of the next iteration For the first The radial coordinate values ​​of the next iteration For the first The seepage matrix obtained in the second iteration For the first The seepage matrix obtained in the second iteration is iterated to... hour, ] represents the steady-state seepage matrix at the boundary.

[0030] Preferably, the one-dimensional higher-order spectral unit adopts Legendre polynomial shape functions, which are suitable for semi-infinite domain steady-state seepage analysis scenarios with complex foundations containing multiple layers, heterogeneity, and anisotropy, such as dam foundations, embankments, and underground engineering. The parallel setting of similar lines and far-field boundaries ensures the continuity of interface head and flow boundary conditions. Combined with the solution process, it realizes accurate prediction and efficient calculation of steady-state head field in engineering practice.

[0031] Compared with the prior art, the technical solution of this application has the following technical effects:

[0032] This invention achieves a significant breakthrough in the accuracy of solving seepage problems in complex foundations. By innovatively introducing the concept of similar lines parallel to the far-field boundary, combined with the accurate mapping of the improved proportional boundary finite element coordinate system, it can effectively capture the heterogeneous, anisotropic, and semi-infinite domain characteristics of the foundation, avoiding the errors introduced by the manual truncation of boundaries in traditional methods. Simultaneously, the governing equations constructed based on Darcy's law and Laplace's equation, combined with the spatial integration weakening treatment of the Galerkin weighted residual method, can accurately describe the continuity of head and flux at different medium interfaces, greatly improving the prediction accuracy of steady-state head field distribution and providing reliable data support for engineering seepage prevention design.

[0033] This invention only requires discretization of the boundary lines and similar lines using one-dimensional high-order spectral units, combined with the efficient iteration of the fourth-order Runge-Kutta method, significantly reducing the number of computational units and degrees of freedom. It can handle semi-infinite domain problems without extending the far-field model, avoiding redundant computation, significantly reducing hardware resource consumption, and enabling rapid response to efficient solution requirements in engineering design such as multi-scheme comparison and parameter optimization.

[0034] This invention can be widely applied to seepage analysis scenarios involving complex foundations, such as dam foundations, embankments, and underground engineering. Whether it is pressurized seepage in multi-level heterogeneous foundations or unpressurized seepage in a combined dam and foundation system, it can achieve accurate solutions by flexibly adjusting the coordinate system mapping relationship and boundary condition processing method. It can adapt to different engineering needs without significantly adjusting the core algorithm, reducing the secondary development cost of method adaptation, and providing a unified and efficient solution path for various complex foundation seepage problems.

[0035] This invention effectively fills the gap in the balance between accuracy, efficiency, and adaptability of traditional methods. Through accurate prediction of the hydraulic head field, it can help engineers optimize the layout of seepage prevention structures, rationally select seepage prevention material parameters, and reduce the risk of seepage damage. Its efficient computational characteristics support rapid iterative verification in the engineering design stage, shortening the design cycle. Its broad scenario adaptability reduces the cost of switching solution methods in different engineering scenarios, providing stable and reliable technical support for seepage analysis of complex foundations, and promoting technological progress and safety assurance in related engineering fields.

[0036] The above description is only an overview of the technical solution of this application. In order to better understand the technical means of this application and implement it in accordance with the contents of the specification, and to make the above and other objects, features and advantages of this application more obvious and understandable, the preferred embodiments of this application are described in detail below with reference to the accompanying drawings.

[0037] The above and other objects, advantages and features of this application will become more apparent to those skilled in the art from the following detailed description of specific embodiments in conjunction with the accompanying drawings. Attached Figure Description

[0038] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In all drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn to scale.

[0039] Based on the description of the figures and their corresponding technical content in the document, the titles of the figures are as follows:

[0040] Figure 1 This is a flowchart illustrating an efficient numerical solution method for the steady-state head field of a dam with complex foundations.

[0041] Figure 2 A schematic diagram of coordinate transformation for an efficient numerical solution method for the steady-state head field of a dam with complex foundations;

[0042] Figure 3 This is a schematic diagram of the rectangular river channel calculation model and unit discretization in an embodiment of the present invention;

[0043] Figure 4 This is a comparison diagram of the seepage field calculated by this invention and other methods in an embodiment of the invention;

[0044] Figure 5 This is a schematic diagram of the calculation model and unit discretization of the dam-foundation-reservoir system according to an embodiment of the present invention;

[0045] Figure 6 This is a comparison chart of the head calculated by this invention and other methods. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. In the following description, specific details such as specific configurations and components are provided merely to help fully understand the embodiments of this application. Therefore, those skilled in the art should understand that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of this application. In addition, for clarity and brevity, descriptions of known functions and structures are omitted in the embodiments.

[0047] It should be understood that the phrase "an embodiment" or "this embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "an embodiment" or "this embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments.

[0048] Furthermore, reference numerals and / or letters may be repeated in different examples within this application. Such repetition is for the purpose of simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or settings discussed.

[0049] In this article, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can mean: A exists alone, B exists alone, and A and B exist simultaneously. The term " / and" in this article describes another type of relationship between related objects, indicating that two relationships can exist. For example, A / and B can mean: A exists alone, and A and B exist alone. In addition, the character " / " in this article generally indicates that the related objects before and after it are in an "or" relationship.

[0050] In this article, the term "at least one" is merely a description of the relationship between related objects, indicating that there can be three relationships. For example, "at least one of A and B" can mean: A exists alone, A and B exist simultaneously, or B exists alone.

[0051] It should also be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion.

[0052] Example 1

[0053] This embodiment details an efficient numerical solution method for the steady-state head field of a dam with complex foundations, such as... Figure 1 As shown, specifically:

[0054] The medium characteristic parameters, boundary conditions, and geometric parameters of complex foundations are obtained. The medium characteristic parameters include permeability coefficient, porosity, and compressibility coefficient. The boundary conditions include head boundary, flow boundary, and infinity boundary.

[0055] A mathematical model of steady-state seepage in a half-space is established in a Cartesian coordinate system. The governing equations are constructed based on Darcy's law and Laplace's equation. The Galerkin weighted residual method is used to weaken the function continuity requirement. The characteristics of the half-infinite domain are handled by Fourier transform.

[0056] A semi-spatial steady-state seepage mathematical model is established in Cartesian coordinates. The governing equations are constructed based on Darcy's law and Laplace's equation, taking into account the anisotropy of the medium. The Galerkin weighted residual method is used for spatial integration weakening to reduce the continuity requirements of the functions. Three types of boundary conditions (head boundary, flow boundary, and far-field boundary) are strictly distinguished, and the semi-infinite domain characteristics are handled through Fourier transform.

[0057] By introducing similar lines parallel to the far-field boundary, a mapping relationship between the Cartesian coordinate system and the improved scaled boundary finite element coordinate system is established, and the differential operator transformation is realized through the Jacobian matrix.

[0058] The concept of similar lines parallel to the far-field boundary is introduced, and an exact mapping is established from the Cartesian coordinate system (x,y) to the improved SBFEM coordinate system (ξ,η): the radial coordinate ξ controls the scaling from the similar line (ξ=0) to the artificial boundary (ξ=1) up to infinity, and the circumferential coordinate η describes the geometric changes of the boundary. Geometric similarity is preserved, and the exact transformation of the differential operator is achieved through the Jacobian matrix.

[0059] In the improved proportional boundary finite element coordinate system, one-dimensional high-order spectral elements are used to discretize the boundary lines and similarity lines, and the variable coefficient steady-state seepage control equation is derived based on the weighted residual method and integral by parts.

[0060] In the improved SBFEM coordinate system, one-dimensional high-order spectral elements are used to discretize the boundary lines and similarity lines, and Legendre polynomial shape functions are used to ensure exponential convergence of the numerical solution. Through the weighted residual method combined with integration by parts, a steady-state seepage control equation with variable coefficient characteristics is obtained, where the coefficient matrix is ​​a function of the radial coordinate ξ.

[0061] like Figure 2 As shown, the top sub-image first anchors the core elements of the physical space, in blue. The connecting line represents the artificial boundary (corresponding to the seepage prevention structure interface of the dam foundation and dike in the project), in red. The connecting lines are similar lines parallel to the far-field interface (corresponding to the far-field boundary of the foundation extending to the semi-infinite domain), and the black dashed lines divide the solution regions from domain 1 to domain 5 (matching the multi-level partitioning of complex foundations). At this point, the parallel relationship between the similar lines and the boundary lines has already provided geometric constraints for the subsequent continuity of interface parameters;

[0062] The intermediate subgraph presents the physical space details in the Cartesian coordinate system: the positions of boundary lines and similarity lines are clearly marked, as well as the distribution of domains 1 to 8 and the nodes. and The arrangement of nodes (these nodes correspond to the interfaces of different media in the foundation and the contact points of the seepage prevention structure); the transformation operation starts here: with similar lines as reference benchmarks, the coordinates of any point in the physical space are associated by radial scaling from similar lines to boundary lines. That is, the position of each point is determined by the corresponding point on the similar line, the corresponding point on the boundary line, and the radial parameters. At the same time, the spatial derivative is transformed through the geometric transformation matrix, transforming the irregular multi-domain morphology in the physical space into a structured form that depends on radial and circumferential parameters.

[0063] The bottom sub-figure is the transformed scaled boundary finite element coordinate system. The originally dispersed subdomains and nodes form a continuous solution domain under the new coordinate system. The radial parameters progress from the similarity lines to the boundary lines, corresponding to the spatial extension of the foundation from the far field to the engineering interface. The circumferential parameters match the circumferential geometric characteristics of the foundation. This not only preserves the heterogeneous and semi-infinite domain properties of complex foundations, but also simplifies the multidimensional discretization of physical space to the one-dimensional discretization of the boundary lines and similarity lines. Subsequent radial iterative calculations can efficiently obtain the steady-state head field, perfectly meeting the seepage analysis needs of dam foundations, levees and other engineering projects.

[0064] By analyzing the asymptotic properties of the coefficient matrix at infinity, the governing equations are simplified to the Riccati equations. The Hamiltonian matrix is ​​introduced and the steady-state seepage matrix at infinity is solved by Schur decomposition.

[0065] Analyzing the asymptotic properties of the coefficient matrix at infinity, the governing equations are simplified to a function relating to the steady-state seepage matrix [K]. ∞ The Riccati equation is used. A Hamiltonian matrix is ​​introduced, and the Schur decomposition technique is employed to solve the corresponding eigenvalue problem. By selecting eigenvector submatrices corresponding to eigenvalues ​​with negative real parts, the steady-state seepage matrix at infinity is obtained, providing initial values ​​for subsequent iterations.

[0066] Using the steady-state seepage matrix at infinity as the initial value, the fourth-order Runge-Kutta method is used to iterate backward along the radial coordinate to obtain the steady-state seepage matrix at the boundary, and then the steady-state head field distribution of the dam is obtained.

[0067] Furthermore, in the Cartesian coordinate system The following structural model is established for the complex half-space steady-state seepage problem. Based on the Laplace equation and Darcy's law, the governing equation for the general seepage problem in soil is: , , ,in It is the source per unit volume. It is the main pressure head. This represents the derivative of the total pressure head with respect to time 167. It is a specific storage coefficient, where, It is the density of the fluid. It is gravitational acceleration. It is the porosity of the medium. and These are the compressibility coefficients of the solid and the fluid, respectively. It is the exhaust velocity vector. It is the permeability matrix, where and respectively along the medium and The permeability coefficient in the direction of the medium, in isotropic media. In anisotropic media, ;

[0068] The percolation matrix at infinity is accurately derived using Fourier transform, and then transformed to the frequency domain: , exist Above, where the symbol , and They are , and Fourier transform;

[0069] For confined seepage problems, the total head, velocity, or flow rate must be specified at the boundary of the domain. The entire boundary is designated as Γ, and the head boundary is designated as Γ. h The velocity boundary is specified as Γ v Boundary conditions can be specified as follows: , , , ;

[0070] Furthermore, a proportional boundary finite element coordinate system is established. The equations for the transformation between the global Cartesian coordinate system and the scaled boundary local coordinate system are expressed as:

[0071]

[0072]

[0073] in, and Represents the value of a point in the problem domain; and Indicates the coordinates of the scaling line; and Represents the coordinates of the boundary line;

[0074] Furthermore, the boundaries and corresponding scaling lines are discretized using the same elements, and points in the boundaries and scaling lines can be represented by the same shape function: , ,in, The shape function representing an n-node line element. and Let be the Cartesian coordinates of n nodes in the line element;

[0075] After substituting the shape function into the transformation relation: , ;

[0076] The derivative transformation between the Cartesian coordinate system and the SBFEM coordinate system is based on the Jacobian matrix: ;

[0077] Its transformation relationship can be expressed as:

[0078]

[0079]

[0080]

[0081] .

[0082] Furthermore, analyze any point within the domain. The field variables can be obtained through interpolation: ,in, The shape function matrix representing the head interpolation of elements in boundary lines and similar lines; The head of the point on the ray that connects to the nodes of similar line elements and boundary line elements is represented by radial coordinates. The function.

[0083] Furthermore, substituting the Laplace operator into the percolation matrix:

[0084]

[0085]

[0086] in, , .

[0087] Applying the weighted residual method to integrate over the analysis domain, we obtain the following equation:

[0088]

[0089]

[0090] in, Represent the analysis domain (bounded or unbounded), and The vector representing the weighting function. It is an infinitesimally small area.

[0091] The weighting function can be formulated using the same shape function as the head, and therefore it can be expressed as:

[0092]

[0093] The first two items combined are represented as Substituting the weighting function, we get the following equation:

[0094]

[0095] Furthermore, applying the product rule of differentials, we derive:

[0096]

[0097] Substitution and The relationship can be simplified to:

[0098]

[0099] For an infinite domain, the local coordinates are ξ ≥ 1 and −1 ≤ η ≤ 1. Separate integration is performed with respect to the directions ξ and η:

[0100]

[0101] Similarly, the remaining part can be represented as:

[0102]

[0103] Introducing the coefficient matrix:

[0104]

[0105] The governing equations in the form of weighted residuals can be expressed as:

[0106]

[0107] Furthermore, assuming and The sum is always zero, taking into account the weighting function. The arbitrariness, for any The following conditions must be met:

[0108]

[0109]

[0110] in, The nodal flux magnitude corresponds to the normal flux magnitude generated by the lines connecting the scaling splice line and the artificial boundary; This corresponds to the internal source amplitude. This condition represents the seepage control equation of the improved SBFEM, which is a second-order linear nonhomogeneous differential equation. It is worth noting that in this equation, , , It is no longer the constant coefficient matrix in the original SBFEM, but it is related to the radial coordinates. The relevant function. This modification is due to the existence of functions related to the Jacobian matrix. The associated additional constant matrix of terms causes the Jacobian matrix to not be represented as ;

[0111] Furthermore, the principle of virtual work is applied to the flux of internal nodes. ,get:

[0112]

[0113]

[0114]

[0115]

[0116] in Defined as a curve radial coordinates The nodal flux components on the surface Let n be the length of any infinitesimal arc with n=const in the domain.

[0117] After combining, we get:

[0118]

[0119] Considering internal node flux Given the arbitrariness of the internal node flux, the expression for the internal node flux is:

[0120]

[0121] Therefore, in an infinite domain, the flux of the outer equivalent nodes is opposite to that of the inner nodes:

[0122]

[0123] Furthermore, within the analysis area, the water head and external node flow The following relationship exists:

[0124]

[0125]

[0126] Considering the case where both the normal flux magnitude and the internal source are zero, it leads to and ,as well as The arbitrariness of the coefficient term, which requires it to be equal to zero, leads to the derivation of the above equation concerning the radial coordinates:

[0127]

[0128] in The expression is as follows:

[0129]

[0130] The above equation is the governing equation of the steady-state seepage matrix of the modified SBFEM, which is a first-order nonlinear variable-coefficient partial differential equation.

[0131] Furthermore, based on the determinant of the Jacobian matrix, the inverse of the Jacobian matrix, and the matrix... Matrix and matrix Matrix, coefficient matrix , and The elements of can be represented in rational fraction form as follows:

[0132]

[0133]

[0134]

[0135] in yes The first subscript of the polynomial indicates that n is in the polynomial. The highest power, the second subscript corresponds to the coefficient matrix. , and The denominator is the determinant of the Jacobian matrix, and the subscript "1" indicates that the determinant is about... A first-order polynomial whose constant term is not zero. At that time, the following relationship exists.

[0136]

[0137]

[0138]

[0139] when When the steady-state seepage matrix is ​​denoted as The steady-state seepage control equation can be written as:

[0140]

[0141] Because of the above equation , and Since they are all constant matrices, it can be deduced that... It is a constant matrix. Therefore, its derivative is... Since the value is 0, the steady-state seepage control equation becomes:

[0142]

[0143] Furthermore, the Riccati equation can be expressed as:

[0144]

[0145] It can be solved by transforming it into an eigenvalue problem of Schur decomposition. First, an intermediate variable is introduced. .

[0146] Its second-order differential equation (with n unknown head functions and nodal fluxes) is transformed into a first-order ordinary differential equation (with 2n unknowns), as shown below:

[0147]

[0148] in, The Hamiltonian matrix is ​​expressed as follows:

[0149]

[0150] Applying Schur decomposition to the Hamiltonian matrix :

[0151]

[0152] in, The eigenvector matrix, It is the eigenvalue matrix.

[0153] By applying Schur decomposition, the equation can be expressed as:

[0154]

[0155] The seepage matrix in an infinite domain is calculated using the following formula:

[0156]

[0157] Furthermore, in order to analyze the seepage field in an infinite domain, the seepage matrix is ​​first... As an initial step, it is then solved iteratively using RK-4 until the seepage matrix at the boundary line is reached. The iterative steps are given by the following formula:

[0158] ;

[0159] ;

[0160] ;

[0161] ;

[0162] in It's the step length. , , , Corresponding to , , , The gradient at that point. It is the gradient at the current point. and It is an estimate based on the intermediate point of the previous gradient. This is the gradient at the next point. By taking a weighted average of these gradients, this method provides a highly accurate estimate of the next solution point.

[0163] This embodiment details how, by acquiring foundation parameters and boundary conditions, a seepage mathematical model is established in a Cartesian coordinate system. An improved proportional boundary finite element coordinate system mapping relationship is constructed by introducing similar lines parallel to the far-field boundary. One-dimensional high-order spectral elements are used to discretize and derive the variable-coefficient control equations. By analyzing the characteristics of the coefficient matrix, the equations are simplified to the Riccati equations. The infinite seepage matrix is ​​solved using Schur decomposition, and then, using this as an initial value, the seepage matrix at the boundary is obtained through iteration using the fourth-order Runge-Kutta method. Finally, the steady-state head field distribution is obtained, ensuring the continuity of the interface head and flow boundary conditions, and significantly improving the accuracy and efficiency of seepage solutions for complex foundations. This provides reliable technical support for engineering seepage prevention design and stability assessment.

[0164] Based on Example 1, this example details the verification of an efficient numerical solution method for the steady-state head field of a dam with complex foundations, specifically as follows:

[0165] Seepage field analysis was performed on a rectangular soil channel, and the calculation results of this invention were compared with those of GeoStudioSEEP / W software (using tens of thousands of elements); the equipotential lines of the rectangular channel were compared, and the two were highly consistent. A comparison of water head values ​​at different depths showed an average relative error of less than 1.5%, demonstrating the high precision of the invention. The seepage field of a rectangular channel was studied; the permeability coefficients of the channel lining and the underlying medium were assumed to be equal. The boundary conditions for the finite and infinite boundaries of the model (assuming they do not penetrate) are specified as zero flux, as follows: Figure 3 As shown, the geometric parameters of the rectangular channel are width ,deep (Water Head) The rectangular river channel boundary line is discretized using 12 two-dimensional third-order line spectrometers based on similarity lines, such as... Figure 3 As shown, the RK-4 parameters are selected as follows: initial radial position Reverse iteration step size Iterated 2000 times until The seepage potential lines in a trapezoidal channel were calculated using the method presented in this paper, and compared with the SEEP / W method of Salmasi and Abraham (2020); the results are as follows. Figure 4 As shown.

[0166] The comparison results from the examples show that the calculation results of this method are basically consistent with the equipotential line orientation of the SEEP / W simulation, exhibiting a high degree of agreement. Furthermore, to more clearly demonstrate the differences between the two calculation results, the horizontal midpoint of the model was taken. (The bottom center point is the origin of the coordinate system). The water head and flow velocity at different depths in the longitudinal direction of the medium are analyzed, as shown in Table 1.

[0167] Table 1 Different water depths in different channels Time head comparison

[0168] As shown in Table 1, the results obtained using the proposed method are in good agreement, with an average error of less than 1.5% from the SEEP / W numerical simulation report. These errors mainly arise from numerical discretization, algorithmic approximation, and boundary condition handling. However, it is worth noting that the SEEP / W numerical model used 19,450 to 51,674 elements, while the SBFEM model used only 10 to 15 elements along the boundary line, combined with 2,000 iterations. Therefore, the proposed method significantly reduces the degrees of freedom (DoF) by 490, highlighting the computational efficiency of the modified SBFEM model.

[0169] To further verify its applicability to practical problems, the proposed method was applied to the seepage problem of a dam-foundation-reservoir system. The model is divided into near-field and far-field, such as... Figure 5As shown, Figure 5 (a) The area enclosed by the red rectangle is the near field. The geometric parameters of the model are: dam height L = 80m, upper base width d = 20m, lower base width W = 60m. The upstream and downstream sides of the dam are water bodies, the foundation is a finite domain, and the far field is the region extending from AEFB to infinity, i.e., an infinite domain foundation. The upstream water level range is AC = 80m, and the water level h u =60m; downstream water level range BD=80m, downstream water level h d =20m, the permeability coefficients of the dam body and foundation are respectively... and It is important to note that the near-field dam body experiences unpressurized seepage, while the far-field foundation experiences pressurized seepage. The far-field boundary conditions are determined by the continuity of head and flow rate at the near-field and far-field interfaces. Figure 5 (b) shows the unpressurized seepage flow, with the following boundary conditions:

[0170]

[0171] In the formula, hi represents the head calculated in each iteration. The iteration stopping criterion for calculating the free surface position in this study is based on the following criteria:

[0172]

[0173] Where N is the number of iterations. Let represent the head of the free surface at the nth iteration, and It is the tolerance error.

[0174] The specific iterative steps are:

[0175] (1) Assume a free surface for seepage, determine its elevation at a specific location, and define the computational domain;

[0176] (2) Apply Neumann boundary conditions to the assumed free surface and calculate the head value on the free surface;

[0177] (3) Compare the head value obtained in step (2) with the assumed elevation to check whether the conditions are met. If the conditions are not met, the calculated head h is used as the new free surface elevation hi, and the computational domain is updated.

[0178] (4) Repeat the above steps until the stopping criteria are met.

[0179] As mentioned earlier, the near field was calculated using a traditional SBFEM, discretized into 42 two-dimensional third-order line elements. The far field was calculated using the method proposed in Section 3, where 22 line elements were selected based on the parameters of the scale line RK-4 as follows: initial radial position The step size of the reverse iteration The process was iterated a total of 2000 times until it reached... .

[0180] The seepage potential curves of the dam-foundation-reservoir system obtained by the finite element method (ANSYS APDL) were compared with those obtained by the proposed method. The results show that the equipotential lines are highly consistent with the overall trend of the flow network, and the proposed method can effectively capture the seepage behavior of the system. To more clearly and intuitively compare the two methods, the hydraulic head values ​​at different vertical distances (x=0) and different horizontal coordinates (y=-20m) were used as indicators. The results are shown in [Figure / Reference]. Figure 6 As shown, the results of the two methods are highly consistent, with relative errors of 1.70% and 1.97%, respectively. This further verifies the accuracy of the proposed method, which uses only 64 line elements and 2000 iterations, significantly reducing computational cost compared to 79653 elements in ANSYS APDL, demonstrating the computational efficiency of this method.

[0181] The above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. For those skilled in the art, the present invention can have various modifications and variations. Any changes, modifications, substitutions, integrations, and parameter changes made to these embodiments within the spirit and principles of the present invention, without departing from the principles and spirit of the present invention, through conventional substitutions or to achieve the same function, fall within the scope of protection of the present invention.

Claims

1. An efficient numerical solution method for the steady-state head field of a dam with complex foundations, characterized in that, include: The medium characteristic parameters, boundary conditions, and geometric parameters of complex foundations are obtained. The medium characteristic parameters include permeability coefficient, porosity, and compressibility coefficient. The boundary conditions include head boundary, flow boundary, and infinity boundary. A mathematical model of steady-state seepage in a half-space is established in a Cartesian coordinate system. The governing equations are constructed based on Darcy's law and Laplace's equation. The Galerkin weighted residual method is used to weaken the function continuity requirement. The characteristics of the half-infinite domain are handled by Fourier transform. By introducing similar lines parallel to the far-field boundary, a mapping relationship between the Cartesian coordinate system and the improved scaled boundary finite element coordinate system is established, and the differential operator transformation is realized through the Jacobian matrix. In the improved proportional boundary finite element coordinate system, one-dimensional high-order spectral elements are used to discretize the boundary lines and similarity lines, and the variable coefficient steady-state seepage control equation is derived based on the weighted residual method and integral by parts. By analyzing the asymptotic properties of the coefficient matrix at infinity, the governing equations are simplified to the Riccati equations. The Hamiltonian matrix is ​​introduced and the steady-state seepage matrix at infinity is solved by Schur decomposition. Using the steady-state seepage matrix at infinity as the initial value, the fourth-order Runge-Kutta method is used to iterate backward along the radial coordinate to obtain the steady-state seepage matrix at the boundary, and then the steady-state head field distribution of the dam is obtained.

2. The efficient numerical solution method for the steady-state head field of dams with complex foundations according to claim 1, characterized in that, The governing equation is: , ,in, It is the source per unit volume. It is the main pressure head. This is the derivative of the total pressure head with respect to time. It is a specific storage coefficient. It is the density of the fluid. It is gravitational acceleration. It is the porosity of the medium. and These are the compressibility coefficients of the solid and the fluid, respectively. It is the exhaust velocity vector. It is a penetration matrix. and respectively along the medium and Permeability coefficient in directional direction, in isotropic media In anisotropic media .

3. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, By transforming the governing equations to the frequency domain using Fourier transform, we obtain: , ,exist Above, among which, , and They are , and Fourier transform; The boundary conditions are specifically as follows: ,exist superior; ,exist Above, among which, For the entire boundary, The water head boundary, For the velocity boundary, To specify the head value, Specify the flow rate value.

4. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The equation for the mapping relationship is: , ,in, and The value of a point in the problem domain. and To scale the coordinates of the splicing line, and The coordinates of the boundary line, To improve the radial coordinates of the scaled boundary finite element coordinate system; The coordinates of nodes in boundaries and similar lines are given by the shape function formula: , ,in, The shape function of the nodal line element. and In line unit Cartesian coordinates of the nodes.

5. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The derivative transformation between the Cartesian coordinate system and the improved scaled boundary finite element coordinate system is based on the Jacobian matrix, the determinant of which is: The derivative transformation relationship is: ,in, To improve the circumferential coordinates of the proportional boundary finite element coordinate system.

6. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The field variables at any point within the analysis domain are obtained through interpolation: ,in, The water head of the points on the ray connecting the nodes of similar line elements and the nodes of the boundary line elements is the radial coordinate. The function; The variable-coefficient steady-state seepage control equation is: ,in, , , These depend on radial coordinates. The coefficient matrix, right The partial derivatives, The transpose of the matrix, right The partial derivatives, right The second-order partial derivative, right The first-order partial derivative, The nodal flux magnitude generated for the normal flux magnitude on the line connecting similar lines and artificial boundaries. The amplitude is the internal source amplitude.

7. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The elements of the coefficient matrix are expressed in rational fraction form as follows: , , ,in , , They are respectively , and The Middle Line number Column elements, For the row and column indices of the matrix elements, It is about polynomial, , , These are the polynomials corresponding to the coefficient matrices, with the first subscript indicating the polynomial's... The highest power, the second subscript corresponds to the coefficient matrix. It is the determinant of the Jacobian matrix, and is about... A first-order polynomial.

8. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The Riccati equation is: ,in, The steady-state seepage matrix at infinity; for The inverse matrix; By introducing intermediate variables Transform the second-order differential equation into a first-order ordinary differential equation: ,in, It is a Hamiltonian matrix. This is the head-related vector. This is the node flux-related vector. right The first-order partial derivative, For Hamiltonian matrices; Schur decomposition yields ,in The eigenvector matrix, for The inverse matrix, For an eigenvalue matrix, select eigenvector submatrices corresponding to eigenvalues ​​with negative real parts. , ,get , for The inverse matrix.

9. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The iterative steps of the fourth-order Runge-Kutta method are as follows: ； ； ； ； ;in, It is the iteration step size. , , , They are respectively , , , Gradient of position, For the first The radial coordinate values ​​of the next iteration For the first The radial coordinate values ​​of the next iteration For the first The seepage matrix obtained in the second iteration For the first The seepage matrix obtained in the second iteration is iterated to... hour, ] represents the steady-state seepage matrix at the boundary.

10. The efficient numerical solution method for the steady-state head field of a dam with complex foundations according to claim 1, characterized in that, The one-dimensional higher-order spectral unit adopts Legendre polynomial shape functions, which are suitable for semi-infinite domain steady-state seepage analysis scenarios with complex foundations containing multiple layers, heterogeneity and anisotropy, such as dam foundations, embankments, and underground engineering. The parallel setting of similar lines and far-field boundaries ensures the continuity of interface head and flow boundary conditions. Combined with the solution process, it realizes accurate prediction and efficient calculation of steady-state head field in engineering practice.

Images