An efficient delaunay surface triangulation method

Abstract

The application belongs to the technical field of three-dimensional fluid mechanics numerical solution, and relates to a high-efficiency Delaunay surface triangular grid generation method. First, the geometric surface, geometric edge and geometric vertex and other information of a to-be-processed model are obtained through a geometric modeling kernel; then, the geometric edge of the model is discretized to construct a one-dimensional discrete grid; then, based on the topological relation of the geometric surface, boundary discrete points are gradually generated and mapped to a parameter plane; subsequently, the initial triangulation grid is generated by using a constrained Delaunay method, and internal points are inserted one by one to avoid searching and multiple traversing of base units. The application can significantly improve the efficiency of grid generation, reduce the calculation cost, and meet the high-quality numerical simulation demand of complex geometric shapes.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of three-dimensional fluid mechanics numerical solution, and relates to a high-efficiency Delaunay surface triangular grid generation method. BACKGROUND

[0002] Computational fluid dynamics (CFD) has been widely used in the fields of automobile manufacturing, civil engineering, environmental engineering, shipbuilding industry and aviation industry, etc. It helps to explain and understand the results of theory and experiment, and is an indispensable method for fluid mechanics analysis. In the early stage, due to the limitation of computer level, CFD was limited to the solution of some simple problems. With the development of computer level, today's CFD can solve a large number of complex three-dimensional flow fields, although a large amount of manpower and computer resources are still needed, but this solving method has been widely used in industrial equipment.

[0003] With the progress of industrial technology, fluid dynamics puts forward higher requirements for the precision of numerical algorithm, so high-precision numerical simulation method is needed. Grid generation is an important part of numerical calculation, especially in the field of computational fluid dynamics (CFD), the requirement for grid generation speed is extremely strict. Most of the total time of CFD calculation is often concentrated in the grid generation link. As the starting point of CFD numerical simulation, grid generation is not only the most basic step, but also the most important step. The degree of automation and quality of grid generation directly affect the period, precision and efficiency of CFD calculation. With the continuous expansion of the application field of CFD and the increase of problem complexity, the existing grid generation technology has been insufficient in dealing with high-quality numerical simulation of complex shape. Therefore, how to develop highly automated and efficient grid generation technology has become a research hotspot.

[0004] At present, the mainstream Delaunay surface triangular grid generation serial program usually includes the following five steps: model input, boundary discretization, surface discretization, grid optimization and grid merging. Among them, surface discretization is the longest time-consuming step. The traditional Delaunay surface triangular grid generation method usually adopts the method of taking the current triangulation as the background grid, generating internal points in batches according to the grid discretization size requirements, and inserting them one by one; this process is looped until the grid density or geometric quality meets the convergence condition. Such a point generation and insertion algorithm inevitably needs to traverse the current triangulation grid multiple times. In the process of internal point insertion, algorithms such as quadtree search algorithm and KD tree search algorithm are often used to speed up the search for cavities. These algorithms locate the base unit containing the point to be inserted, and find the cavity by using the adjacency relationship of the base unit triangle. However, due to the need to repeatedly iterate and search for the base unit, this process consumes a large amount of computing resources, thereby restricting the engineering application of high-precision algorithms. SUMMARY

[0005] In order to solve the problems or deficiencies, the present application provides a high-efficiency Delaunay surface triangular mesh generation method, and develops a more efficient internal point generation and insertion technology, reduces the iteration cycle, and improves the mesh generation efficiency under the premise of ensuring high precision.

[0006] A high-efficiency Delaunay surface triangular mesh generation method comprises the following steps:

[0007] A. A target structure is modeled into a three-dimensional entity, and is added to a model list, and a geometric face list, a geometric edge list, a geometric vertex list and their topological relations of the contour of the target structure are obtained through a geometric modeling kernel.

[0008] B. According to the size requirement of mesh discretization, all geometric edges in the geometric edge list obtained in step A are discretized, and a one-dimensional discrete mesh of the target model is generated according to the discrete information.

