A method for simulating well flow in a three-dimensional finite element unstructured grid

Abstract

The application discloses a kind of well flow simulation methods in three-dimensional finite element unstructured grid, belong to groundwater seepage finite element simulation field, comprising the following steps: S1, determine model calculation basic parameter;S2, based on all basic parameters, construct three-dimensional unstructured grid seepage model and assemble total stiffness matrix;S3, based on the characteristic size of element around well, combined with soil permeability coefficient and pressure-reducing well radius parameter, calculate first round well flow resistance coefficient and solve initial seepage field;S4, based on the spatial position of element around well, layout sampling point around well and interpolation obtains sampling water head;S5, based on each sampling ring average water head, update iteration resistance coefficient and recalculate seepage field;S6, execute iteration convergence determination and output final result.The three-dimensional finite element unstructured grid well flow simulation method described above realizes the original adaptation of arbitrary unstructured grid, and overall considers the calculation efficiency, precision stability and engineering scene universality.

Description

Technical Field

[0001] This invention relates to the field of finite element simulation technology for groundwater seepage, and more particularly to a method for simulating well flow in a three-dimensional finite element unstructured mesh. Background Technology

[0002] In the field of civil engineering, pressure relief wells are a core engineering measure for groundwater seepage control, widely used in anti-buoyancy, anti-seepage, and anti-surge treatment scenarios for various underground structures. From a broader disciplinary perspective, well flow simulation is essentially the simulation of fluid flow within a tubular network embedded in a porous medium. This type of problem is prevalent in multiple engineering and research fields: in the energy sector, the extraction of oil, natural gas, shale gas, and geothermal resources, as well as CO2 formation sealing, all involve seepage coupling between the well casing and the reservoir medium; in the field of bioengineering, the water exchange process between plant roots and soil exhibits similar multi-scale characteristics; in the medical field, the efficiency of tumor thermotherapy, the efficiency of microcirculation drug delivery, and the blood and oxygen transport processes between cerebral blood vessels and brain tissue also fall under the mass transfer coupling problem between tubular channels and the surrounding medium. In these various thermo-hydraulic-mechanical coupling systems, the tubular network is the core channel for mass transfer and the main source of action driving the flow of the entire system. Therefore, accurately simulating the coupling effect between the tubular network and the porous medium is crucial to ensuring the reliability of numerical simulations in these fields.

[0003] All the above problems share common geometric characteristics and mathematical descriptions: the characteristic dimensions of the tubular structure differ from the dimensions of the computational domain of the porous medium by several orders of magnitude, making them typical multi-scale flow problems. If the traditional equal-dimensional finite element method is used for seepage simulation, in order to ensure computational accuracy, the element mesh size around the tubular structure must be divided to be smaller than the characteristic dimensions of the tube itself. This approach leads to a sharp increase in the number of mesh elements, extremely low computational efficiency, and excessive consumption of computational storage and computing power resources. Summary of the Invention

[0004] The purpose of this invention is to provide a method for simulating well flow in a three-dimensional finite element unstructured mesh, thereby solving the aforementioned technical problems.

[0005] To achieve the above objectives, this invention provides a method for simulating well flow in a three-dimensional finite element unstructured mesh, comprising the following steps: S1. Determine the basic parameters for model calculation: Determine the boundary conditions of the three-dimensional seepage calculation domain geometric model, the water head boundary conditions of the pressure relief well, the spatial coordinates of the pressure relief well, and the geometric parameters of the pressure relief well (including the radius of the pressure relief well). ), soil permeability coefficient, soil spatial parameters, and iterative convergence threshold parameters; S2. Constructing a three-dimensional unstructured mesh seepage model and assembling the overall stiffness matrix: Based on all the basic parameters determined in S1, the geometric model of the three-dimensional seepage computation domain is divided into Delaunay meshes. At the same time, the three-dimensional decompression well is geometrically reduced. The well section is discretized into one-dimensional well line elements along the well axis, and the bottom section is discretized into zero-dimensional bottom point elements. The one-dimensional well line elements and zero-dimensional bottom point elements are embedded into the corresponding spatial coordinate positions of the three-dimensional finite element computation mesh of the soil to construct a 3D-1D-0D mixed-dimensional computation mesh. Finally, based on the geometric properties of the three-dimensional finite element computation mesh of the soil and the permeability coefficient parameters of the corresponding soil layer, the three-dimensional seepage overall stiffness matrix is ​​assembled. S3. Calculate the first-round well flow resistance coefficient and solve the initial seepage field: Based on the characteristic dimensions of the well perimeter element obtained in S2, and combined with the soil permeability coefficient and decompression well radius parameters in S1, calculate the initial values ​​of the resistance coefficients of the one-dimensional wellbore line element and the zero-dimensional well bottom point element respectively; construct the corresponding well flow conductivity matrix based on the initial values ​​of the resistance coefficients of the two types of elements respectively; then merge the two types of well flow conductivity matrices into the total stiffness matrix obtained in S2, and solve the linear equation system to obtain the initial well flow rate and the nodal head of the entire field; S4. Deploy sampling points around the well and interpolate to obtain the sampling head: Based on the well perimeter unit and the spatial position of the well unit obtained in S2, deploy multiple layers of annular sampling points along the well axis for the one-dimensional wellbore line unit, and deploy annular sampling points on the bottom plane for the bottom point unit; Based on the initial global nodal head obtained in S3 or the global nodal head obtained in the previous iteration, obtain the head value of each sampling point by linear interpolation of the nodal head within the set range of the sampling points, and calculate the average head of each sampling ring. S5. Update the iterative resistance coefficients and recalculate the seepage field: Based on the average head of each sampling ring obtained in S4, calculate the current iterative flow rate and head difference of the one-dimensional wellbore line element and the zero-dimensional bottom point element respectively using the radial flow well flow analytical formula; update the well flow resistance coefficients of the two types of elements using the relaxation average iterative scheme, and construct the well flow conductivity matrix of the current iteration step; incorporate the well flow conductivity matrix of the current iteration step into the total stiffness matrix obtained in S2, and solve for the well flow rate and the head of all nodes in the field of the current iteration step. S6. Perform iterative convergence judgment and output the final result: Calculate the relative residual of the well flow rate in the current iteration step; if the relative residual is less than the iterative convergence threshold set in S1, the iteration terminates and the final well flow rate and the total head of all nodes are output; if the relative residual is not less than the iterative convergence threshold, return to execute S4 to S5 to enter the next round of iterative calculation.

