Method for constructing groundwater flow field based on irrotational condition and related equipment
By constructing a groundwater flow field based on irrotational conditions, the problem of flow field analysis deviation in heterogeneous aquifers in existing technologies is solved, achieving efficient and automated flow field simulation and generating reasonable groundwater flow field morphology.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHINA UNIV OF MINING & TECH (BEIJING)
- Filing Date
- 2025-11-17
- Publication Date
- 2026-06-02
AI Technical Summary
Existing spatial interpolation methods struggle to satisfy the physical laws of water flow when dealing with complex aquifers with high heterogeneity, leading to deviations in flow field analysis results and reduced reliability.
The groundwater flow field construction method based on the irrotation condition is to construct a groundwater flow model, perform grid division, use a linear interpolation algorithm to determine the nodal head data, and determine the target head data and flow velocity through multiple rounds of iterative operations. Combined with the physical laws of groundwater movement, it satisfies the refraction law of water flow.
It significantly reduces dependence on hard-to-obtain hydrogeological parameters, automates the simulation of head flow field, improves computational efficiency, is suitable for areas with sparse observation data or small-to-medium scale exploration, and generates flow field morphologies with reasonable physical meaning.
Smart Images

Figure CN122133538A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of groundwater flow field analysis technology, and in particular to a method and related equipment for constructing groundwater flow fields based on irrotational conditions. Background Technology
[0002] In the field of hydrogeological research and engineering practice, hydraulic head, as a core physical indicator characterizing the energy state of groundwater, has always played an irreplaceable and crucial role. Currently, the technical approaches for reconstructing regional groundwater flow fields based on discrete water level observation well data mainly focus on two major directions: spatial interpolation methods and numerical simulation methods. Among them, spatial interpolation methods occupy an important position in practical applications due to their ease of operation.
[0003] However, while existing spatial interpolation methods are widely used in constructing head fields for homogeneous or weakly heterogeneous aquifers, they still have inherent limitations that are difficult to overcome when dealing with complex aquifers with high heterogeneity. Summary of the Invention
[0004] In view of this, the purpose of this application is to propose a groundwater flow field construction method and related equipment based on irrotational conditions, so as to solve the problem that the existing spatial interpolation methods have limitations when dealing with complex aquifers with high heterogeneity.
[0005] To achieve the above objectives, the first aspect of this application provides a method for constructing a groundwater flow field based on irrotational conditions, comprising: A groundwater flow model is constructed, and the internal space of the groundwater flow model is divided into a grid to obtain a computational grid that includes multiple triangular elements and multiple nodes. Based on the head data of the known observation nodes in the groundwater flow model, a linear interpolation algorithm is used to determine the head data of each node other than the observation nodes. Based on the head data of each node, the target head data of each node and the target flow velocity corresponding to each triangular unit are determined through multiple rounds of iterative operations.
[0006] Based on the same inventive concept, a second aspect of this application also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable by the processor, wherein the processor implements the method described above when executing the computer program.
[0007] Based on the same inventive concept, a third aspect of this application also provides a non-transitory computer-readable storage medium storing computer instructions for causing a computer to perform the method described above.
[0008] Based on the same inventive concept, a fourth aspect of this application also provides a computer program product, including computer program instructions that, when run on a computer, cause the computer to perform the method described above.
[0009] As described above, the groundwater flow field construction method and related equipment based on irrotational conditions provided in this application include: constructing a groundwater flow model; dividing the internal space of the groundwater flow model into a grid to obtain a computational grid including multiple triangular elements and multiple nodes; determining the head data of each node other than the observed nodes using a linear interpolation algorithm based on the head data of known observation nodes in the groundwater flow model; and determining the target head data of each node and the target flow velocity corresponding to each triangular element through multiple iterations based on the head data of each node. This application uses irrotational conditions as constraints to ensure that the flow field morphology conforms to the groundwater movement law, and satisfies the refraction law of water flow at the boundary of regions where the properties of the aquifer change abruptly. It significantly reduces the dependence on difficult-to-obtain hydrogeological parameters such as storage coefficients and boundary conditions; the core inputs are only the head data of the observed nodes and the element permeability coefficients. It transforms the traditional complex model generalization process into programmable iterative calculation steps, automating the head flow field simulation, making it suitable for integration into lightweight tools, and improving computational efficiency. Even in areas with sparse observation data or in small- to medium-scale surveys, it can still quickly output physically reasonable flow field morphology, demonstrating strong applicability. Attached Figure Description
[0010] To more clearly illustrate the technical solutions in this application or related technologies, the drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0011] Figure 1 This is a flowchart illustrating the groundwater flow field construction method based on irrotational conditions according to an embodiment of this application. Figure 2 This is a schematic diagram of a groundwater flow model according to an embodiment of this application; Figure 3 for Figure 2 Corresponding computational grid diagram; Figure 4 This is a schematic diagram of a triangular unit according to an embodiment of this application; Figure 5 This is a schematic diagram showing two triangular units sharing a side in an embodiment of this application; Figure 6 This is a schematic diagram of a sub-region in an embodiment of this application; Figure 7 This is a schematic diagram of the flow velocity calculation for the triangular element in an embodiment of this application; Figure 8 This is a schematic diagram of multiple sub-regions in an embodiment of this application; Figure 9 This is a schematic diagram illustrating head data correction in an embodiment of this application; Figure 10 This is a schematic diagram of the flow field of the flow around the model in an embodiment of this application; Figure 11 This is a schematic diagram of the head field of the flow around the model in an embodiment of this application; Figure 12 This is a schematic diagram of a groundwater flow model according to another embodiment of this application; Figure 13 for Figure 12 Corresponding computational grid diagram; Figure 14 This is a schematic diagram of the flow field of the source-sink model in an embodiment of this application; Figure 15 This is a schematic diagram of the head field of the source-sink model in an embodiment of this application; Figure 16 This is a schematic diagram of the structure of the groundwater flow field construction device based on the irrotation condition according to an embodiment of this application; Figure 17 This is a schematic diagram of the hardware structure of an electronic device according to an embodiment of this application. Detailed Implementation
[0012] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with specific embodiments and the accompanying drawings.
[0013] It should be noted that, unless otherwise defined, the technical or scientific terms used in the embodiments of this application should have the ordinary meaning understood by one of ordinary skill in the art to which this application pertains. The terms "first," "second," and similar terms used in the embodiments of this application do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are only used to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.
[0014] Kriging, a classic geostatistical interpolation method based on the principle of spatial autocorrelation, uses the variogram as its core analytical tool. It determines the optimal weighting coefficients by quantifying the spatial structural characteristics and random fluctuations of water level data, achieving the optimal unbiased estimate of water level values at unsampled points. This method can deeply mine the spatial variation patterns hidden in water level data, generating high-precision head field results that conform to the theoretical framework of geostatistics, providing fundamental support for flow field analysis. Inverse Distance Weighted (IDW), based on the core principle of "distance decay," performs interpolation by calculating the inverse weighted average of the distances between the point to be estimated and known observation points. It is suitable for groundwater level datasets with significant spatial autocorrelation, reflecting both the spatial continuity of water levels and capturing local variation information. It also boasts advantages such as high computational efficiency, simple implementation process, and strong adaptability to data distribution, demonstrating good robustness in hydrogeological practice. Discrete Smooth Interpolation (DSI) is an interpolation technique developed based on the least squares principle and constrained optimization strategy. By constructing a global optimization model and using an iterative solution algorithm, the interpolation results can strictly fit known water level observations and retain the authenticity of measured data. At the same time, through smoothness constraints, it can reasonably reproduce the spatial gradient characteristics of the water level field, effectively avoid meaningless local fluctuations, and accurately adapt to the complex spatial distribution law of groundwater level controlled by geological structure.
[0015] Although the aforementioned spatial interpolation methods are widely used in constructing hydraulic head fields for homogeneous or weakly heterogeneous aquifers, they still have inherent limitations that are difficult to overcome when dealing with complex aquifers with high degrees of heterogeneity. Especially in areas where the properties of the aquifer medium change abruptly, these methods rely solely on the spatial statistical characteristics of the data for interpolation, failing to couple with the physical mechanisms of groundwater movement and thus failing to satisfy the laws of water flow, such as the law of refraction. This deficiency directly leads to the constructed flow field violating the physical laws of water flow, resulting in a deviation between the constructed flow field and the actual flow field, significantly reducing the reliability of the flow field analysis results.
[0016] In view of this, this application proposes a groundwater flow field construction method based on irrotational conditions, which can deeply couple the physical mechanism of groundwater movement and take into account the computational efficiency, physical realism of the flow field and adaptability to complex geological structures.
[0017] The embodiments of this application will be described in detail below with reference to the accompanying drawings.
[0018] The embodiments of this application provide a method for constructing a groundwater flow field based on irrotational conditions, referencing... Figure 1 This includes the following steps: Step 101: Construct a groundwater flow model and divide the internal space of the groundwater flow model into a grid to obtain a computational grid that includes multiple triangular elements and multiple nodes.
[0019] Specifically, Figure 2 The groundwater flow model in this embodiment is shown; this groundwater flow model is a bypass model. The study area (internal space) of the bypass model has a horizontal and vertical length of 100m. It has a constant head boundary at the top (100m) and a constant head boundary at the bottom (90m), with impermeable boundaries on the left and right sides. The study area has two weakly permeable zones, the locations and dimensions of which are shown below. Figure 2 As shown. The flow velocity in the weakly permeable region is K=0.00001m / d, and the flow velocity in other regions is K=1m / d. Based on the outer boundary of the study area, the boundary of the internal weakly permeable region, and the location of the observation wells, the study area of the flow model is divided into Delaunay triangular meshes as constraints, generating the following... Figure 3 The computational mesh shown is composed of multiple triangular cells and multiple nodes. Each vertex of each triangular cell is a node. Figure 3 The CCP consists of 1579 triangular units and 839 nodes. For example... Figure 3 As shown, observation nodes are represented by solid circles, indicated in red and blue respectively. The head data for each node refers to the head data corresponding to the vertices of the triangular mesh after dividing the study area into triangular meshes. For example, the number of observation nodes can be 10. The head data of the observation nodes corresponding to the observation wells are obtained by running the flow model using the numerical simulation software Feflow. The head data of the observation nodes are considered known data.
[0020] Step 102: Based on the head data of the observation nodes, use a linear interpolation algorithm to determine the head data of each node other than the observation nodes.
[0021] Based on the head data from the observation nodes, the head data for all nodes other than the observation nodes are calculated using the Kriging linear interpolation algorithm. The Kriging linear interpolation algorithm provides an initial head field that is spatially smooth and conforms to the known head data from the observation nodes.
[0022] Step 103: Based on the head data of each node, determine the target head data of each node and the target flow velocity corresponding to each triangular unit through multiple rounds of iterative operations.
[0023] Further, step 103 includes: Based on the node coordinates, the head data of each node, and the preset unit permeability coefficient, the velocity component within each triangular unit is determined.
[0024] Specifically, the velocity components within each triangular unit and It can be expressed by the following formula (1): (1) in, Represents the area of the triangular unit. This represents the preset permeability coefficient of the triangular unit. , and These represent the head data of the three nodes of the triangular element. This represents the x-component of the flow velocity within the triangular unit. This represents the y-component of the flow velocity within the triangular unit.
[0025] Perform multiple rounds of iterative operations, with each round of iterative operations as follows: The flow rate value corresponding to each side of the triangular unit is determined based on the flow velocity component within the triangular unit.
[0026] Specifically, based on the calculated velocity components within the triangular unit... and Determine the flow rate corresponding to each side within the triangular cell. Figure 4 A schematic diagram of a triangular element is shown. The three nodes of the triangular element are e, f, and g. The velocity vector of the triangular element is... The velocity components in the x and y directions are and Given the coordinates of the three nodes e, f, and g, what is the flow rate through the edge fg? It can be calculated using the following formula (2): (2) in, and These are the components of edge fg in the x and y directions. , , ( , ) are the node coordinates of node g, ( , ) are the node coordinates of node f. The flow value corresponding to each side in the triangular cell can be calculated using equation (2).
[0027] It should be noted that if an edge is shared by two triangular units, the flow value corresponding to that edge also needs to be corrected. Figure 5 A schematic diagram showing two triangular units sharing a side is shown, as follows. Figure 5As shown, the left and right triangular units share a common edge. The flow rate calculated through the left triangular unit along the shared edge is... The flow rate of the shared side calculated using the triangular unit on the right is... . and The velocity is calculated according to formula (2). The velocity corresponding to the left triangle is... The flow velocity corresponding to the right triangle is .
[0028] Then the flow value corresponding to the corrected shared edge It can be calculated using the following formula (3): (3) in, and They represent Figure 5 The unit permeability coefficients corresponding to the two triangular units in the diagram.
[0029] Based on the common nodes of the triangular units, the computational grid is divided into multiple sub-regions. For each sub-region, the first flow velocity corresponding to each triangular unit is determined based on the flow value corresponding to each edge within the triangular unit.
[0030] Specifically, in cases where multiple triangular cells share a node in the computational grid, the computational grid is divided into multiple sub-regions based on the situation of multiple triangular cells sharing a node. Figure 6 A schematic diagram of the sub-region is shown. Figure 6 In this example, five triangular units share a single node j. The flow value of the edge opposite node j is... The flow values of the sides corresponding to node j in the five triangular units can be calculated using formula (3). , , , and .
[0031] Furthermore, based on the flow rate value corresponding to each side within the triangular unit, the first flow velocity corresponding to each triangular unit is determined, including: The flow value corresponding to the triangular unit is determined based on the flow value corresponding to each side within the triangular unit.
[0032] Specifically, based on the principle of mass conservation, for the case where the sum of the flows to opposite edges of a node is not zero, the formula for allocating the flow to each triangle as a source and sink based on the area of the triangular units within the sub-region is as follows: (4) in, This represents the flow rate allocated to the i-th triangular cell within the sub-region, which is also the flow rate value corresponding to the triangular cell. This represents the preset cell permeability coefficient of the i-th triangular cell within the sub-region. Let represent the area of the i-th triangular unit within the subregion. The flow value of the edge corresponding to node j in the i-th triangular cell within the sub-region.
[0033] Based on the flow value corresponding to each edge within the triangular unit and the flow value corresponding to the triangular unit itself, the flow value of the inner boundary within each sub-region is determined by calculation, including: Based on the flow value corresponding to each side within the triangular unit and the flow value corresponding to the triangular unit, an effective equation for the flow value of the inner boundary within the sub-region is constructed using the principle of mass conservation and the irrotation condition. The effective equation is then solved to determine the flow value of the inner boundary within each sub-region.
[0034] Specifically, by utilizing the mass conservation condition of each triangle within the sub-region, effective equations regarding the flow of the triangle edges within the sub-region can be established with respect to the number of triangles *n*. The number of effective equations is *n-1*. According to the principle of mass conservation, the total inflow to each triangular unit equals the total outflow. Figure 6 For each flow direction shown, establish an effective equation for the flow along the triangular edges within the sub-region, and express it as follows: (5) in, This represents the flow of each edge opposite node j. This represents the flow rate allocated to each triangular cell. This represents the flow (belonging to the flow to be calculated) around the inner boundary of node j within the sub-region. Figure 6 middle, The preceding minus sign indicates The flow direction is negative.
[0035] Based on the irrotational condition, i.e., the groundwater hydraulic gradient within the subregion is irrotational, an effective equation for the triangular side flow within the subregion is established, expressed by the following formula: (6) Where, when i takes the value n, The value is 1. n is the total number of triangular units within the subregion. Let represent the side length of the opposite side of node j within the i-th triangular cell in the subregion. This represents the area of the i-th triangular unit in the subregion.
[0036] Based on formulas (5) and (6), a system of equations is constructed to obtain the flow rates of the inner boundaries of all triangular units within the subregion. This yields the flow value for each edge within each triangular cell of the subregion. Figure 6 The first triangular cell in the array has flow values for each side as follows: , and .
[0037] The first flow velocity corresponding to each triangular unit is determined based on the flow value of the inner boundary of each sub-region and the flow value corresponding to each side within the triangular unit.
[0038] Specifically, based on the flow rate value of the inner boundary within each sub-region and the flow rate value corresponding to each edge within the triangular unit, the first flow velocity corresponding to each triangular unit is calculated according to the Raviart-Thomas space. Figure 7 A schematic diagram illustrating the flow velocity calculation for a triangular element is shown. Figure 7 middle, , and This is a vector representing the direction and magnitude of the flow. When the flow q flows into the triangle along one edge, the defined vector is the vector pointing from the centroid of the triangle to the vertex opposite that edge. When the flow q flows out of the triangle, the defined vector is the vector pointing from the vertex opposite that edge to the centroid of the triangle. , and The direction and flow direction of flow q are as follows Figure 7 As shown. The first flow velocity corresponding to the triangular element is calculated based on the Raviart-Thomas space. It can be expressed by the following formula: (7) The first flow velocity is corrected based on the division of multiple sub-regions to obtain the initial flow velocity corresponding to each triangular unit.
[0039] It should be noted that triangular units can be shared, such as... Figure 8 As shown, the same triangular unit can be shared by a maximum of three sub-regions. Figure 8 In this example, nodes e, f, and g each correspond to a sub-region. Nodes e, r, and g form a shared triangular unit. For each shared triangular unit, the first velocity of that unit is calculated for each sub-region. The average of all first velocities is then taken to obtain the initial velocity corresponding to the shared triangular unit. For example, if three sub-regions share a triangular unit, the first velocities calculated for each sub-region are as follows: , and Then the initial flow velocity of the triangular unit = ( + ) / 3.
[0040] Based on the initial flow velocity corresponding to each triangular unit and the head data of each node, the head data of each node other than the observed node is corrected to obtain the corrected head data.
[0041] For each triangular cell within each sub-region, based on the initial flow velocity corresponding to the triangular cell and the head data of the nodes within the triangular cell, the corrected sub-head data of the shared nodes within the sub-region corresponding to the triangular cell are determined; after weighted averaging of all sub-head data within the sub-region, the corrected head data of the shared nodes is obtained.
[0042] Specifically, Figure 9 The diagram illustrates head data correction. The central node is 'a', and other nodes include b, c, d, e, and f. The sub-region comprises five triangular units, numbered 1, 2, 3, 4, and 5. The corrected sub-head data for the shared node 'a' within the sub-region corresponding to each triangular unit can be expressed by the following formula: (8) in, This represents the sub-head data of the first triangular cell within the sub-region relative to the center node a. This represents the head data for node b. This represents the head data for node c. This represents the initial flow velocity of the first triangular unit. and Let represent the edge vector. The sub-head data of other triangular elements within the sub-region relative to the central node a can be calculated according to formula (8), which will not be elaborated here. After calculating the sub-head data of all triangular elements relative to the central node a, a weighted average is performed using the following formula to obtain the corrected head data of the shared node a. (9).
[0043] For each node, determine whether the difference between the corrected head data and the original head data is less than a preset threshold. If not, update the velocity component within the triangular unit for the next iteration based on the initial velocity of the triangular unit, and execute the next iteration. If yes, exit the multi-round iteration, and use the corrected head data obtained in the current round as the target head data and the initial velocity obtained in the current round as the target velocity.
[0044] Specifically, after calculating the corrected head data for each node, compared with the uncorrected head data for that node, if the difference between the corrected and uncorrected head data is greater than a preset threshold, it indicates that the difference between the corrected and uncorrected head data is large and does not meet the convergence condition. Therefore, the next iteration is performed, updating the velocity components within the triangular unit for the next iteration based on the initial velocity of the triangular unit. + ) / 3 is decomposed into the next round of iteration operations. and Starting from the above formula (2), the initial flow velocity and the corrected head data are calculated again until the difference between the corrected head data and the uncorrected head data is less than a preset threshold. Then, the multi-round iteration operation is exited, and the corrected head data obtained in the current round is taken as the target head data, and the initial flow velocity obtained in the current round is taken as the target flow velocity. The converged groundwater flow field is obtained, including the head field and the velocity field. The head field is composed of the target head data, and the velocity field is composed of the target flow velocity.
[0045] In some embodiments, it also includes: Based on the target head data of each node and the target flow velocity corresponding to each triangular unit, the head data corresponding to each coordinate point in the internal space is calculated and determined.
[0046] Specifically, after obtaining the converged groundwater flow field through the aforementioned embodiments, the hydraulic head data at any coordinate point in the study area (internal space) is calculated based on the linear distribution of hydraulic head within the triangle. Each triangular unit consists of three nodes e, f, and g, with corresponding hydraulic head data as follows: The head data at any point is calculated using linear interpolation and is expressed by the following formula: (10) in, The weighting coefficient function for node e is represented by... The weighting coefficient function represents node f. The weighting coefficient function represents node g. The calculation formula is as follows:
[0047]
[0048]
[0049] in, Represents the area of the triangular unit, ( , ) represents the node coordinates of node e, ( , ) represents the node coordinates of node f, ( , ) represents the node coordinates of node g.
[0050] The flow field of the flow model obtained by the method in the foregoing embodiments is as follows: Figure 10 As shown, the head field is as follows Figure 11 As shown. Figure 10 The legend on the right side of the middle figure shows the different flow rates represented by different colors. Figure 11 The legend on the right side of the middle figure uses different colors to represent different head values. From Figure 10 and Figure 11 As can be seen, the flow field clearly demonstrates the flow around the groundwater when it encounters a weakly permeable area. The water flows around the two weakly permeable areas, and at the boundary of the weakly permeable areas, the flow pattern conforms to the law of refraction, which is highly consistent with the physical laws of groundwater movement and the results of traditional numerical simulations.
[0051] In another example, Figure 12 A groundwater flow model according to another embodiment of this application is shown, which is a source-sink model. The study area of the source-sink model has a constant head boundary of 100 m on the left, a constant head boundary of 90 m on the right, and impermeable boundaries on the top and bottom. The study area has one weakly permeable zone and one strongly permeable zone, the locations of which are as follows... Figure 12 As shown. The flow velocity in the poorly permeable zone is 0.1 m / d, and the flow velocity in the highly permeable zone is 10 m / d. The flow velocity in all other zones (excluding the poorly and highly permeable zones) is 1 m / d. Figure 12 In the study, five points were equidistantly selected at the left boundary of the fixed head, with a head of 100 m, and five points were equidistantly selected at the right boundary of the fixed head, with a head of 90 m. The point at coordinates (50, 50) was designated as the source-sink point with a head of 91 m. Based on the outer boundary of the study area, the boundaries of the internal weakly permeable and strongly permeable regions, and the locations of the observation wells, local mesh refinement was performed near the source-sink points as constraints. The study area was then divided into Delaunay triangular meshes, and the final generated computational mesh is shown below. Figure 13 As shown, it contains 1913 triangular units and 1008 nodes, with observation nodes represented by red and blue solid dots. Figure 12 The groundwater flow model shown is consistent with Figure 2 The groundwater flow field construction method of the groundwater flow model shown is basically the same. The difference is that in this embodiment, the first flow velocity of the triangular unit in the sub-region where the source and sink are located is not calculated, and the variable flow rate of the triangular unit in the sub-region where the source and sink are located is set to zero. Figure 14 It shows the relationship with Figure 12 The corresponding source-sink model flow field diagram, Figure 15 It shows the relationship with Figure 12 The corresponding water head field diagram. Figure 14 and Figure 15 The flow field clearly demonstrates the characteristics of the flow field under the condition of source and sink at the observation well. The flow velocity originates from the constant head boundary on the left and converges towards the low head observation well. The head field forms a distinct drop funnel at the source and sink. Furthermore, when the water flows through media with different permeability, its path is deflected according to the law of refraction. It converges when entering the low-permeability zone and moves away when entering the high-permeability zone, forming a guiding effect similar to a lens. The method in this embodiment successfully incorporates point-like source and sink terms as constraints into the iterative calculation framework, ultimately obtaining a flow field that simultaneously satisfies the law of refraction, mass conservation, and certain source and sink conditions.
[0052] The method described in this application can effectively handle different hydrogeological conditions, such as abrupt changes in permeability coefficients and situations where the flow field has sources and sinks. The groundwater flow field construction method based on irrotational conditions in this application, driven by the physical laws of groundwater movement and based on the water level at observation well points, can quickly and automatically generate physically reasonable flow fields through efficient iterative calculations, providing a powerful new tool for groundwater simulation and analysis.
[0053] It should be noted that the method in this embodiment can be executed by a single device, such as a computer or server. The method can also be applied in a distributed scenario, where multiple devices cooperate to complete the task. In such a distributed scenario, one of these devices may execute only one or more steps of the method in this embodiment, and the multiple devices will interact with each other to complete the method described.
[0054] It should be noted that some embodiments of this application have been described above. In some cases, the actions or steps described in the above embodiments can be performed in a different order than that shown in the above embodiments and the desired result can still be achieved. In addition, the processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.
[0055] Based on the same inventive concept, corresponding to any of the above embodiments, this application also provides a groundwater flow field construction device based on irrotational conditions.
[0056] refer to Figure 16 The groundwater flow field construction device based on the irrotation condition includes: The partitioning module 10 is configured to construct a groundwater flow model and perform mesh partitioning on the internal space of the groundwater flow model to obtain a computational mesh including multiple triangular elements and multiple nodes. The first determining module 20 is configured to determine the head data of each node other than the observed nodes based on the head data of the known observation nodes in the groundwater flow model using a linear interpolation algorithm. The second determining module 30 is configured to determine the target head data of each node and the target flow velocity corresponding to each triangular unit through multiple rounds of iterative operations based on the head data of each node.
[0057] In some embodiments, the second determining module 30 is further configured to determine the velocity component within each triangular unit based on the node coordinates of the nodes, the head data of each node, and the preset unit permeability coefficient; and to perform multiple rounds of iterative operations, each round of iterative operations being as follows: The flow rate value corresponding to each side of the triangular unit is determined based on the flow velocity component within the triangular unit. Based on the common nodes of the triangular units, the computing network is divided into multiple sub-regions. For each sub-region, the first flow velocity corresponding to each triangular unit is determined based on the flow value corresponding to each edge within the triangular unit. The first flow velocity is corrected based on the division of multiple sub-regions to obtain the initial flow velocity corresponding to each triangular unit; Based on the initial flow velocity corresponding to each triangular unit and the head data of each node, the head data of each node other than the observed node is corrected to obtain the corrected head data. For each node, determine whether the difference between the corrected head data and the uncorrected head data is less than a preset threshold. If not, update the velocity component within the triangular unit for the next iteration based on the initial velocity of the triangular unit, and execute the next iteration. If yes, exit the multi-round iteration, take the corrected head data obtained in the current round as the target head data, and take the initial velocity obtained in the current round as the target velocity observation node.
[0058] In some embodiments, after determining the flow rate value corresponding to each side of the triangular unit based on the flow velocity component within the triangular unit, the second determining module 30 is further configured to, in response to determining that the side is a shared side of two triangular units, correct the flow rate value of the shared side based on the preset unit permeability coefficient of each of the two triangular units and the flow rate value corresponding to the shared side in the two triangular units respectively. In some embodiments, the second determining module 30 is further configured to determine the flow value corresponding to the triangle unit based on the flow value corresponding to each side within the triangle unit; Based on the flow value corresponding to each side within the triangular unit and the flow value corresponding to the triangular unit, the flow value of the inner boundary within each sub-region is determined by calculation. The first flow velocity corresponding to each triangular unit is determined based on the flow value of the inner boundary of each sub-region and the flow value corresponding to each side within the triangular unit.
[0059] In some embodiments, the second determining module 30 is further configured to construct an effective equation for the flow value of the inner boundary within a sub-region based on the flow value corresponding to each side within the triangular unit and the flow value corresponding to the triangular unit, using the principle of mass conservation and the irrotation condition, and to solve the effective equation to determine the flow value of the inner boundary within each sub-region.
[0060] In some embodiments, the second determining module 30 is further configured to, for each triangular cell in each sub-region, determine the corrected sub-head data of the common nodes in the sub-region corresponding to the triangular cell, based on the initial flow velocity corresponding to the triangular cell and the head data of the nodes within the triangular cell; After performing a weighted average of all sub-head data within the sub-region, the corrected head data of the shared node is obtained.
[0061] In some embodiments, a calculation module is also included, configured to calculate and determine the head data corresponding to each coordinate point in the internal space based on the target head data of each node and the target flow velocity corresponding to each triangular unit.
[0062] For ease of description, the above devices are described in terms of function, divided into various modules. Of course, in implementing this application, the functions of each module can be implemented in one or more software and / or hardware.
[0063] The apparatus described above is used to implement the corresponding groundwater flow field construction method based on irrotational conditions in any of the foregoing embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be repeated here.
[0064] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this application also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the groundwater flow field construction method based on irrotational conditions described in any of the above embodiments.
[0065] Figure 17This embodiment illustrates a more specific hardware structure of an electronic device. The device may include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040, and a bus 1050. The processor 1010, memory 1020, input / output interface 1030, and communication interface 1040 are interconnected internally via the bus 1050.
[0066] The processor 1010 can be implemented using a general-purpose CPU (Central Processing Unit), microprocessor, application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.
[0067] The memory 1020 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 can store the operating system and other applications. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program code is stored in the memory 1020 and is called and executed by the processor 1010.
[0068] The input / output interface 1030 is used to connect input / output modules to realize information input and output. Input / output modules can be configured as components within the device (not shown in the figure) or externally connected to the device to provide corresponding functions. Input devices may include keyboards, mice, touchscreens, microphones, various sensors, etc., while output devices may include displays, speakers, vibrators, indicator lights, etc.
[0069] The communication interface 1040 is used to connect a communication module (not shown in the figure) to enable communication between this device and other devices. The communication module can communicate via wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.).
[0070] Bus 1050 includes a pathway for transmitting information between various components of the device, such as processor 1010, memory 1020, input / output interface 1030, and communication interface 1040.
[0071] It should be noted that although the above-described device only shows the processor 1010, memory 1020, input / output interface 1030, communication interface 1040, and bus 1050, in specific implementations, the device may also include other components necessary for normal operation. Furthermore, those skilled in the art will understand that the above-described device may only include the components necessary for implementing the embodiments of this specification, and not necessarily all the components shown in the figures.
[0072] The electronic devices described above are used to implement the corresponding groundwater flow field construction method based on irrotational conditions in any of the foregoing embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.
[0073] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this application also provides a non-transitory computer-readable storage medium storing computer instructions for causing the computer to execute the groundwater flow field construction method based on irrotational conditions as described in any of the above embodiments.
[0074] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.
[0075] The computer instructions stored in the storage medium of the above embodiments are used to cause the computer to execute the groundwater flow field construction method based on the irrotation condition as described in any of the above embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.
[0076] Based on the same concept, corresponding to any of the above embodiments, this application also provides a computer program product, including computer program instructions, which, when run on a computer, cause the computer to perform the method described in any of the above embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.
[0077] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of this application is limited to these examples; under the concept of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of the embodiments of this application as described above, which are not provided in detail for the sake of brevity.
[0078] Additionally, to simplify the description and discussion, and to avoid obscuring the embodiments of this application, the well-known power / ground connections to integrated circuit (IC) chips and other components may or may not be shown in the provided drawings. Furthermore, the apparatus may be shown in block diagram form to avoid obscuring the embodiments of this application, and this also takes into account the fact that the details of the implementation of these block diagram apparatuses are highly dependent on the platform on which the embodiments of this application will be implemented (i.e., these details should be fully understood by those skilled in the art). While specific details (e.g., circuits) have been set forth to describe exemplary embodiments of this application, it will be apparent to those skilled in the art that the embodiments of this application can be implemented without these specific details or with variations thereof. Therefore, these descriptions should be considered illustrative rather than restrictive.
[0079] Although this application has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed.
[0080] The embodiments of this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of this application. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the embodiments of this application should be included within the protection scope of this application.
Claims
1. A method for constructing a groundwater flow field based on irrotational conditions, characterized in that, include: A groundwater flow model is constructed, and the internal space of the groundwater flow model is divided into a grid to obtain a computational grid that includes multiple triangular elements and multiple nodes. Based on the head data of the known observation nodes in the groundwater flow model, a linear interpolation algorithm is used to determine the head data of each node other than the observation nodes. Based on the head data of each node, the target head data of each node and the target flow velocity corresponding to each triangular unit are determined through multiple rounds of iterative operations.
2. The method according to claim 1, characterized in that, The head data based on each node is used to determine the target head data for each node and the target flow velocity corresponding to each triangular unit through multiple rounds of iterative operations, including: Based on the node coordinates, head data of each node, and preset unit permeability coefficient, the velocity component within each triangular unit is determined; multiple iterations are performed, with each iteration proceeding as follows: The flow rate value corresponding to each side of the triangular unit is determined based on the flow velocity component within the triangular unit. Based on the common nodes of the triangular units, the computing network is divided into multiple sub-regions. For each sub-region, the first flow velocity corresponding to each triangular unit is determined based on the flow value corresponding to each edge within the triangular unit. The first flow velocity is corrected based on the division of multiple sub-regions to obtain the initial flow velocity corresponding to each triangular unit; Based on the initial flow velocity corresponding to each triangular unit and the head data of each node, the head data of each node other than the observed node is corrected to obtain the corrected head data. For each node, determine whether the difference between the corrected head data and the original head data is less than a preset threshold. If not, update the velocity component within the triangular unit for the next iteration based on the initial velocity of the triangular unit, and execute the next iteration. If yes, exit the multi-round iteration, and use the corrected head data obtained in the current round as the target head data and the initial velocity obtained in the current round as the target velocity.
3. The method according to claim 2, characterized in that, After determining the flow rate value corresponding to each side within the triangular unit based on the flow velocity components within the triangular unit, the process includes: In response to determining that the edge is a shared edge of two triangular units, the flow rate value of the shared edge is corrected according to the preset unit permeability coefficient of each of the two triangular units and the flow rate value corresponding to the shared edge in the two triangular units respectively.
4. The method according to claim 2, characterized in that, The step of determining the first flow velocity corresponding to each triangular unit based on the flow value corresponding to each side within the triangular unit includes: The flow value corresponding to the triangle unit is determined based on the flow value corresponding to each side within the triangle unit; Based on the flow value corresponding to each side within the triangular unit and the flow value corresponding to the triangular unit, the flow value of the inner boundary within each sub-region is determined by calculation. The first flow velocity corresponding to each triangular unit is determined based on the flow value of the inner boundary of each sub-region and the flow value corresponding to each side within the triangular unit.
5. The method according to claim 4, characterized in that, The process of calculating and determining the flow value of the inner boundary of each sub-region based on the flow value corresponding to each edge within the triangular unit and the flow value corresponding to the triangular unit itself includes: Based on the flow value corresponding to each side within the triangular unit and the flow value corresponding to the triangular unit, an effective equation for the flow value of the inner boundary within the sub-region is constructed using the principle of mass conservation and the irrotation condition. The effective equation is then solved to determine the flow value of the inner boundary within each sub-region.
6. The method according to claim 2, characterized in that, The method involves correcting the head data of each node (excluding the observed node) based on the initial flow velocity corresponding to each triangular unit and the head data of each node, resulting in corrected head data, including: For each triangular cell within each sub-region, based on the initial flow velocity corresponding to the triangular cell and the head data of the nodes within the triangular cell, the corrected sub-head data of the shared nodes within the sub-region corresponding to the triangular cell are determined; After performing a weighted average of all sub-head data within the sub-region, the corrected head data of the shared node is obtained.
7. The method according to claim 1, characterized in that, Also includes: Based on the target head data of each node and the target flow velocity corresponding to each triangular unit, the head data corresponding to each coordinate point in the internal space is calculated and determined.
8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the method as described in any one of claims 1 to 7.
9. A non-transitory computer-readable storage medium storing computer instructions, characterized in that, The computer instructions are used to cause the computer to perform the method according to any one of claims 1 to 7.
10. A computer program product comprising computer program instructions, characterized in that, When the computer program instructions are executed on a computer, the computer causes the computer to perform the method as described in any one of claims 1-7.