[0009] C. According to the geometric face, geometric edge, geometric vertex and their topological relations obtained in step A, and in combination with the one-dimensional discrete mesh obtained in step B, the discrete points of the geometric edges contained in a single geometric face are combined into the boundary discrete points of the geometric face, and the connection relations between the boundary discrete points are obtained.

[0010] D. In combination with the parameter domain information of the current geometric face, a parameter plane is constructed, and the boundary discrete points obtained in step C are mapped to the parameter plane to obtain the boundary discrete points on the parameter plane of the current geometric face and the mapping relation.

[0011] E. Based on the boundary discrete points on the parameter plane of the current geometric face obtained in step D, an incremental point insertion algorithm is adopted to insert the boundary points and perform boundary recovery to construct an initial triangulation mesh of the current geometric face.

[0012] F. Based on the initial triangulation mesh of the current geometric face obtained in step E, the initial triangulation mesh is traversed, the internal points are generated and inserted, and the initial triangulation mesh is updated to obtain a triangulation mesh of the current geometric face meeting the mesh discretization size requirement.

[0013] G. Based on the triangulation mesh of the single geometric face generated in step F, mesh quality optimization is performed to obtain an optimized triangulation mesh of the single geometric face.

[0014] H. Steps C to G are repeated until all single geometric faces in the geometric face list have been traversed, and the triangulation meshes of all single geometric faces in the geometric face list are obtained.

[0015] I. Based on the topological relationship obtained in step A, the mapping relationship in step D and the triangular mesh obtained in step H, the triangular mesh of all single geometric faces is reflected from the parameter plane back to the three-dimensional real space to obtain the surface triangular mesh of the target structure contour.

[0016] The application firstly reads a CAD model to be meshed, obtains the faces, edges and vertices of the model and establishes an index relationship, then discretizes the edges, performs Delaunay triangular meshing on each surface, optimizes the quality, and finally combines the whole surface mesh.

[0017] In summary, in the process of generating the mesh on the parameter plane of a single geometric face, the application generates the insertion points one by one, uses the triangle of the insertion points as the base unit, thereby avoiding the search of the base unit and avoiding multiple traversals of the current triangular mesh, and in the process of generating the mesh of a single geometric face, each triangular mesh is only traversed once. Therefore, a large amount of calculation overhead is saved, the memory overhead is saved, and the efficiency of mesh generation is improved. BRIEF DESCRIPTION OF DRAWINGS

[0018] Figure 1 is a flowchart of the application;

[0019] Figure 2 is a CAD model diagram built in step A of the embodiment;

[0020] Figure 3 is a comparison diagram of the calculation performance of the embodiment and the MeshSim software of Simmetrix company. DETAILED DESCRIPTION

[0021] The technical solutions of the application will be described in detail below with reference to the drawings and embodiments.

[0022] An efficient Delaunay surface triangular mesh generation method, with reference to Figure 1 , includes the following steps:

[0023] A. A target aircraft structure is modeled as a three-dimensional entity, and is added to a model list. A geometric modeling kernel is used to obtain a geometric face list, a geometric edge list, a geometric vertex list and a topological relationship therebetween of the aircraft contour.

[0024] This embodiment establishes a X-51 aircraft model as shown in Figure 2 by using an ACIS geometric kernel. An API function is used to add the X-51 aircraft model to the model list. A geometric modeling kernel is used to obtain a geometric face list, a geometric edge list, a geometric vertex list and a topological relationship therebetween of the X-51 aircraft contour in the model list.

[0025] In step A, based on the information of geometric faces, geometric edges, and geometric vertices in the geometric face list, geometric edge list, and geometric vertex list, the geometric edges and geometric vertices contained in each geometric face and the geometric vertices contained in each geometric edge are obtained through the ACIS geometric kernel.

[0026] B. According to the size requirement of the grid discretization, all geometric edges in the geometric edge list obtained in step A are discretized, and a one-dimensional discrete grid of the target model is generated according to the discretization information.