[0006] Preferably, step S2 specifically includes the following steps: S21. Read all the basic parameters determined in step S1 as the input basis for mesh generation, well element construction and stiffness matrix assembly; S22. Based on the model boundary and soil layer spatial parameters loaded in S21, the geometric model of the three-dimensional seepage computation domain is divided into Delaunay meshes; and the characteristic dimensions of the well perimeter elements are extracted. ; S23. Based on the geometric parameters of the decompression well loaded in S21, the three-dimensional decompression well is discretized into multi-dimensional units. S231. Discretize the three-dimensional wellbore segment along the wellbore axis into several one-dimensional wellbore linear elements. The length of a single wellbore linear element is denoted as... Each line segment unit contains two nodes: the beginning and the end. S232. Discretize the three-dimensional bottom section into zero-dimensional bottom point elements, corresponding to a single node at the bottom of the well; S233. For incomplete wells, the bottom point unit is configured to simultaneously bear the dual role of seepage resistance in the wellbore section and seepage resistance in the bottom section; for complete wells, the bottom point unit is configured to only bear the role of seepage resistance in the bottom section. S24. The one-dimensional wellbore line element and zero-dimensional well bottom point element obtained in S23 are embedded into the corresponding spatial position of the three-dimensional finite element calculation mesh of soil generated in S22 according to the spatial coordinates of the decompression well loaded in S21, so as to construct a mixed-dimensional calculation mesh system of 3D soil-1D wellbore-0D well bottom. S25. Based on the geometric properties of the Delaunay mesh elements obtained in S22, combined with the corresponding soil permeability coefficient loaded in S21. A three-dimensional total stiffness matrix of seepage is assembled based on the finite element control equations, which serves as the basis matrix for subsequent well-flow coupling iterative solutions.

[0007] Preferably, step S3 specifically includes the following steps: S31. Extract the soil permeability coefficient from step S1. Radius of relief well Extract the feature size of the well perimeter unit from step S2. Length of each one-dimensional wellbore line unit The feature angle is set based on empirical values. ; S32. Based on the parameters extracted in S31, calculate the initial values ​​of the resistance coefficients for the one-dimensional wellbore line element and the zero-dimensional bottom point element. For each one-dimensional wellbore line element, the initial value of the drag coefficient of the one-dimensional wellbore element in the first iteration is calculated using the axisymmetric radial flow analytical formula. : ; For zero-dimensional bottom-hole element, the initial value of the drag coefficient of the zero-dimensional bottom-hole element in the first iteration is calculated using the analytical formula for the hemispherical flow convergence at the bottom of the well. : ; S33. Based on the initial values ​​of the resistance coefficient of the one-dimensional wellbore unit and the initial values ​​of the resistance coefficient of the zero-dimensional well bottom unit obtained in S32, construct the first-round well flow water conduction matrix of the corresponding unit. Initial value of the resistance coefficient of the one-dimensional wellbore element For each one-dimensional wellbore line element, a fourth-order symmetric well flow conductivity matrix is ​​constructed by introducing a double-sided coupling form between the wellbore nodes and the additional seepage path into the element nodes: ; in, ; ; In the formula, Let be the well flow conductivity of a one-dimensional wellbore linear element node, during the first iteration. = ; The flow conductivity of the element nodes is introduced to account for the additional seepage diameter of the one-dimensional wellbore element, where... ; and These represent the water heads at both ends of a one-dimensional wellbore linear element; and The additional seepage diameters at both ends of the one-dimensional wellbore line element introduce water head into the element node; and These represent the seepage flow rates at the upper and lower nodes of a one-dimensional wellbore unit, respectively. and These represent the seepage flow rates introduced into the unit nodes by the additional seepage paths at the upper and lower ends of the one-dimensional wellbore line unit, respectively. For the first The wellbore node corresponds to the sampling circle on the [number]th [number]th [number]. The head value at each sampling point; For the first The water head at each wellhead node; The radius of the well perimeter sampling zone; The total number of sampling points is 12; For the first The wellbore node corresponds to the sampling circle on the [number]th [number]th [number]. The head value at each sampling point; for The water head at each wellhead node; Initial value of the resistance coefficient of the zero-dimensional bottom hole element For zero-dimensional bottom-point elements, a second-order symmetric well flow conductivity matrix is ​​constructed that couples the bottom-point nodes with the nodes of the additional seepage path introduction elements: ; in, ; In the formula, The well flow conductivity is the coefficient of the zero-dimensional bottom-point element during the first iteration. ; The water head at the bottom node of the well; Introducing head into the unit node for additional seepage path; and These are the seepage flow rates corresponding to the bottom node and the seepage flow rates corresponding to the nodes of the additional seepage path introduction unit, respectively. The water head at the bottom node of the well on the side of the soil; S34. The one-dimensional wellbore line unit flow water conduction matrix and the zero-dimensional well bottom point unit flow water conduction matrix obtained in S33 are superimposed and incorporated into the three-dimensional seepage total stiffness matrix assembled in step S2 according to the correspondence of node spatial coordinates to obtain the first round of total stiffness matrix. S35. Initial seepage field solution: Based on the first-round total stiffness matrix containing well flow coupling terms obtained in S34, and combined with the model boundary conditions and decompression well head boundary conditions determined in step S1, solve the three-dimensional stable seepage linear equation system to obtain the full-field nodal head distribution and the initial well flow rate of each well unit corresponding to the first round of iteration.

