Finite Element Mesh Adaptive Subdivision Method for Solving Convection-Based Pollutant Migration

Abstract

The present invention provides a finite element grid adaptive subdivision method for solving convection-based pollutant migration, which relates to the field of finite element numerical simulation technology, including: identifying convection-based finite element grids through a grid adaptive subdivision algorithm, and subdividing the convection-based finite element grids to obtain subdivided subunits; for each unsubdivided unit shared by an active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits; for the units where all the remaining active edges are located, the remaining active edges are connected to the diagonal vertices of the remaining active edges to form 2 subdivided subunits. After the subdivided subunits are divided, the subdivided subunits have the same unit number as the parent unit, so that the subdivided subunits can reference the parameters of the parent unit. The present invention can effectively solve the problems of vibration of numerical solutions and large numerical dispersion in traditional numerical models when solving convection-based solute migration problems.

Description

Technical Field

[0001] The invention relates to the technical field of finite element numerical simulation, and in particular to a finite element grid adaptive subdivision method for solving convection-based pollutant migration. Background Art

[0002] In the finite element numerical method, for the pollutant migration problem dominated by diffusion (the judgment standard is Peclet number < 2), the finite element differential equation has the characteristics of a self-harmonic differential equation, and the standard Galerkin method is optimal at this time. The typical groundwater numerical simulation software FEFLOW mostly uses the standard Galerkin method, so when solving the pollutant migration problem dominated by diffusion, the numerical solution has good stability.

[0003] However, when solving pollutant migration problems dominated by convection (the judgment standard is Peclet number ≥ 2), the standard Galerkin method will cause vibration of the numerical solution. The usual solution is to use the upstream weighted method instead. This algorithm can effectively overcome the vibration problem of the numerical solution, but at the expense of the accuracy of the solution (i.e., large numerical dispersion). In general, when solving pollutant migration problems dominated by convection (the judgment standard is Peclet number ≥ 2) in traditional groundwater models, the standard Galerkin method and upstream weighted method used will bring about problems such as vibration of the numerical solution and large numerical dispersion. Summary of the invention

[0004] In order to overcome the deficiencies of the prior art, the object of the present invention is to provide a finite element mesh adaptive subdivision method for solving convection-based pollutant migration.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A finite element mesh adaptive subdivision method for solving convection-based pollutant migration, comprising:

[0007] The finite element mesh dominated by convection is identified by a mesh adaptive subdivision algorithm, and the finite element mesh dominated by convection is subdivided to obtain subdivided sub-units;

[0008] For each unsubdivided unit shared by the active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits;

[0009] For the cells where all the remaining active edges are located, connect the remaining active edges and the diagonal vertices of the remaining active edges to form two subdivided sub-cells;

[0010] After the subdivision unit is divided, the subdivision unit has the same unit number as the parent unit, so that the subdivision unit can reference the parameters of the parent unit.

[0011] Preferably, the identification method of the finite element mesh based on convection includes:

[0012] It is determined whether the Peclet number is greater than or equal to 2. If so, the mesh is determined to be the convection-dominated finite element mesh.

[0013] Preferably, the process of subdividing the finite element mesh mainly based on convection further includes:

[0014] Check the properties of the three edges of the finite element mesh dominated by convection, assuming that the three edges of the finite element mesh dominated by convection are ij = 1, jk = 2, ki = 3. Consider the ij edge, if the ij edge is a boundary edge, create a new node at the midpoint of ij; if the ij edge has been subdivided, the ij edge is an active edge, connect the existing nodes at the midpoint of the ij edge, and set the ij edge as an inactive edge; if the ij edge has not been subdivided, create a new node at the midpoint of ij, set the ij edge as an active edge, find another unit Eto that shares the ij edge, and determine the number of the ij edge in the Eto unit.

[0015] Preferably, for each unsubdivided unit shared by the active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits, including:

[0016] Check each active edge. There is only one unit among the two units sharing the active edge that is not subdivided, which is represented by Eto.

[0017] If the unit Eto needs to be subdivided into 4 subunits, subdivide the unit Eto;

[0018] Assuming that the ij edge of unit Eto is an active edge, the conditions under which unit Eto needs to be subdivided into 4 subunits include:

[0019] The jk edge and ki edge of unit Eto are both active edges. The middle nodes on the jk edge and ki edge have been established and connected to form 4 sub-units. The jk edge and ki edge are transformed into inactive edges.