[0027] Further, based on the geometric edge list obtained in step A, the distance of each geometric edge from the starting point of the curve is calculated based on the sampling size using a numerical integration algorithm to calculate the arc length and correct the sampling distance to obtain the distance of the sampling point from the starting point of the curve:

[0028]

[0029] where s i is the distance of the i-th sampling point from the starting point of the curve, u i is the parameter coordinate of the i-th sampling point, r(t) is the curve equation, m is the total number of sampling points, and l is the arc length of the curve.

[0030] Solving equation (1) integral, the parameter coordinate of the sampling point is obtained (set the sampling point), and then the geometric edge is discretized in one-dimensional parameter space according to the size information of the sampling point to obtain the one-dimensional parameter coordinate of the discrete point.

[0031] Then, through the geometric modeling kernel, based on the one-dimensional parameter coordinate of the discrete point and the information of the geometric edge, the three-dimensional real space coordinates of the discrete point are obtained through the ACIS geometric kernel.

[0032] After all the geometric edges are discretized, a one-dimensional grid is generated according to the three-dimensional real space coordinates of the discrete points of all the geometric edges, thereby obtaining the connection relationship between the discrete points of all the geometric edges, and a one-dimensional discrete grid is obtained.

[0033] C. According to the geometric face, geometric edge, geometric vertex, and their topological relationship obtained in step A, and combining the one-dimensional discrete grid obtained in step B, the discrete points of the geometric edges contained in a single geometric face are combined into the boundary discrete points of the geometric face, and the connection relationship between the boundary discrete points is obtained.

[0034] D. Combining the parameter domain information of the current geometric face, a parameter plane is constructed, and the boundary discrete points obtained in step C are mapped to the parameter plane to obtain the boundary discrete points on the parameter plane of the current geometric face and the mapping relationship.

[0035] Further, according to the parameter domain information of the single geometric face, the parameter range of the parameter domain of the geometric face is obtained through the ACIS geometric kernel, and a parameter plane is constructed. Then, according to the boundary discrete points of the single geometric face obtained in step C and the geometric face information, the boundary discrete points of the single geometric face are mapped onto the parameter plane through the ACIS geometric kernel, and it is recorded that the boundary discrete points on the parameter plane are mapped from which boundary discrete points, so as to obtain the boundary discrete points on the parameter plane of the single geometric face and the mapping relationship.

[0036] E. Based on the boundary discrete points on the parameter plane of the current geometric face obtained in step D, an initial triangulation grid of the current geometric face is constructed by using a constrained Delaunay method.

[0037] Further, according to the parameter plane constructed in step D, a super-triangle containing all the boundary discrete points on the parameter plane of the single geometric face is constructed. The boundary discrete points on the parameter plane of the single geometric face are inserted one by one by using the Bowyer-Watson incremental point insertion algorithm with the super-triangle as the initial triangulation. Then, the boundary is recovered according to the connection relationship between the discrete points obtained in step B, so as to obtain a triangulation grid containing external cells.

[0038] Based on the triangulation grid containing external cells, a triangle containing the vertices of the super-triangle is searched, and the attribute of the triangle is set as invalid. Then, taking the triangle as a starting point, the attributes of the surrounding triangles are set according to the triangle adjacency relationship: if both vertices of the common edge of two triangles are boundary discrete points on the parameter plane and they are connected, then the attribute is inverted; otherwise, the attribute is the same. Until the effective attributes of all the triangles in the triangulation grid containing external cells are set, the triangles with invalid attributes are deleted, and an initial triangulation grid of the single geometric face is obtained.

[0039] F. Based on the initial triangulation grid of the current geometric face obtained in step E, the initial triangulation grid is traversed to generate and insert internal points, and the initial triangulation grid is updated, so as to obtain a triangulation grid of the current geometric face meeting the grid discrete size requirement.

[0040] Further, according to the initial triangulation grid of the single geometric face obtained in step E, the triangles in the initial triangulation grid are traversed to determine whether the current triangle meets the grid discrete size requirement: if not, the parameter coordinates of the to-be-inserted point P are calculated according to the parameter coordinates of the three vertices of the current triangle and the Riemann metric of the three vertices:

[0041] (P u ,P v )(2)

[0042] P u =(size 1u ×P1u +size 2u ×P 2u +size 3u ×P 3u ) / (size 1u +size 2u +size 3u ) (3)

[0043] P v =(size 1v ×P 1v +size 2v ×P 2v +size 3v ×P 3v ) / (size 1v +size 2v +size 3v ) (4)

[0044] where P u is the u-direction coordinate of the insertion point P, where P v is the v-direction coordinate of the insertion point P, P 1u , P 2u , P 3u are the u-direction coordinates of the three vertices, respectively, P 1v , P 2v , P 3v are the v-direction coordinates of the three vertices, respectively, size 1u size 2u , size 3u are the u-direction grid discrete sizes of the three vertices of the current triangle, respectively, size 1v size 2v , size 3v are the v-direction grid discrete sizes of the three vertices of the current triangle, respectively.

[0045] According to the combined Delaunay criterion applicable to the parametric plane:

[0046]

[0047] where a p is the Delaunay measure of the insertion point, a vi is the Delaunay measure of the three vertices v1, v2, v3 of the current triangle.

[0048] Take the current triangle as the base unit, search all triangles containing the insertion point according to the adjacency relationship of the triangles, form a cavity, insert the to-be-inserted point, and generate new triangles; the destroyed triangle sets its attribute to invalid, and updates the current triangulation grid and the adjacency relationship of the internal triangles.

[0049] The updated new triangle is traversed to determine whether it meets the grid discrete size requirement, and if not, the generation and insertion of internal points are continued until the current triangulation grid is completed and meets the grid discrete size requirement, and all triangles with invalid attributes are deleted.

[0050] In the above process, after traversal is completed, it can be ensured that all triangles meet the grid discrete size requirement, so the triangulation grid at this time is the triangulation grid of a single geometric surface that meets the grid discrete size requirement.

[0051] G. Based on the single geometric surface triangulation grid finally generated in step F, a gradient descent method is used to optimize the grid quality to obtain an optimized single geometric surface triangulation grid.

[0052] Further, according to the single geometric surface triangulation grid generated in step F, the square of the difference between the actual grid length and the ideal grid length required by the grid discrete size requirement is taken as the energy function:

[0053]

[0054] Where E is the energy function, N is the total number of edges of the single geometric surface triangulation grid generated in step F, L i is the length of the i-th edge, and L i0 is the ideal grid length of the i-th edge required by the grid discrete size requirement.

[0055] The gradient of the energy function with respect to the actual grid length is calculated, and the grid point positions except the boundary discrete points are moved in the gradient direction until the energy function meets the set convergence condition to obtain the optimized single geometric surface triangulation grid.

[0056] H. Repeat steps C to G until all single geometric surfaces in the geometric surface list have been traversed to obtain the triangulation grid of all single geometric surfaces in the geometric surface list.

[0057] I. Based on the topological relationship obtained in step A, the mapping relationship in step D, and the triangulation grid obtained in step H, all single geometric surface triangulation grids are reflected from the parameter plane back to the three-dimensional real space to obtain the surface triangulation grid of the target structure contour.

[0058] Based on the topological relationship obtained in step A and the triangular mesh obtained in step H, the internal vertex parameter coordinates of the triangular mesh of the single geometric surface generated in the parameter plane are mapped into the three-dimensional real space coordinates by the ACIS geometry kernel. The boundary discrete points of the triangular mesh are reflected to the three-dimensional real space coordinates by the mapping relationship in step D, and the triangular mesh of the aircraft contour surface is obtained.

[0059] Figure 3 The embodiment and the comparative example generate the same number of grids, and the time consumption comparison results show that the application can greatly improve the grid generation efficiency.

