A hexahedral mesh adaptation method based on a posteriori error estimation

Abstract

The present invention provides a hexahedral mesh adaptation method based on a posteriori error estimation, which includes the following steps: (1) Determination of the area that needs to be refined based on the Riemann metric. Based on the finite element analysis results of the initial mesh, the Riemann metric field is determined, and then the local area that needs to be refined is determined. (2) Accurate local mesh refinement. By solving the objective function, all face sets that need to be refined are iteratively solved based on the graph segmentation method, and finally a layer insertion operation is performed. (3) Size-preserving geometric optimization. While ensuring quality optimization, the size is kept as consistent as possible with the requirements of the Riemann metric field. The present invention can perform physical-based high-quality mesh adaptation for the coarse-grained original hexahedral mesh, so that the final adaptive mesh can meet the needs of simulation analysis.

Description

Technical Field

[0001] The invention relates to a hexahedral grid self-adaptation method based on a posteriori error estimation. Background Art

[0002] With the development of computer technology and numerical methods, finite element analysis (FEA) has become increasingly widely used in engineering and plays an increasingly important role. Finite element analysis (FEA) first requires finite element meshing, and the efficiency and accuracy of the analysis directly depend on the type and quality of the finite element mesh. Therefore, the automatic generation of high-quality finite element meshes has long been an important research direction in the field of FEA. Compared with tetrahedral meshes, hexahedral meshes with the same number of meshes have the advantages of high computational accuracy and fast convergence. Therefore, they have always been more popular in the industry and are the ideal mesh type for FEA.

[0003] To improve analysis efficiency and quality, researchers have proposed adaptive meshing, which significantly enhances analysis efficiency while maintaining finite element analysis accuracy. Mesh adaptation methods can be broadly categorized into two types: geometry-based mesh adaptation and physics-based mesh adaptation. Compared to geometry-based mesh adaptation, physics-based mesh adaptation is theoretically more ideal, but also more difficult.

[0004] Compared with the research on triangular and tetrahedral mesh adaptation, there is currently less research on hexahedral mesh adaptation. The main existing method is the template encryption method for octree meshes, which performs local encryption on the points and edges marked as requiring encryption according to the template. However, since it is only applicable to octree meshes, it has great limitations. Other methods mainly use layer insertion operations for local encryption. Borden et al. proposed the Cleave And Fill algorithm, which performs layer insertion operations by finding a shortest path on the surface and extending it to the interior, but this method is only applicable to swept volume meshes. It should be pointed out that due to the strong structural constraints of the hexahedral mesh, the above methods cannot flexibly control the density of the local mesh. Summary of the Invention

[0005] In view of the deficiencies in the background technology, the present invention aims to provide a hexahedral mesh adaptation method based on a posteriori error estimation, which can accurately implement physics-based mesh adaptation and improve the efficiency of finite element analysis.

[0006] In order to achieve the above-mentioned purpose, the present invention adopts the following technical solutions:

[0007] The following steps are included:

[0008] A hexahedral mesh adaptation method based on a posteriori error estimation is characterized by comprising the following steps:

[0009] (1) Determination of the area that needs to be refined based on the Riemannian metric,

[0010] Based on the results of finite element analysis of the initial mesh, the Riemannian metric field is determined, thereby identifying the local areas that need to be refined;

[0011] (2) Accurate local mesh refinement,

[0012] By solving the objective function, all the face sets that need to be refined are iteratively solved based on the graph segmentation method, and finally the layer insertion operation is performed;

[0013] (3) Dimension-preserving geometric optimization,

[0014] While ensuring quality optimization, the size is kept as consistent as possible with the requirements of the Riemann metric field.

[0015] On the basis of adopting the above technical solutions, the present invention may also adopt or combine the following further technical solutions:

[0016] The step of determining the area to be refined based on the Riemannian metric includes the following steps:

[0017] Calculate the Riemann metric field based on finite element analysis results;

[0018] Identify areas that require refinement;

[0019] The formula for calculating the Riemann metric field based on the finite element analysis results is as follows:

[0020]

[0021] Where P is a grid point in the hexahedral grid, P i represents the vertices adjacent to P, H u Represents the second-order matrix at point P, namely the Riemann metric, u represents the displacement value at point P, u i Indicates P i The displacement value at Represents the first-order matrix at point P.

[0022] The step of determining the area to be refined comprises the following steps:

[0023] Determine the mesh edges that need to be refined and their target values;

[0024] Determine the elements that need to be refined and their refinement values;

[0025] In the step of determining the refined mesh edges and their target values, the length of each edge in the mesh in Riemann space is calculated using the following formula:

[0026]