[0020] One of the jk and ki edges of unit Eto is an active edge. Assuming that the jk edge is an active edge, the ki edge is an inactive edge. The midpoint of the jk edge already exists. Convert the ik edge to an inactive edge, create a new node at the midpoint of the ki edge, set the ki edge as an active edge, and connect the nodes on the three edges to form four sub-units.

[0021] The jk edge and ki edge of unit Eto are both inactive edges and the diagonal of the ij edge cannot be divided into two. A new node is established at the midpoint of each jk edge and ki edge, and they are connected to form 4 sub-units. The jk edge and ki edge are set as active edges.

[0022] Preferably, the calculation formula of the Peclet number is:

[0023]

[0024] Where Pe is the Peclet number, Δx and Δy are the characteristic lengths in the x and y directions respectively, Δt is the time step, and u x 、u y are the first groundwater velocity component and the second groundwater velocity component respectively; D x , D y : are the first hydrodynamic diffusion coefficient component and the second hydrodynamic diffusion coefficient component respectively.

[0025] Preferably, the time step restriction conditions include:

[0026] Among them, C r is the Courant number, q x is the single width flow in the x direction, n is the effective porosity, b is the thickness of the aquifer; the time step is:

[0027] Preferably, the time step restriction condition also includes:

[0028] in, β=1,2,…N e ; N e is the number of finite element meshes, are the first hydrodynamic diffusion coefficient component and the second hydrodynamic diffusion coefficient component corresponding to a certain β grid respectively; are the first groundwater velocity component and the second groundwater velocity component corresponding to a certain β grid respectively.

[0029] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0030] The present invention provides a finite element grid adaptive subdivision method for solving convection-based pollutant migration, including: identifying the convection-based finite element grid through a grid adaptive subdivision algorithm, and subdividing the convection-based finite element grid to obtain subdivided subunits; for each unsubdivided unit shared by an active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits; for the units where all the remaining active edges are located, the remaining active edges are connected to the diagonal vertices of the remaining active edges to form 2 subdivided subunits. The present invention can effectively solve the problems of vibration of numerical solutions and large numerical dispersion in traditional numerical models when solving convection-based solute migration problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0032] Figure 1 A flow chart of a method provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0033] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0034] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0035] Figure 1 A flow chart of a method provided by an embodiment of the present invention, such as Figure 1 As shown, the present invention provides a finite element grid adaptive subdivision method for solving convection-based pollutant migration, comprising:

[0036] Step 100: identifying a finite element mesh dominated by convection through a mesh adaptive subdivision algorithm, and subdividing the finite element mesh dominated by convection to obtain subdivided sub-units;

[0037] Step 200: For each unsubdivided unit shared by each active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits;

[0038] Step 300: For the cells where all the remaining active edges are located, connect the remaining active edges and the diagonal vertices of the remaining active edges to form two subdivided sub-cells;

[0039] Step 400: After the subdivision unit is divided, the subdivision unit has the same unit number as the parent unit, so that the subdivision unit references the parameters of the parent unit.

[0040] After the central difference format is used in this embodiment to discretize the time domain, theoretically, the difference format is unconditionally convergent with respect to the time step Δt, but the size of Δt affects the calculation accuracy and is related to the oscillation of the solution and numerical dispersion. The Peclet number (Pe) is defined as follows:

[0041]

[0042] Among them, u x 、u y : Groundwater velocity component (m·d -1 );D x , D y : Hydrodynamic dispersion coefficient component (m 2 ·d -1 )

[0043] It can be seen that the size of Pe is not only related to the spatial step size Δx and Δy, but also to the time step size Δt. In PRB-Trans, a grid adaptive subdivision algorithm is used to control the size of the spatial step size. In order to ensure the condition of Pe < 2, Δt must also be limited to be small enough.

[0044] The first restriction on the time step is to satisfy the Courant number condition. The Courant number is defined as:

[0045]

[0046] In the formula, C r is the Courant number, C r =1 is the most commonly used standard, q x is the single width flow in the x direction (m 2 ·d -1 ), n is the effective porosity, and b is the thickness of the aquifer (m). In order to meet this condition, the time step of the two-dimensional problem is:

[0047]

[0048] Wherein, Δx and Δy are characteristic lengths (maximum lengths in the directions) in the x and y directions, respectively.

[0049] The second restriction on the time step is to satisfy the Peclet number restriction. The Peclet number for the time step is:

[0050]

[0051] Usually Pe=2. In order to meet this condition, the time step is:

[0052]

[0053] Among them, N e is the number of finite element meshes, is the hydrodynamic diffusion coefficient component corresponding to a certain β grid (m 2 ·d -1 ); is the groundwater velocity component corresponding to a certain β grid (m·d -1 ).