[0060] As can be seen from the above embodiments and the comparative example, the application generates the insertion points one by one, generates the triangle of the insertion point as the base unit, avoids the search of the base unit, and avoids multiple traversals of the triangulation mesh of the current geometric surface, thereby ensuring that the mesh meets the size requirements and the mesh quality, and improving the efficiency of the Delaunay surface triangular mesh generation. The application provides a basis for the high-precision algorithm engineering application of computational fluid dynamics and computational electromagnetics in the fields of automobile manufacturing, civil engineering, environmental engineering, shipbuilding industry and aviation industry, and helps the interpretation, understanding theory, experiment and design of related industries.

Claims

1. An efficient Delaunay surface triangulation mesh generation method, characterized in that, The method comprises the following steps: A. Three-dimensional entity modeling is performed on the target aircraft structure, and is added to a model list, an aircraft model is established through a geometric modeling kernel, a geometric face list, a geometric edge list, a geometric vertex list of the target structure contour profile, and a topological relationship therebetween are obtained; B. According to the size requirement of grid discretization, all geometric edges in the geometric edge list obtained in step A are discretized, and a one-dimensional discrete grid of the target model is generated according to the discretization information; C. According to the geometric face, geometric edge, geometric vertex and topological relationship thereof obtained in step A, and in combination with the one-dimensional discrete grid obtained in step B, the discrete points of the geometric edges contained in a single geometric face are combined into the boundary discrete points of the geometric face, and the connection relationship between the boundary discrete points is obtained; D. In combination with the parameter domain information of the current geometric face, a parameter plane is constructed, and the boundary discrete points obtained in step C are mapped to the parameter plane to obtain the boundary discrete points on the parameter plane of the current geometric face and the mapping relationship; E. Based on the boundary discrete points on the parameter plane of the current geometric face obtained in step D, an incremental point insertion algorithm is adopted to insert and perform boundary recovery on the boundary points to construct an initial triangulation grid of the current geometric face; F. Based on the initial triangulation grid of the current geometric face obtained in step E, the initial triangulation grid is traversed, internal points are generated and inserted, and the initial triangulation grid is updated to obtain a triangulation grid of the current geometric face meeting the grid discretization size requirement; Specifically, according to the initial triangulation grid of a single geometric face obtained in step E, the triangles in the initial triangulation grid are traversed, and it is judged whether the current triangle meets the grid discretization size requirement: if not, the parameter coordinates of the to-be-inserted point P are calculated according to the parameter coordinates of the three vertices of the current triangle and the Riemann metric of the three vertices: (P u ,P v )(2) P u = (size 1u × P 1u + size 2u × P 2u + size 3u × P 3u ) / (size 1u + size 2u + size 3u ) (3) P v = (size 1v × P 1v + size 2v × P 2v + size 3v × P 3v ) / (size 1v + size 2v + size 3v ) (4) where P u is the u-direction coordinate of the insertion point P, where P v is the v-direction coordinate of the insertion point P, where P 1u , P 2u , P 3u , are the u-direction coordinates of the three vertices, respectively, where P 1v , P 2v , P 3v , are the v-direction coordinates of the three vertices, respectively, where size 1u size 2u size 3u are the u-direction grid discrete sizes of the three vertices of the current triangle, respectively, where size 1v size 2v size 3v are the v-direction grid discrete sizes of the three vertices of the current triangle, respectively. According to the combined Delaunay criterion applicable to the parameter plane: where α p is the Delaunay metric of the insertion point, is the Delaunay metric of the three vertices v1, v2, v3 of the current triangle. Taking the current triangle as a base unit, in combination with the adjacency relationship of the triangle, all triangles containing the to-be-inserted point are searched, a cavity is formed, and the to-be-inserted point is inserted to generate a new triangle; the destroyed triangle sets its attribute to invalid, and the current triangulation grid and the adjacency relationship of the internal triangle are updated; The updated new triangle is traversed to judge whether it meets the grid discretization size requirement, if not, the generation and insertion of the internal point are continued, and when the current triangulation grid is completely traversed and meets the grid discretization size requirement, all triangles with the invalid attribute are deleted; G. Based on the triangulation grid of the single geometric face finally generated in step F, grid quality optimization is performed to obtain the optimized triangulation grid of the single geometric face; H. Steps C to G are repeated until all single geometric faces in the geometric face list have been traversed, and the triangulation grid of all single geometric faces in the geometric face list is obtained; I. Based on the topological relationship obtained in step A, the mapping relationship in step D and the triangulation grid obtained in step H, the triangulation grid of all single geometric faces is reflected from the parameter plane back to the three-dimensional real space to obtain the surface triangulation grid of the target structure contour profile.