[0008] Preferably, during the superposition process in step S34, the diagonal elements of the matrix corresponding to the well node are superimposed with the corresponding well flow conductivity coefficient, the intersection of the well node and the node of the additional seepage path introduced by the unit node is filled with the corresponding negative conductivity coefficient, and the remaining matrix elements keep the original three-dimensional seepage total stiffness matrix unchanged, thus assembling to form the first round total stiffness matrix containing the well flow coupling term.

[0009] Preferably, step S4 specifically includes the following steps: S41. Sampling and Layout Basic Parameter Extraction: Extract the spatial coordinates and element length of each one-dimensional wellbore line element, the spatial coordinates of the zero-dimensional well bottom point element, and the characteristic dimensions of the well perimeter element from step S2. The first iteration uses the initial total head obtained from step S3, while subsequent iterations use the total head obtained from the previous iteration. S42. Layout of the annular sampling circle: Based on the spatial coordinates and element size parameters of the one-dimensional wellbore line element extracted by S41, for each one-dimensional wellbore line element segment, annular sampling rings are arranged at three positions along the wellbore axis: the upper end, the midpoint, and the lower end of the element. The center of each sampling ring coincides with the wellbore axis, and the radius of the sampling ring is... Values ​​and characteristic dimensions of well perimeter units Equal; each sampling circle is evenly distributed along the circumference. One sampling point; Simultaneously, based on the spatial coordinates of the zero-dimensional well bottom point unit extracted by S41, a ring of sampling points is arranged on the horizontal plane where the well bottom is located for the zero-dimensional well bottom point unit; the center of the sampling ring coincides with the axis of the well body, and the radius of the sampling ring is... Values ​​and characteristic dimensions of well perimeter units Equal; sampling circles are evenly distributed along the circumference. One sampling point; S43. Calculation of head interpolation for each sampling point: Based on the head data of all nodes in the field determined in S41, and combined with the spatial coordinates of all sampling points obtained in S42, the head value corresponding to each sampling point is obtained by linear interpolation of the head of the nodes within the set range of the sampling points. S44. Calculation of Average Head for Each Sampling Loop: Based on the head values ​​of all sampling points obtained in S43, calculate the average head for each wellbore sampling loop and the bottom sampling loop; for any sampling loop in the wellbore section, its average head... The calculation formula is: ; In the formula, For the sampling circle The water head value at each sampling point.

[0010] Preferably, step S5 specifically includes the following steps: S51. Extract the average water head of each sampling loop from step S4, and extract the soil permeability coefficient from step S1. Radius of relief well , well water head Extract the length of each one-dimensional wellbore line unit from step S2. well perimeter sampling radius Extract the one-dimensional wellbore line element resistance coefficient from the previous iteration results. Zero-dimensional well bottom point element resistance coefficient The first iteration after the first iteration uses the initial value of the resistance coefficient obtained from S3.