[0054] In actual calculation, the smallest time step among all restrictions is taken as the calculation time step.

[0055] Specifically, the specific subdivision steps of this embodiment include:

[0056] ① Check whether each grid is mainly convection (Pe ≥ 2). For cells that are mainly convection, subdivide them into 4 sub-cells.

[0057] Check each unit ke, if its Pe ≥ 2, then subdivide it into 4 subunits. In the subdivision, it is necessary to check the properties of the three edges of ke, assuming that the three edges of the ke unit are ij = 1, jk = 2, ki = 3. Considering the ij edge, there are three cases:

[0058] 1. If edge ij is a boundary edge, create a new node at the midpoint of ij. End.

[0059] 2. If the ij edge has been subdivided, Sides.Active of that edge must be True. Connect the existing node at the midpoint of the ij edge and set Sides.Active of the ij edge to False. End

[0060] 3. If the ij edge is not subdivided, create a new node at the midpoint of ij, set Sides.Active = True for the edge, find another unit Eto that shares the ij edge, and determine the number of the ij edge in the Eto unit. End.

[0061] Therefore, among all the units, the ones that have not been subdivided are definitely not convection-dominated units.

[0062] After a unit is subdivided into subunits, the subunits have the same unit number as the parent unit in order to reference the parameters of the parent unit, such as permeability, dispersivity, etc.

[0063] ② Check the unsubdivided cells shared by each active side (Sides.Active = True), and only subdivide the shared cells that need to be subdivided into 4 sub-cells

[0064] Check each active edge (Sides.Active = True). There is only one unit among the two units sharing the edge that is not subdivided, represented by Eto. If Eto needs to be subdivided into 4 subunits, subdivide the unit, otherwise end.

[0065] Assuming that the ij edge of the Eto unit is an active edge, the conditions under which Eto needs to be subdivided into 4 subunits are as follows:

[0066] 1. The jk and ki edges of Eto are both active edges. At this point, the middle nodes on the jk and ki edges have been established, connecting to form 4 sub-units, and the jk and ki edges are transformed into inactive edges, Sides.Active = False.

[0067] Finish

[0068] 2. One of the jk and ki edges of Eto is an active edge. Assuming that the jk edge is an active edge, the ki edge is an inactive edge. At this time, the midpoint of the jk edge already exists, and the ik edge is converted to an inactive edge Sides.Active = False. Create a new node at the midpoint of the ki edge, set the ki edge to an active edge (Sides.Active = True), and connect the nodes on the three edges to form 4 sub-units. End

[0069] 3. The jk and ki edges of Eto are both inactive edges and the diagonal of the ij edge cannot be divided into two. At this time, a new node is established at the midpoint of the jk and ki edges, and they are connected to form 4 sub-units. The jk and ki edges are set as active edges. End.

[0070] If Eto is subdivided into 4 subunits, the ij edge is converted into an inactive edge, otherwise it remains active.

[0071] This process will generate new active edges, which will also be checked and subdivided as necessary during this process.

[0072] ③The remaining active edges are connected to their diagonal vertices to form two subdivision units

[0073] After the second step, all the remaining cells where the active edges are located are cells that do not need to be divided into 4 parts. Just connect the nodes in the active edges with the diagonal vertices to generate 2 subdivision cells. End.

[0074] After the above three steps, the first subdivision of the grid is completed and a new grid is formed. Repeat the above three steps until all units meet Pe<2.

[0075] The structure of the new grid formed by the above process is usually not as good as that of the base grid. This is because the vertex angles of some cells of the base grid are divided into two angles with smaller angles. We know that the structure is best when the three internal angles of the cell are all 60 degrees. If the original vertex angle is divided into two angles close to 60 degrees, the grid structure is improved. However, for the grid divided by the Delaunay triangulation method, it has an optimized topological structure, and the internal angles of most cells are close to 60 degrees. After subdivision, some 60-degree internal angles are divided into two angles of nearly 30 degrees, making the grid structure worse. We can also refine the grid while keeping the original grid structure unchanged. At this time, as long as there is one cell that does not meet the condition of Pe<2, each grid needs to be divided into 4 sub-grids, which increases the number of cells in the new grid by 4 times, which is obviously not economical.

[0076] Optionally, during mesh refinement, to prevent memory fragmentation, do not request memory one cell at a time, but request memory for a certain number of cells and nodes at a time, for example, 100 cells and 100 nodes at a time.

[0077] The beneficial effects of the present invention are as follows:

[0078] (1) The present invention takes reducing the Peclet number as a starting point, and autonomously identifies the finite element grid dominated by convection (the judgment standard is Pe ≥ 2) through the grid adaptive subdivision algorithm. Under the premise of controlling the time step, the spatial step is automatically subdivided to ensure that Pe of each unit is less than 2, thereby converting the pollutant migration problem dominated by convection into a pollutant migration problem dominated by diffusion (the judgment standard is Pe < 2) that is better solved, thereby improving the stability of the numerical solution;

[0079] (2) The finite element mesh is initially divided manually, so the Pe of each unit mesh is unknown. The present invention automatically identifies and subdivides the mesh based on the initial division.

[0080] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0081] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A finite element mesh adaptive subdivision method for solving convection-based pollutant migration, characterized in that: include: The finite element mesh dominated by convection is identified by a mesh adaptive subdivision algorithm, and the finite element mesh dominated by convection is subdivided to obtain subdivided sub-units; For each unsubdivided unit shared by the active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits; For the cells where all the remaining active edges are located, connect the remaining active edges and the diagonal vertices of the remaining active edges to form two subdivided sub-cells; After the subdivision unit is divided, the subdivision unit has the same unit number as the parent unit, so that the subdivision unit can reference the parameters of the parent unit; The identification methods for convection-dominated finite element meshes include: Determine whether the Peclet number is greater than or equal to 2, and if so, determine the mesh as the convection-dominated finite element mesh; The calculation formula of the Peclet number is: Where Pe is the Peclet number, Δx and Δy are the characteristic lengths in the x and y directions respectively, Δt is the time step, and u x 、u y are the first groundwater velocity component and the second groundwater velocity component respectively; D x , D y : are the first hydrodynamic diffusion coefficient component and the second hydrodynamic diffusion coefficient component respectively.

2. The finite element mesh adaptive subdivision method for solving convection-based pollutant migration according to claim 1 is characterized in that: The process of subdividing a convection-dominated finite element mesh also includes: Check the properties of the three edges of the finite element mesh dominated by convection, assuming that the three edges of the finite element mesh dominated by convection are ij=1, jk=2, ki=3 respectively; consider the ij edge, if the ij edge is a boundary edge, create a new node at the midpoint of ij; if the ij edge has been subdivided, the ij edge is an active edge, connect the existing nodes on the midpoint of the ij edge, and set the ij edge as an inactive edge; if the ij edge has not been subdivided, create a new node at the midpoint of ij, set the ij edge as an active edge, find another unit Eto that shares the ij edge, and determine the number of the ij edge in the Eto unit.

3. The finite element mesh adaptive subdivision method for solving convection-based pollutant migration according to claim 1 is characterized in that: For each unsubdivided unit shared by the active edge, only the shared unit that needs to be subdivided into 4 subunits is subdivided to obtain subdivided subunits, including: Check each active edge. There is only one unit among the two units sharing the active edge that is not subdivided, which is represented by Eto. If the unit Eto needs to be subdivided into 4 subunits, subdivide the unit Eto; Assuming that the ij edge of unit Eto is an active edge, the conditions under which unit Eto needs to be subdivided into 4 subunits include: The jk edge and ki edge of unit Eto are both active edges. The middle nodes on the jk edge and ki edge have been established and connected to form 4 sub-units. The jk edge and ki edge are transformed into inactive edges. One of the jk and ki edges of unit Eto is an active edge. Assuming that the jk edge is an active edge, the ki edge is an inactive edge. The midpoint of the jk edge already exists. Convert the ik edge to an inactive edge, create a new node at the midpoint of the ki edge, set the ki edge as an active edge, and connect the nodes on the three edges to form four sub-units. The jk edge and ki edge of unit Eto are both inactive edges and the diagonal of the ij edge cannot be divided into two. A new node is established at the midpoint of each jk edge and ki edge, and they are connected to form 4 sub-units. The jk edge and ki edge are set as active edges.

4. The finite element mesh adaptive subdivision method for solving convection-based pollutant migration according to claim 1 is characterized in that: The time step constraints include: Among them, C r is the Courant number, q x is the single width flow in the x direction, n is the effective porosity, b is the thickness of the aquifer; the time step is:

5. The finite element mesh adaptive subdivision method for solving convection-based pollutant migration according to claim 4 is characterized in that: The time step constraints also include: in, β=1,2,…N e ; N e is the number of finite element meshes, are the first hydrodynamic diffusion coefficient component and the second hydrodynamic diffusion coefficient component corresponding to a certain β grid respectively; are the first groundwater velocity component and the second groundwater velocity component corresponding to a certain β grid respectively.

Images