[0027] in, Indicates PP i The length in the Riemannian metric field, M(P) represents the Riemannian metric at point P, PP i Represents point P and its adjacent point P i connected vectors.

[0028] The precise local refinement step includes the following steps:

[0029] Determine the initial refined face set based on graph partitioning;

[0030] Complete the refined face sets on the model surface, group the surface face sets according to geometry, and add the face sets that need the most refinement to the final refined face set;

[0031] Implement layer insertion operations on all refined face sets;

[0032] In the step of determining the initial refined face set based on graph segmentation, an optimization solution is performed based on the graph segmentation method [Boykov Y, Veksler O, Zabih R. Fast approximate energy minimization via graphcuts [J]. IEEE Transactions on pattern analysis and machine intelligence, 2001, 23(11): 1222-1239.], wherein the optimization objective function is as follows:

[0033]

[0034]

[0035]

[0036] Where P represents the set of nodes in the dual graph, f represents the ownership of each node p∈P in the dual graph, and it is marked as f p , D p represents the node p’s ownership cost, Q p，q represents the topological quality cost between nodes p and q, T v Indicates the target value of the unit to be refined, T v The value range is greater than or equal to 1.5 and less than or equal to max(T v ), Dif represents the square of the topological difference between the adjacent surface and the adjacent cut surface between the two units p and q.

[0037] In the step of determining the initial refined face set based on graph partitioning, the iteration is terminated when the solution obtained by the graph partitioning method contains only empty sets.

[0038] In the step of dimension-preserving geometric optimization, the optimal position of each grid point is obtained by minimizing the objective function, which is as follows:

[0039]

[0040] Among them, M is the number of hexahedral units associated with point v, and N is the number of tetrahedrons of the hexahedral decomposition associated with point v. Here, the value of N is 8. represents the shape metric function, Represents the length of the vector formed by points v and vi under the Riemannian metric.

[0041] The present invention provides a hexahedral mesh adaptation method based on a posteriori error estimation, enabling high-quality mesh adaptation for general hexahedral meshes. By solving the optimization objective function and employing a graph-based segmentation approach, the set of quadrilateral faces requiring refinement can be optimally determined. By incorporating a topological quality term into the optimization objective function, the quality of the insertion layer is improved, thereby balancing mesh quality requirements during the mesh adaptation process. A dimension-preserving geometric optimization method ensures that the final geometric optimization achieves both improved mesh quality and a robust mesh adaptation effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 Flowchart of the present invention; wherein, (a) input hexahedral mesh, (b) stress distribution of analysis results, (c) area to be encrypted, (d) quadrilateral face set to be inserted, and (e) optimized adaptive hexahedral mesh.

[0043] Figure 2 This is the interpretation of the Riemannian metric in geometric space; among them, (a) the Riemannian metric measures the lengths of the three directions of the local coordinate system in Cartesian space, and (b) when the metric in (a) is converted to the Riemannian space, the three directions become unit lengths.

[0044] Figure 3 Figure 1 is a flowchart of the process of determining the set of faces that need to be refined; (a) the area that needs to be refined, (b) the areas that need more refinement and the areas that need less refinement determined based on the graph segmentation algorithm, (c) the initial set of refined faces that are completely inside the mesh, and (d) the final set of refined faces that includes the boundary faces.

[0045] Figure 4 It is a schematic diagram of the topological difference between adjacent segmentation surfaces.

[0046] Figure 5Schematic diagram of the completion of the refined face set; (a) the surface mesh of the area to be refined is divided into four areas, and the average target value to be refined in each area is calculated; (b) the refined face set after completion. Specific implementation plan

[0047] The various steps of the present invention are further described below with reference to the accompanying drawings and examples.

[0048] 1. Calculation of Riemannian metric fields

[0049] In order to measure the size of the grid unit, there are certain constraints on its stretching degree in different directions. To this end, we introduce the Riemann metric. In three-dimensional space, each point has a 3×3 tensor matrix, represented by M, as follows Figure 2 As shown in (a). Let P be a point in space, M(P) be the Riemannian metric matrix at point P, then the ideal length PX in this space should be close to 1:

[0050]

[0051] A grid that meets the size requirements of the Riemannian metric is a unit grid in the Riemannian space, such as Figure 2 As shown in Figure (b).

[0052] The present invention obtains the Riemann metric by posterior error estimation. Let P represent a vertex in the hexahedral grid, u represent the numerical solution at point P, P i represents the vertex adjacent to P, u i Indicates P i The displacement value at . According to the Taylor expansion formula, we can get the following formula:

[0053]

[0054] The above formula is expanded as follows:

[0055]

[0056] You can get:

[0057]