[0011] S52. Based on the parameters extracted in S51, and combined with the axisymmetric radial flow analytical formula, calculate the comprehensive flow-head difference ratio of a single well section unit. : ; In the formula, This represents the equivalent flow rate of the wellbore unit in the current iteration step. This is the difference between the water head at the node of the well unit and the constant water head inside the well; The analytical flow rate at the midpoint of a one-dimensional wellbore element; , , These are the nodal heads at the top, bottom, and midpoint of a one-dimensional well unit, respectively. Simultaneously, by combining the analytical formula for the hemispherical flow at the bottom of the well, the ratio of flow rate to head difference in the bottom unit of the current iteration step is calculated. : ; In the formula, This represents the equivalent flow rate of the bottom cell in the current iteration step. This is the difference between the water head at the bottom unit node and the constant water head inside the well; S53. Iterative resistance coefficient update: The current iteration step resistance coefficients of the one-dimensional wellbore line element and the zero-dimensional well bottom point element are updated using the relaxed average iterative format. The update expression for the resistance coefficient of the one-dimensional wellbore element is as follows: ; In the formula, For the first The resistance coefficient of the one-dimensional wellbore line element in the next iteration; For the first The resistance coefficient of the one-dimensional wellbore line element in the next iteration; The update expression for the drag coefficient of the zero-dimensional bottom-hole element is as follows: ; In the formula, For the first The drag coefficient of the zero-dimensional bottom-point element in the next iteration; For the first The drag coefficient of the zero-dimensional bottom-point element in the next iteration; S54. Based on the resistance coefficient obtained in S53 for the current iteration step, construct the well flow conductivity matrix for one-dimensional wellbore line elements and zero-dimensional well bottom point elements, respectively. The matrix structure is consistent with step S33, only the conductivity coefficient values ​​are replaced with the calculation results of the current iteration step; where the conductivity coefficient of the wellbore element... Bottom unit water conductivity ; S55. The water conduction matrices of the one-dimensional wellbore line element and the zero-dimensional well bottom point element constructed in S54 are superimposed and incorporated into the three-dimensional total seepage stiffness matrix obtained in step S2 according to the correspondence of node spatial coordinates and the superposition rules in step S34: ; In the formula, For the first The node and the first The seepage stiffness coefficient between nodes, and ; For a one-dimensional wellbore line element, the in-well pipe flow resistance coefficient is given. and All values ​​are water head values ​​for ordinary soil nodes within the aquifer calculation domain; and These represent the water heads at the upper and lower well nodes of a one-dimensional well unit; and The additional seepage paths at the top and bottom of the one-dimensional well unit introduce water head into the unit node; and These represent the seepage flow rates at ordinary nodes in the soil region; S56. Based on the total stiffness matrix of the current iteration step containing the well flow coupling term obtained in S55, and combined with the model boundary conditions and decompression well head boundary conditions determined in step S1, solve the three-dimensional stable seepage linear equation system to obtain the full field node head distribution and well flow rate of each well unit corresponding to the current iteration step.

[0012] Preferably, the relative residual mentioned in step S6 is the ratio of the residual flow rate to the boundary total flow rate.

[0013] Therefore, the present invention employs the above-mentioned method for simulating well flow in a three-dimensional finite element unstructured mesh, which has the following beneficial effects: 1. Significantly improved computational efficiency and greatly reduced resource input: By reducing the dimensionality of the three-dimensional solid decompression well to one-dimensional line elements and zero-dimensional point elements, the total number of model elements is effectively reduced. There is no need to perform extreme densification of the well perimeter mesh, which can significantly shorten the computation time and reduce the consumption of computing storage and computing power resources. 2. Stable calculation accuracy with weak correlation to grid size: The calculation error does not deteriorate sharply with the increase of the grid size around the well. Even if the grid size reaches hundreds of times the well radius, the flow calculation error of complete wells and incomplete wells can still be kept within a controllable range, which solves the problem of significant decrease in accuracy under coarse grids in traditional methods. 3. Strong adaptability to mesh and working conditions: It natively supports 3D Delaunay unstructured meshes and is applicable to various working conditions such as complete wells, incomplete wells, single wells and multi-well groups. It can effectively simulate inter-well interference effects and has a wide range of engineering applications. 4. The method is robust and easy to implement in engineering: The initial value of the resistance coefficient is calculated using analytical formulas, and the iterative update format of relaxation averaging is used, resulting in a stable iterative convergence process; the well flow conductivity matrix can be directly superimposed on the original three-dimensional seepage total stiffness matrix without the need to reconstruct the finite element solution framework, making it easy to embed into existing seepage calculation programs.

[0014] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0015] Figure 1 This is a flowchart of a well flow simulation method in a three-dimensional finite element unstructured mesh according to the present invention; Figure 2This is a schematic diagram of three-dimensional wellbore dimensionality reduction for a three-dimensional finite element unstructured mesh well flow simulation method according to the present invention; Figure 3 This is a schematic diagram of the sampling points around the well in a three-dimensional finite element unstructured mesh well flow simulation method according to the present invention, wherein (a) is a schematic diagram of the sampling point layout of the zero-dimensional bottom point element, and (b) is a schematic diagram of the sampling point layout of the one-dimensional wellbore line element. Figure 4 This is a schematic diagram of a well flow simulation method in a three-dimensional finite element unstructured mesh according to the present invention, wherein (a) is a schematic diagram of a zero-dimensional bottom point element structure, and (b) is a schematic diagram of a one-dimensional wellbore line element structure. Figure 5 This is a schematic diagram of the well point distribution in the multi-well example described in Example 2; Figure 6 This is a comparison chart of flow rate errors in the single-well complete well calculation example described in Example 1; Figure 7 This is a comparison chart of flow rate errors in the single-well incomplete well calculation example described in Example 1; Figure 8 This is a comparison chart of flow rate errors in the multi-well complete well calculation example described in Example 2; Figure 9 This is a comparison chart of flow rate errors in the multi-well incomplete well calculation example described in Example 2; Figure 10 This is a comparison chart of flow rate errors in the well spacing sensitivity analysis described in Example 3; Figure 11 This is a comparison chart of flow error in the boundary distance sensitivity analysis described in Example 4. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the embodiments of the present invention and are not intended to limit the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of this application. Examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout.

[0017] It should be noted that the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion, such as a process, method, system, product, or server that includes a series of steps or units, not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such process, method, product, or device.