2. The efficient Delaunay surface triangular mesh generation method of claim 1, wherein: the discretization of the geometric edge in step B is specifically: based on the sampling size, the arc length is calculated by using numerical integral algorithm, and the sampling distance is corrected to obtain the distance of the sampling point from the starting point of the curve; where s i is the distance of the i-th sampling point from the start of the curve, u i is the parametric coordinate of the i-th sampling point, r(t) is the curve equation, m is the total number of sampling points, and l is the arc length of the curve. the integral of formula (1) is solved to obtain the parameter coordinate of the sampling point, and the geometric edge is discretized in one-dimensional parameter space according to the size information of the sampling point to obtain the one-dimensional parameter coordinate of the discrete point; then, based on the one-dimensional parameter coordinate of the discrete point and the information of the geometric edge, the three-dimensional real space coordinate of the discrete point is obtained through the ACIS geometric kernel; after the discretization of all geometric edges is completed, a one-dimensional mesh is generated according to the three-dimensional real space coordinates of the discrete points of all geometric edges, so as to obtain the connection relationship between the discrete points of all geometric edges, and a one-dimensional discrete mesh is obtained.

3. The method of claim 1, wherein the method is efficient. step D is specifically: based on the parameter domain information of the single geometric face, the parameter range of the parameter domain of the geometric face is obtained through the ACIS geometric kernel to construct a parameter plane; then, based on the boundary discrete points and the geometric face information of the single geometric face obtained in step C, the boundary discrete points are mapped onto the parameter plane through the ACIS geometric kernel, and it is recorded that the boundary discrete points on the parameter plane are mapped from which boundary discrete points, so as to obtain the boundary discrete points on the parameter plane of the single geometric face and the mapping relationship.

4. The method of claim 1, wherein the method is efficient. step E adopts the constrained Delaunay method to insert the boundary points.

5. The method of claim 1, wherein the method is efficient. step E is specifically: based on the parameter plane constructed in step D, a super-triangle containing all the boundary discrete points on the parameter plane of the single geometric face is constructed; the boundary discrete points on the parameter plane of the single geometric face are inserted one by one by using the Bowyer-Watson incremental point insertion algorithm with the super-triangle as the initial triangulation; then, the boundary is recovered according to the connection relationship between the discrete points obtained in step B to obtain a triangulation mesh containing external elements; based on the triangulation mesh containing external elements, a triangle containing the vertices of the super-triangle is searched, the attribute of the triangle is set as invalid; and taking the triangle as the starting point, the attributes of the surrounding triangles are set according to the triangle adjacency relationship: if both vertices of the common edge of two triangles are boundary discrete points on the parameter plane and they are connected, the attribute is reversed; otherwise, the attribute is the same; until the valid attributes of all triangles in the triangulation mesh containing external elements are set, the triangles with invalid attributes are deleted, and the initial triangulation mesh of the single geometric face is obtained.

6. The method of claim 1, wherein the method is efficient. step G adopts the gradient descent method to optimize the mesh quality: based on the triangular mesh of the single geometric face generated in step F, the square of the difference between the actual mesh edge length and the ideal mesh edge length required by the mesh discretization size is taken as the energy function: where E is an energy function, N is the total number of edges of the triangulated mesh of the individual geometric face generated by step F, L i is the length of the i-th edge, L i0 is the ideal mesh edge length required by the mesh discretization size requirement of the i-th edge; the gradient of the energy function with respect to the actual mesh edge length is calculated, and the grid point positions except the boundary discrete points are moved in the gradient direction until the energy function meets the set convergence condition, and the optimized triangular mesh of the single geometric face is obtained.

Images