[0058] Convert the above formula into the form of AX=B and use the least squares method to find the Riemann metric H u (P).

[0059] 2. Determine the units that need to be refined and their refinement values

[0060] First, the information about the mesh edges that need to be refined is transferred to the mesh faces that need to be refined. The refinement value of the mesh face that needs to be refined is equal to the average of the refinement values of the mesh edges that are topologically perpendicular to it. Next, the cells that need to be refined are determined, and their refinement value is equal to the average of the refinement values of all mesh faces in the cell. If the refinement value of a cell is greater than 1.5, it is considered to need to be refined.

[0061] 3. Determination of the initial refined face set based on graph segmentation

[0062] Because there are cells within a region that require varying degrees of refinement, we need to insert multiple layers with varying structures and positions to meet the adaptive density requirements. To this end, we adopt an iterative approach to graph segmentation, determining in each round the face set that needs to be refined and the number of times the face set needs to be refined.

[0063] During each iteration, we define an energy function to ensure that the cut sets obtained from the graph segmentation meet the requirements. By solving the problem of minimizing the energy function, we can obtain an optimal set of segmentation faces. The segmentation face set divides the area requiring refinement into two parts: area A that requires more refinement and area B that requires less refinement. The objective function corresponding to each round of graph segmentation is as follows:

[0064]

[0065] Where P represents the set of nodes in the dual graph, f represents the ownership of each node p∈P in the dual graph, and it is marked as f p There are two types of attributes here, namely area A which needs more refinement and area B which does not need more refinement. p Indicates that node p belongs to f p The price, Q p，q Measures the topological quality cost when a split surface is formed between nodes p and q.

[0066] By node attribution cost D p , we can distinguish between the area A that needs more refinement and the area B that does not need refinement. The attribution cost for each unit can be calculated according to the following formula:

[0067]