[0018] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings. Specific Implementation This embodiment comprehensively verifies, through four sets of numerical examples, the following four dimensions: single-well accuracy, multi-well applicability, well spacing sensitivity, and boundary sensitivity: Figures 1-4 The computational accuracy, computational efficiency, and applicable boundaries of the well flow simulation method in a three-dimensional finite element unstructured mesh described in this invention are illustrated. All examples use traditional equal-dimensional finite element solutions or theoretical analytical solutions under a fine mesh with a size of 0.1 times the well radius as reference benchmarks to evaluate the performance of this method under different coarse mesh conditions.

[0020] Example 1: Single Well Example This example uses a cylindrical aquifer calculation model, with the soil layer being a single homogeneous layer and the soil permeability coefficient being... The aquifer thickness is 10m; the pressure relief well is located at the center of the cylinder, with a well radius of... The wellhead was set to 0m; the model's outer boundary radius was 100m, and a constant water head of 1m was applied to the outer boundary. During the calculation, various Delaunay tetrahedral unstructured meshes with different element sizes were used to discretize the aquifer computational domain. Simulations were conducted for two conditions: a complete well and a non-complete well with a penetration depth of 1 / 2. The results are as follows: Figure 6 and Figure 7 The results show the comparison of flow error.

[0021] Calculation Example 2: Multi-well Calculation Example This example uses a cylindrical aquifer calculation model, with the soil layer being a single homogeneous layer and the soil permeability coefficient being... The aquifer thickness is 40m; five pressure relief wells are set up in the model, each with a radius of [missing information]. The water head inside the well is 0m; a water head of 1m is applied to the outer boundary of the model. The specific arrangement of the well points is as follows: Figure 5 As shown. During the calculation, various Delaunay tetrahedral unstructured meshes with different element sizes were used to discretize the aquifer computational domain. Simulations were conducted for two conditions: a complete well and a non-complete well with a penetration depth of 1 / 2, to verify the applicability of this method under seepage interference conditions in a group of wells. The results are as follows. Figure 8 and Figure 9 The results show the comparison of flow error.

[0022] Example 3: Well Spacing Sensitivity Analysis Example This example is used to determine the applicable well spacing range for the flow error convergence of the method of the present invention. A dual-well cylindrical aquifer model is used, and the soil layer is a single homogeneous soil layer with a permeability coefficient of [missing information]. The aquifer thickness is 10m; the radius of the pressure relief well is... The water head inside the well is set to 0m; a constant water head of 1m is applied to the outer boundary of the model. Calculations are based on the well spacing. With well radius The ratio is a variable, and multiple sets of comparison calculations are carried out using incomplete wells with a penetration depth of 1 / 2. The results are as follows: Figure 10 The results show the comparison of flow error.

[0023] Example 4: Boundary Distance Sensitivity Analysis Example This example is used to determine the influence range of the model boundary distance on the calculation accuracy of this method. A cuboid aquifer calculation model with dimensions of 100m × 100m × 10m is used. The soil layer is a single homogeneous soil layer with a permeability coefficient of [missing value]. The aquifer thickness is 1m; the radius of the pressure relief well is... The water head inside the well is set to 0m; a water head of 1m is applied to the outer boundary of the model. Calculate the distance from the center of the well to the boundary of the model. With well radius The ratio is a variable, and multiple sets of comparison calculations are carried out using incomplete wells with a penetration depth of 1 / 2. The results are as follows: Figure 11 The results show the comparison of flow error.

[0024] Example Results and Analysis Single-well accuracy performance: For complete well conditions, under coarse grid conditions with a grid size increased to 100 times the well radius (unit feature size 10m), the relative error between the calculated well flow rate and the analytical solution can be stabilized within ±5%. Under the same grid conditions, the calculation accuracy is significantly better than simulation methods that only perform geometric dimensionality reduction without resistance correction. For incomplete well conditions with a penetration depth of 1 / 2, this method also has good adaptability. The flow rate error is stable within ±5% at different grid scales, and the error does not increase significantly with the increase of grid size, fully demonstrating the weak correlation between the calculation accuracy of this method and grid size.

[0025] Multi-well adaptability: In a well group system containing five wells, whether in complete well or incomplete well conditions, the relative error between the total flow rate of each well calculated by this method and the reference value does not exceed ±6%, and the flow distribution ratio between each well is consistent with the reference result. This indicates that this method can effectively simulate the seepage interference effect between well groups and can be applied to engineering scenarios with complex multi-well layouts.

[0026] Sensitivity Analysis and Applicable Boundaries: The well spacing sensitivity analysis results show that when the distance between two wells is greater than 5 times the well radius, the flow rate calculation error of this method converges and stabilizes within ±5%. Boundary distance sensitivity analysis reveals a consistent pattern: when the distance from the well center to the model boundary is greater than 5 times the well radius, the boundary effect on calculation accuracy is negligible, and the flow rate error remains within ±5%.

[0027] Improved computational efficiency: Since this method adopts a geometric dimensionality reduction approach, the three-dimensional wellbore is reduced to one-dimensional line elements and zero-dimensional point elements. There is no need to perform extreme densification of the wellbore mesh. The total number of model elements is significantly reduced compared with traditional fine mesh equal-dimensional simulation methods, the overall calculation time is significantly shortened, and the improvement in computational efficiency is obvious.

[0028] The results of the four examples show that the core advantage of the method of this invention lies in the weak correlation between its accuracy and the mesh size. Its calculation error fluctuates very little within a huge range of mesh sizes, from 10 to 100 times the well radius, while the error of traditional methods deteriorates sharply as the mesh size becomes coarser. In terms of computational efficiency, by reducing the three-dimensional wellbore to one-dimensional line elements and zero-dimensional point elements, the number of model elements is reduced, and the corresponding computation time is shortened compared to traditional fine-mesh methods.

[0029] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for simulating well flow in a three-dimensional finite element unstructured mesh, characterized in that: Includes the following steps: S1. Determine the basic parameters for model calculation: Determine the boundary conditions of the geometric model of the three-dimensional seepage calculation domain, the water head boundary conditions of the pressure relief well, the spatial coordinates of the pressure relief well, the geometric parameters of the pressure relief well, the soil permeability coefficient, the soil spatial parameters, and the iterative convergence threshold parameters. S2. Constructing a three-dimensional unstructured mesh seepage model and assembling the overall stiffness matrix: Based on all the basic parameters determined in S1, the geometric model of the three-dimensional seepage computation domain is divided into Delaunay meshes. At the same time, the three-dimensional decompression well is geometrically reduced. The well section is discretized into one-dimensional well line elements along the well axis, and the bottom section is discretized into zero-dimensional bottom point elements. The one-dimensional well line elements and zero-dimensional bottom point elements are embedded into the corresponding spatial coordinate positions of the three-dimensional finite element computation mesh of the soil to construct a 3D-1D-0D mixed-dimensional computation mesh. Finally, based on the geometric properties of the three-dimensional finite element computation mesh of the soil and the permeability coefficient parameters of the corresponding soil layer, the three-dimensional seepage overall stiffness matrix is ​​assembled. S3. Calculate the first-round well flow resistance coefficient and solve the initial seepage field: Based on the characteristic dimensions of the well perimeter element obtained in S2, and combined with the soil permeability coefficient and decompression well radius parameter in S1, calculate the initial values ​​of the resistance coefficients of the one-dimensional wellbore line element and the zero-dimensional well bottom point element respectively; construct the corresponding well flow conductivity matrix based on the initial values ​​of the resistance coefficients of the two types of elements respectively. Then, the two types of well flow conductivity matrices are incorporated into the total stiffness matrix obtained by S2, and the initial well flow rate and the total field nodal head are obtained by solving the linear equations. S4. Deploy sampling points around the well and interpolate to obtain the sampling head: Based on the well perimeter unit and the spatial position of the well unit obtained in S2, deploy multiple layers of annular sampling points along the well axis for the one-dimensional wellbore line unit, and deploy annular sampling points on the bottom plane for the bottom point unit; Based on the initial global nodal head obtained in S3 or the global nodal head obtained in the previous iteration, obtain the head value of each sampling point by linear interpolation of the nodal head within the set range of the sampling points, and calculate the average head of each sampling ring. S5. Update the iterative resistance coefficient and recalculate the seepage field: Based on the average head of each sampling ring obtained in S4, calculate the current iterative flow rate and head difference of the one-dimensional wellbore line element and the zero-dimensional bottom point element respectively by combining the radial flow well flow analytical formula. The well flow resistance coefficients of the two types of cells are updated using a relaxed average iterative scheme to construct the well flow conductivity matrix for the current iteration step; The well flow conductivity matrix of the current iteration step is incorporated into the total stiffness matrix obtained by S2, and the well flow rate and the total field head of the nodes are solved to obtain the current iteration step. S6. Perform iterative convergence determination and output the final result: calculate the relative residual of the well flow rate in the current iteration step; If the relative residual is less than the iteration convergence threshold set by S1, the iteration terminates and the final well flow rate and total node head calculation results are output. If the relative residual is not less than the iteration convergence threshold, then return to execute S4 to S5 to enter the next round of iteration calculation.

2. The well flow simulation method in a three-dimensional finite element unstructured mesh according to claim 1, characterized in that: Step S2 specifically includes the following steps: S21. Read all the basic parameters determined in step S1 as the input basis for mesh generation, well element construction and stiffness matrix assembly; S22. Based on the model boundary and soil layer spatial parameters loaded in S21, the geometric model of the three-dimensional seepage computation domain is divided into Delaunay meshes; and the characteristic dimensions of the well perimeter elements are extracted. ; S23. Based on the geometric parameters of the decompression well loaded in S21, the three-dimensional decompression well is discretized into multi-dimensional units. S231. Discretize the three-dimensional wellbore segment along the wellbore axis into several one-dimensional wellbore linear elements. The length of a single wellbore linear element is denoted as... Each line segment unit contains two nodes: the beginning and the end. S232. Discretize the three-dimensional bottom section into zero-dimensional bottom point elements, corresponding to a single node at the bottom of the well; S233. For incomplete wells, the bottom point unit is configured to simultaneously bear the dual role of seepage resistance in the wellbore section and seepage resistance in the bottom section; for complete wells, the bottom point unit is configured to only bear the role of seepage resistance in the bottom section. S24. The one-dimensional wellbore line element and zero-dimensional well bottom point element obtained in S23 are embedded into the corresponding spatial position of the three-dimensional finite element calculation mesh of soil generated in S22 according to the spatial coordinates of the decompression well loaded in S21, so as to construct a mixed-dimensional calculation mesh system of 3D soil-1D wellbore-0D well bottom. S25. Based on the geometric properties of the Delaunay mesh elements obtained in S22, combined with the corresponding soil permeability coefficient loaded in S21. A three-dimensional total stiffness matrix of seepage is assembled based on the finite element control equations, which serves as the basis matrix for subsequent well-flow coupling iterative solutions.