[0068] Where T v Indicates the target value of the unit to be refined, and its value range is (1.5, max(T v When the unit belongs to region A, the attribution cost decreases as the unit refinement value increases; when the unit belongs to region B, the attribution cost decreases as the unit refinement value decreases.

[0069] On the other hand, considering that the layer insertion operation on the refined face set will affect the quality of the mesh, we add the topological quality cost Q to the objective function p，q , to ensure mesh quality. Specifically, the topological quality cost between two adjacent hexahedral elements can be calculated according to the following formula:

[0070]

[0071] Dif measures the square of the topological difference between the adjacent faces and the adjacent cut faces between the two units p and q, as follows: Figure 4 As shown in the figure. When two units belong to the same area, the quadrilateral face between them is not a partition face set, so the topological quality cost is 0. When two units belong to different areas, the adjacent face between the current units is connected to a partition face at edge e. The two quadrilateral faces form a local partition for all hexahedral mesh units adjacent to e. This hexahedral unit set is partitioned on both sides of the two quadrilateral faces. When the number of units on both sides differs too much, the degree of the mesh edge will increase after the layer insertion operation and the mesh quality will deteriorate. Therefore, when performing local partitioning, we hope that the number of units on both sides is as equal as possible so that the partitioned quadrilateral face set is topologically smooth. For example Figure 4 As shown in , the number of cells on both sides of the split surface are 2 and 3 respectively, so the Dif value is 1. The adjacent surface has four edges, and the topological quality cost Q p，q It is equal to the sum of the squares of the Dif values of the cut faces adjacent to the four edges.

[0072] 4. Completion of refined face sets

[0073] The specific steps are as follows:

[0074] (1) Partition the surface quadrilateral mesh of area A that needs to be refined. Based on the geometric affiliation of the four vertices of the quadrilateral, we can determine the geometric affiliation of all quadrilateral faces. According to the geometric affiliation of the face, we can divide the quadrilaterals belonging to the same geometric face into the same area, such as Figure 5 As shown in part (a) of .

[0075] (2) Determine the final refined face set. Calculate the average target value for refinement of each surface quadrilateral mesh area, and select the quadrilateral face set with the highest average refinement value to be included in the final refined face set, because this surface area is the most in need of refinement, such as Figure 5 As shown in part (b) of .

[0076] (3) Determine the number of times the final refined face set needs to be refined. Calculate the average value T of the sum of the refinement values of all quadrilateral faces in the final refined face set. f , the number of refinements of this face set is equal to [T f]-1.

[0077] (4) Update the target value of mesh elements that need to be refined. In order to accurately solve the refined face set in the next round, we update the target value of mesh elements that need to be refined after the previous round of refinement. Since the refinement value on the mesh surface is directly affected by the layer operation, we first update the refinement value on the mesh surface. The affected mesh surface is the refined face set solved in this round. The new refinement value on the face set is equal to the original refinement value minus the number of times the face set needs to be refined. Finally, the refinement target value of all elements in the refined area is recalculated according to the specific implementation step 2.

[0078] 5. Dimension-preserving geometry optimization

[0079] The geometry optimization algorithm for the entire hexahedral mesh is shown below. Mesh points are divided into three categories: those located on geometric edges, those located on geometric faces, and those located within the body. For the first two categories of points, constraints are solved based on the geometric equations they represent, and the points are projected onto the nearest location on the corresponding geometric element. The geometry optimization algorithm terminates when the point's movement distance falls below a threshold or the number of iterations reaches a threshold.

[0080] When performing the geometric position optimization algorithm for each point, it is necessary to solve the optimal value of the objective function. The numerical optimization methods often used here are gradient descent method, Newton method and quasi-Newton method. The gradient descent method uses the minimum descent direction of the current position as the search direction and approaches the minimum value. The disadvantage is that the convergence speed slows down when approaching the extreme value, and multiple iterations are required to solve the problem. The Newton method is an iterative algorithm. Each iteration requires that the objective function is second-order differentiable and the Hessian matrix is reversible. Therefore, the requirements for the objective function are relatively high and the calculation is complex. The essential idea of the quasi-Newton method is to improve the defect of the Newton method that the inverse matrix of the Hessian matrix needs to be solved during iteration, thereby adopting the approximate matrix of the Hessian matrix to overcome the shortcomings of the Newton method such as the large amount of calculation. The present invention adopts the BFGS algorithm, which is a better numerical effect among the quasi-Newton methods.

Claims

1. A hexahedral mesh adaptation method based on a posteriori error estimation, characterized in that The following steps are involved: (1) Determination of the area that needs to be refined based on the Riemannian metric, Based on the results of finite element analysis of the initial mesh, the Riemannian metric field is determined, thereby identifying the local areas that need to be refined; (2) Accurate local mesh refinement, By solving the objective function, all the face sets that need to be refined are iteratively solved based on the graph segmentation method, and finally the layer insertion operation is performed; (3) Dimension-preserving geometric optimization, Ensure quality optimization while keeping the size as close as possible to the requirements of the Riemann metric field; The steps to determine the area that needs to be refined based on the Riemannian metric include the following steps: Calculate the Riemann metric field based on finite element analysis results; Identify areas that require refinement; The formula for calculating the Riemann metric field based on the finite element analysis results is as follows: Where P is a grid point in the hexahedral grid, P i represents the vertices adjacent to P, H u Represents the second-order matrix at point P, namely the Riemann metric, u represents the displacement value at point P, u i Indicates P i The displacement value at represents the first-order matrix at point P; The step of determining the area to be refined comprises the following steps: Determine the mesh edges that need to be refined and their target values; Determine the elements that need to be refined and their refinement values; In the step of determining the refined mesh edges and their target values, the length of each edge in the mesh in Riemann space is calculated using the following formula: in, Indicates PP i The length in the Riemannian metric field, M(P) represents the Riemannian metric at point P, PP i Represents point P and its adjacent point P i The connected vectors of ; The precise local refinement step includes the following steps: Determine the initial refined face set based on graph partitioning; Complete the refined face sets on the model surface, group the surface face sets according to geometry, and add the face sets that need the most refinement to the final refined face set; Implement layer insertion operations on all refined face sets; In the step of determining the initial refined face set based on graph segmentation, an optimization solution is performed based on the graph segmentation method, wherein the optimization objective function is as follows: Where P represents the set of nodes in the dual graph, f represents the ownership of each node p∈P in the dual graph, and it is marked as f p , D p represents the node p’s ownership cost, Q p,q represents the topological quality cost between nodes p and q, T v Indicates the target value of the unit to be refined, T v The value range is greater than or equal to 1.5 and less than or equal to max(T v ), Dif represents the square of the topological difference between the adjacent surface and the adjacent cut surface between the two units p and q.

2. A hexahedral mesh adaptation method based on a posteriori error estimation according to claim 1, characterized in that: In the step of determining the initial refined face set based on graph partitioning, the iteration is terminated when the solution of the graph partitioning method contains only empty sets.

3. A hexahedral mesh adaptation method based on a posteriori error estimation according to claim 1, characterized in that: In the step of dimension-preserving geometric optimization, the optimal position of each grid point is obtained by minimizing the objective function, which is as follows: Among them, M is the number of hexahedral units associated with point v, and N is the number of tetrahedrons of the hexahedral decomposition associated with point v. Here, the value of N is 8. represents the shape metric function, Indicates v and v i The length of the vector formed by the points in the Riemannian metric.

Images