3. A method for simulating well flow in a three-dimensional finite element unstructured mesh according to claim 1 or 2, characterized in that: Step S3 specifically includes the following steps: S31. Extract the soil permeability coefficient from step S1. Radius of relief well Extract the feature size of the well perimeter unit from step S2. Length of each one-dimensional wellbore line unit The feature angle is set based on empirical values. ; S32. Based on the parameters extracted in S31, calculate the initial values ​​of the resistance coefficients for the one-dimensional wellbore line element and the zero-dimensional bottom point element. For each one-dimensional wellbore line element, the initial value of the drag coefficient of the one-dimensional wellbore element in the first iteration is calculated using the axisymmetric radial flow analytical formula. : ； For zero-dimensional bottom-hole element, the initial value of the drag coefficient of the zero-dimensional bottom-hole element in the first iteration is calculated using the analytical formula for the hemispherical flow convergence at the bottom of the well. : ； S33. Based on the initial values ​​of the resistance coefficient of the one-dimensional wellbore unit and the initial values ​​of the resistance coefficient of the zero-dimensional well bottom unit obtained in S32, construct the first-round well flow water conduction matrix of the corresponding unit. Based on the initial value of the resistance coefficient of the one-dimensional wellbore element For each one-dimensional wellbore line element, a fourth-order symmetric well flow conductivity matrix is ​​constructed by introducing a double-sided coupling form between the wellbore nodes and the additional seepage path into the element nodes: ； in, ； ； In the formula, Let be the well flow conductivity of a one-dimensional wellbore linear element node, during the first iteration. = ; The flow conductivity of the wellbore is introduced into the element nodes to represent the additional seepage diameter on the wellbore side of a one-dimensional wellbore line element, where... ; and These represent the water heads at both ends of a one-dimensional wellbore linear element; and The additional seepage diameters at both ends of the one-dimensional wellbore line element introduce water head into the element node; and These represent the seepage flow rates at the upper and lower nodes of a one-dimensional wellbore unit, respectively. and These represent the seepage flow rates introduced into the unit nodes by the additional seepage paths at the upper and lower ends of the one-dimensional wellbore line unit, respectively. For the first The wellbore node corresponds to the sampling circle on the [number]th [number]th [number]. The head value at each sampling point; For the first The water head at each wellhead node; The radius of the well perimeter sampling zone; This represents the total number of sampling points; For the first The wellbore node corresponds to the sampling circle on the [number]th [number]th [number]. The head value at each sampling point; for The water head at each wellhead node; Initial value of the resistance coefficient of the zero-dimensional bottom hole element For zero-dimensional bottom-point elements, a second-order symmetric well flow conductivity matrix is ​​constructed that couples the bottom-point nodes with the nodes of the additional seepage path introduction elements: ； in, ； In the formula, The well flow conductivity is the coefficient of the zero-dimensional bottom-point element during the first iteration. ; The water head at the bottom node of the well; Introducing head into the unit node for additional seepage path; and These are the seepage flow rates corresponding to the bottom node and the seepage flow rates corresponding to the nodes of the additional seepage path introduction unit, respectively. The water head at the bottom node of the well on the side of the soil; S34. The one-dimensional wellbore line unit flow water conduction matrix and the zero-dimensional well bottom point unit flow water conduction matrix obtained in S33 are superimposed and incorporated into the three-dimensional seepage total stiffness matrix assembled in step S2 according to the correspondence of node spatial coordinates to obtain the first round of total stiffness matrix. S35. Initial seepage field solution: Based on the first-round total stiffness matrix containing well flow coupling terms obtained in S34, and combined with the model boundary conditions and decompression well head boundary conditions determined in step S1, solve the three-dimensional stable seepage linear equation system to obtain the full-field nodal head distribution and the initial well flow rate of each well unit corresponding to the first round of iteration.

4. The well flow simulation method in a three-dimensional finite element unstructured mesh according to claim 3, characterized in that: During the superposition process in step S34, the diagonal elements of the matrix corresponding to the well node are superimposed with the corresponding well flow conductivity coefficient. The intersection of the well node and the node of the additional seepage path introduced by the unit is filled with the corresponding negative conductivity coefficient. The remaining matrix elements keep the original three-dimensional seepage total stiffness matrix unchanged, and the first round total stiffness matrix containing the well flow coupling term is assembled.

5. The well flow simulation method in a three-dimensional finite element unstructured mesh according to claim 4, characterized in that: Step S4 specifically includes the following steps: S41. Sampling and Layout Basic Parameter Extraction: Extract the spatial coordinates and element length of each one-dimensional wellbore line element, the spatial coordinates of the zero-dimensional well bottom point element, and the characteristic dimensions of the well perimeter element from step S2. The first iteration uses the initial total head obtained from step S3, while subsequent iterations use the total head obtained from the previous iteration. S42. Layout of the annular sampling circle: Based on the spatial coordinates and element size parameters of the one-dimensional wellbore line element extracted by S41, for each one-dimensional wellbore line element segment, annular sampling rings are arranged at three positions along the wellbore axis: the upper end, the midpoint, and the lower end of the element. The center of each sampling ring coincides with the wellbore axis, and the radius of the sampling ring is... Values ​​and characteristic dimensions of well perimeter units Equal; each sampling circle is evenly distributed along the circumference. One sampling point; Simultaneously, based on the spatial coordinates of the zero-dimensional well bottom point unit extracted by S41, a ring of sampling points is arranged on the horizontal plane where the well bottom is located for the zero-dimensional well bottom point unit; the center of the sampling ring coincides with the axis of the well body, and the radius of the sampling ring is... Values ​​and characteristic dimensions of well perimeter units Equal; sampling circles are evenly distributed along the circumference. One sampling point; S43. Calculation of head interpolation for each sampling point: Based on the head data of all nodes in the field determined in S41, and combined with the spatial coordinates of all sampling points obtained in S42, the head value corresponding to each sampling point is obtained by linear interpolation of the head of the nodes within the set range of the sampling points. S44. Calculation of Average Head for Each Sampling Loop: Based on the head values ​​of all sampling points obtained in S43, calculate the average head for each wellbore sampling loop and the bottom sampling loop; for any sampling loop in the wellbore section, its average head... The calculation formula is: ； In the formula, For the sampling circle The water head value at each sampling point.

6. The well flow simulation method in a three-dimensional finite element unstructured mesh according to claim 5, characterized in that: Step S5 specifically includes the following steps: S51. Extract the average water head of each sampling loop from step S4, and extract the soil permeability coefficient from step S1. Radius of relief well , well water head Extract the length of each one-dimensional wellbore line unit from step S2. well perimeter sampling radius Extract the one-dimensional wellbore line element resistance coefficient from the previous iteration results. Zero-dimensional well bottom point element resistance coefficient The first iteration after the first iteration uses the initial value of the resistance coefficient obtained from S3. S52. Based on the parameters extracted in S51, and combined with the axisymmetric radial flow analytical formula, calculate the comprehensive flow-head difference ratio of a single well section unit. : ； In the formula, This represents the equivalent flow rate of the wellbore unit in the current iteration step. This is the difference between the water head at the node of the well unit and the constant water head inside the well; The analytical flow rate at the midpoint of a one-dimensional wellbore element; , , These are the nodal heads at the top, bottom, and midpoint of a one-dimensional well unit, respectively. Simultaneously, by combining the analytical formula for the hemispherical flow at the bottom of the well, the ratio of flow rate to head difference in the bottom unit of the current iteration step is calculated. : ； In the formula, This represents the equivalent flow rate of the bottom cell in the current iteration step. This is the difference between the water head at the bottom unit node and the constant water head inside the well; S53. Iterative resistance coefficient update: The current iteration step resistance coefficients of the one-dimensional wellbore line element and the zero-dimensional well bottom point element are updated using the relaxed average iterative format. The update expression for the resistance coefficient of the one-dimensional wellbore element is as follows: ； In the formula, For the first The resistance coefficient of the one-dimensional wellbore line element in the next iteration; For the first The resistance coefficient of the one-dimensional wellbore line element in the next iteration; The update expression for the drag coefficient of the zero-dimensional bottom-hole element is as follows: ； In the formula, For the first The drag coefficient of the zero-dimensional bottom-point element in the next iteration; For the first The drag coefficient of the zero-dimensional bottom-point element in the next iteration; S54. Based on the resistance coefficient obtained in S53 for the current iteration step, construct the well flow conductivity matrix for one-dimensional wellbore line elements and zero-dimensional well bottom point elements, respectively. The matrix structure is consistent with step S33, only the conductivity coefficient values ​​are replaced with the calculation results of the current iteration step; where the conductivity coefficient of the wellbore element... Bottom unit water conductivity ; S55. The water conduction matrices of the one-dimensional wellbore line element and the zero-dimensional well bottom point element constructed in S54 are superimposed and incorporated into the three-dimensional total seepage stiffness matrix obtained in step S2 according to the correspondence of node spatial coordinates and the superposition rules in step S34: ； In the formula, For the first The node and the first The seepage stiffness coefficient between nodes, and ; For a one-dimensional wellbore line element, the in-well pipe flow resistance coefficient is given. and All values ​​are water head values ​​for ordinary soil nodes within the aquifer calculation domain; and These represent the water heads at the upper and lower well nodes of a one-dimensional well unit; and The additional seepage paths at the top and bottom of the one-dimensional well unit introduce water head into the unit node; and These represent the seepage flow rates at ordinary nodes in the soil region; S56. Based on the total stiffness matrix of the current iteration step containing the well flow coupling term obtained in S55, and combined with the model boundary conditions and decompression well head boundary conditions determined in step S1, solve the three-dimensional stable seepage linear equation system to obtain the full field node head distribution and well flow rate of each well unit corresponding to the current iteration step.

7. The well flow simulation method in a three-dimensional finite element unstructured mesh according to claim 6, characterized in that: The relative residual mentioned in step S6 is the ratio of the residual flow rate to the total boundary flow rate.

Images