CN114510775B - A 3D Space Curved Meshing Method for Complex Models - Google Patents

Abstract

The invention belongs to the technical field of aircraft design, aircraft numerical simulation and numerical calculation technology and grid division, and particularly relates to a three-dimensional space curve grid division method of a complex model. The invention provides a novel strategy for selecting control points and space curved grid points, which can well meet the grid division requirement of a complex model under the condition of not lifting the grid magnitude. The selection of control points and deformation points is reduced while the grid quality is ensured, a large number of position calculations in the grid division process are reduced, and the grid division strategy is optimized. And a basis function selection mode based on an included angle cosine value is provided, and on the premise of ensuring the grid quality, calculation required by deformation quantity can be greatly reduced by using a low-order polynomial basis function with a simple form, so that grid division time consumption is reduced, and calculation force requirement is reduced. The method effectively solves the problems of poor universality of grid division and poor grid quality in the existing complex model grid division process.

Description

A 3D Space Curved Meshing Method for Complex Models

technical field

The invention belongs to the technical fields of aircraft design, aircraft numerical simulation and numerical calculation technology, and grid division, and in particular relates to a complex model three-dimensional space curved grid division method.

Background technique

With the continuous development of computer technology, numerical simulation technology and numerical calculation technology are more and more widely used in the related research of aerospace vehicles. Compared with experimental methods such as wind tunnel tests, numerical simulation technology has the advantages of low experimental cost, strong repeatability, and strong applicability. Through the numerical simulation of aerospace vehicles, the flow characteristics, thermal characteristics and electromagnetic characteristics of aerospace vehicles can be obtained, and then the design of relevant parameters of aerospace vehicles can be optimized according to the results of numerical simulation. In the process of numerical simulation, various numerical calculation methods (such as finite element method, finite volume method, discontinuous Galerkin method, etc.) are widely used in aircraft design, and efficient and accurate meshing process is a key to these numerical calculation methods step.

Mesh division is the basis for high-precision numerical simulation and numerical calculation, and the quality of the grid directly affects the accuracy of numerical simulation and numerical calculation. And in the numerical simulation process, the meshing process often occupies most of the workload of the entire numerical simulation process, so meshing, especially the meshing of complex models, is a very important research part in numerical simulation. There are higher requirements for grid division in complex models and numerical simulations under complex conditions. Complex models usually contain parts with high curvature or local small geometric adjacent features. In order to ensure the discretization accuracy and grid quality, the linear tetrahedral grid needs to continuously reduce the grid size to match its geometric features in such geometric feature areas, which further leads to the continuous increase of the grid level, which leads to an increase in the amount of calculation and an increase in computation time. For complex shapes, poor mesh quality may also lead to divergence in the calculation process, which in turn leads to calculation failure. Therefore, it is necessary to improve the grid quality through some new technologies under complex conditions such as supersonic speed or multi-scale, and the curved grid is a better way to solve the existing grid quality problems.

The existing curved meshing methods basically include direct method and indirect method. The direct method is rarely used in numerical calculations due to its complex calculation and poor adaptability. The indirect method can be regarded as mesh deformation based on linear mesh division, and it is a commonly used method at present. In the indirect method, mesh deformation mainly includes extrapolation method, spring method and interpolation method. The extrapolation method is difficult to apply in complex models; the theory of the spring method is relatively complicated, and the calculation process takes a long time and has poor stability. The radial basis function method is an interpolation method that constructs an interpolation base through the distance between interpolation points. Because of its good versatility and good quality of deformed grids, it has been widely used in dynamic grids and grid deformation fields in recent years. However, the radial basis function method often requires a large number of calculations, and the solution process is very time-consuming and unstable.

Contents of the invention

In view of the above-mentioned problems or deficiencies, in order to solve the problems of poor grid division universality and poor grid quality in the existing complex model grid division process, the present invention provides a complex model three-dimensional space curved grid division method. By proposing a new control point and deformation point selection strategy, the selection of control points and deformation points is reduced while ensuring the grid quality, reducing a large number of position calculations in the meshing process, and optimizing the meshing strategy. It reduces the time-consuming mesh division and reduces the computing power demand.

The specific scheme is as follows, a complex model three-dimensional space curved meshing method, including the following steps:

Step 1. Input the target 3D physical model and grid parameters.

Step 2. Convert the target 3D physical model into 3D spatial geometric information, mark and store the calculation domain boundary information.

Step 3. According to the three-dimensional spatial geometric information in step 2, T initial linear tetrahedral elements are generated based on the front advancing method and the three-dimensional constrained Delaunay method in sequence, and at the same time, the grid topological relationship of the T initial linear tetrahedral elements is constructed The hierarchical structure of tetra (volume) → face (face) → edge (edge) → node (point) stores grid information.

Step 4. According to the grid information of the initial linear tetrahedral unit generated in step 3 and the computational domain boundary information in step 2, retrieve the T initial linear tetrahedral units generated in step 3 one by one.

If the tth unit is a boundary unit, attach a boundary mark to the unit and count it as the ѱth boundary unit, t∈{1,2,3,…, T }. By traversing the T initial linear tetrahedral elements in step 3, a total of Ψ boundary elements can be obtained, ѱ∈{1,2,3,…,Ψ}.

Step 5. According to the Ψ boundary units marked in step 4, retrieve the Ψ boundary units one by one, and take out the boundary surface, boundary edge and boundary point information of the ѱth boundary unit.

The boundary surface in the tetrahedral unit contains three boundary edges and three vertices, and the edges on all the boundary surfaces in the Ψ boundary unit are cut off to generate boundary surface truncation points. And project the generated boundary surface truncation point to the calculation domain boundary in step 2 to obtain the projection point, and there is a one-to-one correspondence between the projection point and the boundary surface truncation point.

Step 6. According to the Ψ boundary cell information obtained in step 5, retrieve the surface, edge and point information of the non-boundary surfaces of Ψ boundary cells one by one.

Step 7. According to the information of Ψ boundary units obtained in step 5, the adjacent units of Ψ boundary units are retrieved one by one, and the information of the faces, edges and points of the non-adjacent surfaces in the adjacent units is taken out.

The adjacent units: two units sharing one surface in the tetrahedral unit are defined as adjacent units, there are four adjacent units for non-boundary units, and there are less than four adjacent units for boundary units but at least one.

Step 8. Traverse all boundary surface truncation points and their corresponding projection points, and calculate the cosine of the angle between the boundary surface truncation points and their corresponding projection points.

Step 9. Substitute all the truncation points obtained by truncation in steps 5, 6 and 7 into the radial basis function calculation formula, and the value obtained through the radial basis function calculation formula is the deformation of the truncation point.

The control point is the truncation point of the boundary surface in step 5, that is, a subset of all boundary surface truncation points; first calculate the weight through the control point and the projection point; then take the global numbering sequence of the outer layer of the boundary surface first, and then take the inner layer of the boundary edge The sequence of the local numbers of the cycle traverses the truncation point of the boundary surface to calculate the displacement of the control point, and skips the repeated point.

Step 10. Substitute the deformation calculated in step 9 into the truncation point coordinates to calculate the space curvature grid coordinates (x _c , y _c , z _c ).

Step 11. The spatial curve grid points calculated in step 10 are used as new grid nodes, inserted into the grid data, and the grid topological relationship is updated. Then output the new spatial curve grid data in the corresponding grid data file format, and import the output grid data file into visualization.

Further, the specific process of step 3 is as follows: firstly, adopting the front advancing method, taking the surface of the calculation domain as the initial front to advance to the inside to generate a hierarchical grid of W layers, where W≥1, and W can be selected by the user according to the actual situation (W The larger the space curvature is, the higher the accuracy of the grid is, but the greater the computational complexity, the value of W is W∈{1,2,3,...,10}).

Then, in the remaining calculation domain, the initial linear tetrahedral element is generated by using the three-dimensional constrained Delaunay method with the current hierarchical grid surface as the initial constraint condition.

Further, when calculating the position of the truncation point in step 5, the scale parameter λ is set to 1, which is the midpoint of the edge.

Further, if a higher-precision curved mesh is to be constructed in step 5, new truncation points are continuously generated on the basis of the first truncation of each boundary unit. Then project the obtained new first-order truncation points to the boundary of the computational domain to obtain corresponding projection points and store them in the truncation point data class and the projection point data class, so as to obtain more higher-order truncation points and higher-order projection points in a cycle.

Further, the basis function of the radial basis function in step 9 is a compactly supported radial basis function.

Further, the radial basis function in step 9 uses a quadratic polynomial.

Further, the grid data file format output in the step 11 adopts the CGNS format.

Further, in step 6, the same method as in step 5 is used to obtain the truncation point of the non-boundary surface.

Further, the step 7 adopts the same method in the step 5 to truncate each side to obtain the truncation point of the adjacent unit.

Further, in the step 7, for the multi-layer spatial curvature grid, the adjacent units of the adjacent units are further cyclically retrieved.

The complex model three-dimensional space curved grid division method proposed by the present invention is realized by truncating the original linear tetrahedron to generate a curved tetrahedron unit, which can well meet the grid division requirements of the complex model without increasing the grid level , the initial linear grid is based on the front-advancing method and the three-dimensional constrained Delaunay method. In the process of generating space curve grid points, only the truncation point of the boundary surface in step 5 is selected as the control point, and only the truncation point of the boundary unit and the adjacent unit of the boundary unit is selected as the deformation point, which greatly reduces the number of control points and deformation points, thus The calculation related to the deformation point, that is, the grid point of the space curve is greatly reduced. Compared with the existing technology, step 8 innovatively proposes a basis function selection method based on the cosine value of the included angle; on the premise of ensuring the quality of the grid, the use of simple low-degree polynomial basis functions can greatly reduce the amount of deformation. required calculations. The calculation reduction of deformation points and the simplification of basis function form greatly improve the stability and speed of grid division.

In summary, the present invention reduces the selection of control points and deformation points while ensuring grid quality by proposing a new strategy for selecting control points and deformation points, and reduces a large number of position calculations in the grid division process. The grid division strategy is optimized to reduce the time-consuming grid division and the computing power demand. It satisfies the grid division requirements of complex models well without increasing the grid level.

Description of drawings

Fig. 1 is a schematic diagram of a three-dimensional space curved grid of an embodiment target model;

Figure 2 is a comparison diagram of blunt cone (10°) linear/curved grid;

Figure 3 is a schematic diagram of the cross-section (two-dimensional) of adjacent units;

Figure 4 is a schematic diagram of the X-51A linear grid;

Figure 5 is a schematic diagram of the X-51A spatial curvature grid;

Fig. 6 is a flowchart of the present invention.

Detailed ways

The technical solution of the present invention will be described in detail below in conjunction with the drawings and embodiments.

With reference to the accompanying drawings, a complex model three-dimensional curved mesh division method, comprising the following steps:

step 1. Input the target 3D physical model and grid parameters.

Draw a digital model based on the ACIS kernel according to the 3D physical model (dimensional parameters, formulas, drawings, objects, etc.), and read the grid parameter settings such as grid size and boundary type through the program.

Step 2. Convert the target 3D physical model into 3D spatial geometric information, mark and store the calculation domain boundary information.

Establish the Model class, import the three-dimensional space geometric information in the digital model into the model Model class, and establish the boundary mark according to the set boundary type (the member BCtype of the Model class in the program).

Step 3. According to the three-dimensional spatial geometric information in step 2, T initial linear tetrahedral elements are generated sequentially based on the front advancing method and the three-dimensional constrained Delaunay method, and a tetra body is constructed according to the grid topology of the T initial linear tetrahedral elements →face→edge→node hierarchy stores grid information.

To read the grid parameters and the three-dimensional spatial geometric information in the Model class, first use the front advance method to advance from the surface of the calculation area to the interior to generate multi-layer (two-layer in the spherical embodiment) grid units, and then use the three-dimensional constraint Deloitte The inner (Delaunay) method uses the current front as the constraint condition to generate the grid of the remaining area. At the same time, class tetraArray, faceArray, edgeArray and nodeArray are constructed to store unit information.

Step 4. According to the grid information of the initial linear tetrahedral unit generated in step 3 and the computational domain boundary information in step 2, retrieve the T initial linear tetrahedral units generated in step 3 one by one.

If the tth unit is a boundary unit, attach a boundary mark to the unit and count it as the ѱth boundary unit, t∈{1,2,3,…, T }; traverse the T initial linear tetrahedral units in step 3 , a total of Ψ boundary elements can be obtained, ѱ∈{1,2,3,…,Ψ}.

In this embodiment: the boundary marker BCtype is inserted into the tetraArray class, and a corresponding relationship is established. Then insert the BCtype in the tetraArray class (the member BCtype of the tetraArray class in the program) into the faceArray class, edgeArray class, and nodeArray class to mark the boundary surface, boundary edge and boundary point in the unit.

Step 5. According to the Ψ boundary units marked in step 4, retrieve the Ψ boundary units one by one, and take out the boundary surface, boundary edge and boundary point information of the ѱth boundary unit;

A boundary surface in a tetrahedral element consists of three boundary edges and vertices P ₁ (x ₁ ,y ₁ ,z ₁ ), P ₂ (x ₂ ,y ₂ ,z ₂ ) and P ₃ (x ₃ ,y ₃ ,z ₃ ), x, y and z respectively represent the coordinate values of points in three dimensions in the three-dimensional Cartesian coordinate system; the first side corresponds to vertices P ₁ and P ₂ , the second side corresponds to P ₁ and P ₃ , and the third side The edges correspond to P ₂ and P ₃ . Cut off the edges on all the boundary faces of Ψ boundary cells (tetrahedral cell faces contain three edges), and generate boundary face truncation points.

The position of the interception point of the first boundary edge (the second and third boundary edges are the same) is calculated by the vertices P ₁ (x ₁ ,y ₁ ,z ₁ ) and P ₂ (x ₂ ,y ₂ ,z ₂ ) , whose calculation formula is:

And project the generated boundary surface truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) ( H represents the label of the boundary surface truncation point) to the calculation domain boundary in step 2 to obtain the projected point P _sub ( x _n,1 ,y _n,1 ,z _n,1 ) ( sub represents the projection point label, and there is a one-to-one correspondence between the projection point and the boundary surface truncation point), λ is the proportional parameter (in most cases, the calculation step is omitted You can directly take λ=1 as the midpoint of the edge. In special complex models such as high curvature, it can be adjusted based on the curvature of the physical model and the radial basis function. The adjustment range is λ∈[0.5,2].

Every time a truncation point is generated, the truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) is stored in the truncation point data class P _H (x _N,1 ,y _N,1 ,z _{N, 1} ) and count (n in PH represents the nth truncation point, and _the maximum count value is N ). At the same time, all the projection points P _sub (x _n,1 ,y _n,1 ,z _n,1 ) obtained by projecting onto the boundary surface are also stored in the projection point data class P _sub (x _N,1 ,y _N,1 ,z _{N ,1} ) and count (n in P _sub represents the nth projection point, and the maximum count value is also N ).

Furthermore, if a higher-precision curved mesh is to be constructed, new truncation points are continuously generated on the basis of the first truncation of each boundary unit. Taking points P _H (x _n,1 ,y _n,1 ,z _n,1 ) and P ₁ (x ₁ ,y ₁ ,z ₁ ), and P _H (x _n,1 ,y _n,1 ,z _{n ,1} ) and P ₂ (x ₂ ,y ₂ ,z ₂ ) are endpoints respectively and continue according to

Generate a truncation point, at this time ( _xi , y _i , _zi ) is a truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ), (x _j ,y _j ,z _j ) is are the values of vertices P ₁ (x ₁ , y ₁ , z ₁ ) and P ₂ (x ₂ , y ₂ , z ₂ ), respectively.

Then project the obtained truncation points PH (x _n,21 , _{y n ,21} ,z _n,21 ) and PH (x _n,22 ,y _n,22 ,z _n,22 ) to the computational domain boundary Get P _sub (x _n,21 ,y _n,21 ,z _n,21 ) and P _sub (x _n,22 ,y _n,22 ,z _n,22 ) and store them in the corresponding truncation point data class and projection point data class. In this way, more high-order truncation points and high-order projection points are obtained. Generally, each edge is truncated once to obtain a first-order truncation point or twice to obtain three second-order truncation points, which can meet the grid accuracy requirements. Although more truncation points can be cyclically intercepted to improve grid accuracy, the more truncation times, the greater the amount of calculation, the slower the calculation speed and the greater the demand for computing power. Therefore, calculation schemes higher than the second-order cut-off point are not recommended in meshing.

In this embodiment, the boundary mark BCtype is searched in the tetraArray class, and when the unit is a boundary unit, the edgeArray class and nodeArray class corresponding to the unit are searched. Then further retrieve the boundary mark BCtype of the edge. When the boundary edge is retrieved, the truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) and the projection point P are calculated by the point coordinates corresponding to the boundary edge and the formula _sub (x _n,1 ,y _n,1 ,z _n,1 ). Obtain the truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) and projection point P _sub (x _n,1 ,y _n,1 ,z _n,1 ) and store them in the truncation point data class P _H (x _{N ,1} ,y _{N ,1} ,z _{N ,1} ) and P _sub (x _{N ,1} ,y _{N ,1} ,z _{N ,1} ), and count (the maximum count value of n is N ).

Step 6. According to the Ψ boundary cell information obtained in step 5, retrieve the surface, edge and point information of the non-boundary surface of Ψ boundary cells one by one; use the same method in step 5 to obtain the non-boundary surface truncation point P _inr (x _{m ,1} ,y _m,1 ,z _m,1 ), stored in the truncation point data class P _inr (x _M,1 ,y _M,1 ,z _M,1 ), and counted, inr represents the label of the non-boundary surface truncation point , the maximum count value of m in P _inr is denoted as M.

Step 7. According to the information of Ψ boundary units obtained in step 5, retrieve the adjacent units of Ψ boundary units one by one, and take out the information of the faces, edges and points of the non-adjacent surfaces in the adjacent units. If a multi-layer spatial curved grid is required, then A further loop retrieves neighbors of neighbors.

The adjacent units: two units sharing one surface in the tetrahedral unit are defined as adjacent units, there are four adjacent units for non-boundary units, and there are less than four adjacent units for boundary units but at least one.

Use the same method in step 5 to truncate each edge to obtain the truncation point P _adj (x _k,1 ,y _k,1 ,z _k,1 ) of the adjacent unit and store it in the truncation point data class P _adj (x _{K, 1} ,y _K,1 ,z _K,1 ), and counting, the maximum count value of k in P _adj is K, and adj represents the label of the truncation point of the adjacent unit.

In this embodiment, the class member tetra[A,B] of the faceArray class is retrieved, and A and B are the global numbers of the two units to which the face belongs. Then, according to the global number A or B, retrieve the adjacent units in the tetraArray class, and take out the information of the edges and points of the non-adjacent faces.

Step 8. Traversing the boundary surface truncation point data class P _H (x _N,1 ,y _N,1 ,z _N,1 ) and projection point data class P _sub (x _N,1 ,y _N,1 ,z _N,1 ) point, calculate the boundary surface truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) and its corresponding projection point P _sub (x _n,1 ,y _n,1 ,z _n,1 ) The cosine value of the included angle, denoted as cos( θ ) _H-sub , is stored in the corresponding position of each point in the nodeArray class.

Among them, x _{n, H} , y _{n, H} , z _{n, H} are the three-dimensional coordinates of the truncation point of the nth boundary surface, which means x _{n, sub} , y _{n, sub} , z _{n, sub} are the nth projection point 3D coordinate value.

Step 9. Substitute all the truncation points obtained by truncation in steps 5, 6 and 7 into the radial basis function calculation formula, and the value of f ( x, y, z ) obtained through the radial basis function calculation formula is the truncation point deformation; save Enter the temporary variable mesh-distancecurved (total N+M+K ). The calculation formula of radial basis function is:

为权重，I _N 为控制点总数且I _N =N；/>

为径向基函数，

表示全部截断点到对应控制点的欧式距离，控制点/>

为步骤5中边界面截断点P _H （x_N,1,y_N,1,z_N,1），即全部边界面截断点的子集；where α is the coefficient,

is the weight, I _N is the total number of control points and I _N = N ;/>

is the radial basis function,

Indicates the Euclidean distance from all truncation points to the corresponding control point, control point />

is the boundary surface truncation point P _H (x _N,1 ,y _N,1 ,z _N,1 ) in step 5, which is a subset of all boundary surface truncation points;

In the embodiment, based on the cosine value of the maximum included angle of the truncation point (that is, the spatial curvature grid point), the calculation formula of the radial basis function used in the present invention is:

Radial basis functions only use quadratic polynomials in most cases, and the use of simple low-degree polynomial basis functions can greatly reduce the calculation required for deformation; only in complex cases such as high curvature in complex models, high-degree polynomials are used. Mesh quality requirements are met.

，计算权重时/>

中X代表边界面截断点（即为步骤5中边界面截断点P _H （x_N,1,y_N,1,z_N,1），所以总数为I _N =N；计算权重时f (x,y,z)为已知量，代表控制点位移量，f (x,y,z)取同一边界面上截断点到投影点的欧式距离的最小值，以避免有平凡解的情况出现；First calculate the weight through the control point (the truncation point of the boundary surface) and the projection point

, when calculating weights />

where X represents the truncation point of the boundary surface (that is, the truncation point of the boundary surface P _H (x _N,1 ,y _N,1 ,z _N,1 ) in step 5, so the total number is I _N = N ; when calculating the weight f ( x , y, z ) are known quantities, representing the displacement of control points, f ( x, y, z ) take the minimum value of the Euclidean distance from the truncation point to the projected point on the same boundary surface to avoid the occurrence of trivial solutions;

f ( x,y,z )=min{║ X _b - X _sub ║}

In the case of the first-order truncation point, X _b , b∈{1,2,3} are the three truncation points corresponding to the three boundary edges on the same boundary surface, and X _sub is the coordinate value of the projection point corresponding to the truncation point; then Traverse the truncation points of the boundary surface to calculate the displacement of the control points according to the order of the global numbering of the outer loop of the boundary surface first, and then the local numbers of the inner loop of the boundary surface, and skip the repeated points;

后代入截断点计算f (x,y,z)，此时/>

为已知量，f (x,y,z)为待求量；计算网格点位移时/>

中X=X _j ，X _j 代表步骤5中边界面截断点、步骤6中非边界面截断点和步骤7中相邻单元截断点求得的全部截断点，截断点总数I _SUM =N+M+K；get the weight

Substituting into the truncation point to calculate f ( x,y,z ), at this time />

is the known quantity, f ( x,y,z ) is the quantity to be sought; when calculating the grid point displacement />

Where X = X _j , X _j represents all the truncation points obtained from the boundary surface truncation point in step 5, the non-boundary surface truncation point in step 6 and the adjacent unit truncation point in step 7, the total number of truncation points I _SUM = N+ M + K;

Step 10. Substitute the deformation variable f ( x,y,z ) calculated in step 9 into the truncation point coordinates (x _J,1 ,y _J,1 ,z _J,1 ) to calculate the spatial curvature grid coordinates (x _c ,y _c , z _c ), its calculation formula is:

(x _c , y _c , z _c ) = (x _J,1 ,y _J,1 ,z _J,1 ) + f ( x,y,z ),J∈{ N } , {M} , {K}

Traverse the temporary variable mesh-distancecurved to take out the deformation value corresponding to each truncation point, and add it to the initial position of the truncation point to obtain the position of the truncation point after deformation (the position of the truncation point after addition is the space curve grid point), which is Space curve grid point coordinates. Then store the coordinates of the obtained space curve grid points into the truncation point class to replace the original truncation point coordinates.

Step 11. The high-order grid points (x _c , y _c , z _c ) calculated in step 10 are used as new grid nodes, that is, the space curve grid point coordinates calculated in step 10 are read from the truncation point class, and used as new grid nodes Grid nodes node5 to node10 are inserted into the grid data, and the grid topological relationship is updated; then the new spatial curvature grid data is output in the corresponding grid data file format (such as CGNS format), and the output grid data File import visualization.

In this embodiment, three target objects are used to make an example to verify the effectiveness of the present invention; the space curvature grid results of the final grid division of each embodiment are shown in Figure 1 and the right figure of Figure 2, and Figure 5.

Figure 1 takes a sphere as the target object. The left picture is the global surface view of the space curved grid, the middle picture is the section view of the space curved grid, and the right picture is the partial enlarged view.

Figure 2 is the results of blunt cone (10°) space curvature mesh. The left picture is the initial linear tetrahedral mesh obtained in step 3, and the right picture is the schematic diagram of the space curved mesh after completion.

Figures 4 and 5 are resulting plots for the X-51A aircraft. Figure 4 is the initial linear tetrahedral grid obtained in step 3, and Figure 5 is a schematic diagram of the completed space curved grid.

From the comparison of the left and right diagrams in Fig. 2, and the comparison of Fig. 4 and Fig. 5, it can be seen that, under the condition that the grid quality and the basic topological structure remain unchanged, the space curved mesh division method proposed by the present invention realizes High-precision space curved grid. The existing curved mesh technology basically only realizes the curved mesh of the object surface. Although the curved mesh of the object surface improves the calculation accuracy at the boundary, the numerical calculation accuracy is often limited by the accuracy of the internal mesh. The present invention adopts a brand-new selection strategy of control points and deformation points (that is, space curved grid points), and extends the curved grid to the space inner layer grid to realize high-precision spatial curved grid division of complex models. In summary, it can be seen that the present invention is practical and effective.

It can be seen from the above that, by proposing a new strategy for selecting control points and deformation points, the present invention satisfies the grid division requirements of complex models well without increasing the grid level. Compared with the existing technology: 1. The present invention realizes space curved grid division, and the existing technology is basically a curved surface grid; 2. Adopting a new control point and deformation point selection strategy and calculation method, without upgrading the grid In the case of an order of magnitude, it can well meet the grid division requirements of complex models; 3. The radial basis function is selected according to the cosine value of the included angle. The present invention reduces the selection of control points and deformation points while ensuring the grid quality, reduces a large number of position calculations in the grid division process, optimizes the grid division strategy, reduces the time-consuming grid division and reduces the computing power need. The calculation reduction of deformation points and the simplification of basis function form greatly improve the stability and speed of grid division.

Claims (5)

1. A complex model three-dimensional space curved meshing method, is characterized in that, comprises the following steps: Step 1. Input the target 3D physical model and grid parameters; Step 2. Convert the target 3D physical model into 3D spatial geometric information, mark and store the computational domain boundary information; Step 3. According to the three-dimensional spatial geometric information in step 2, T initial linear tetrahedral elements are generated based on the front advancing method and the three-dimensional constrained Delaunay method in turn, and at the same time, according to the grid topology relationship of the T initial linear tetrahedral elements Construct a hierarchical structure of tetra body → face → edge → node to store grid information; Step 4. According to the grid information of the initial linear tetrahedral unit generated in step 3 and the computational domain boundary information in step 2, retrieve the T initial linear tetrahedral units generated in step 3 one by one; If the tth unit is a boundary unit, attach a boundary mark to the unit and count it as the ψth boundary unit, t∈{1,2,3,…,T}; traverse the T initial linear tetrahedral units in step 3 , a total of Ψ boundary elements can be obtained, ψ∈{1,2,3,…,Ψ}; Step 5. According to the Ψ boundary units marked in step 4, retrieve the Ψ boundary units one by one, and take out the boundary surface, boundary edge and boundary point information of the ψth boundary unit; The boundary surface contains three boundary edges and vertices P ₁ (x ₁ ,y ₁ ,z ₁ ), P ₂ (x ₂ ,y ₂ ,z ₂ ) and P ₃ (x ₃ ,y ₃ ,z ₃ ), x, y and z respectively represent the coordinate values of points in three dimensions in the three-dimensional Cartesian coordinate system; the first side corresponds to vertices P ₁ and P ₂ , the second side corresponds to P ₁ and P ₃ , and the third side corresponds to P ₂ and P ₃ ; the edges on all boundary surfaces in the Ψ boundary cells are truncated to generate boundary surface truncation points, and the tetrahedron element faces include three edges; The position of the truncation point of the first boundary edge is calculated by the vertices P ₁ (x ₁ ,y ₁ ,z ₁ ) and P ₂ (x ₂ ,y ₂ ,z ₂ ), and the same is true for the second and third boundary edges. Its calculation formula is:

And project the generated boundary surface truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) to the calculation domain boundary in step 2 to obtain the projected point P _sub (x _n,1 ,y _n,1 ,z _n,1 ), where H represents the label of the truncation point of the boundary surface, sub represents the label of the projection point, and there is a one-to-one correspondence between the projection point and the truncation point of the boundary surface; λ is a proportional parameter, based on the physical model curvature, radial basis function to adjust, the adjustment range is λ∈[0.5,2]; Every time a truncation point is generated, the truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) is stored in the truncation point data class P _H (x _N,1 ,y _N,1 ,z _{N, 1} ), and counting, n in P _H represents the nth truncation point, and the maximum count value is N; at the same time, project all projections to the boundary surface to obtain the projection point P _sub (x _n,1 ,y _n,1 ,z _{n, 1} ) Also stored in the projection point data class P _sub (x _N,1 ,y _N,1 ,z _N,1 ), and counted, n in P _sub represents the nth projection point, and the maximum count value is also N; Step 6. According to the Ψ boundary cell information obtained in step 5, retrieve the surface, edge and point information of the non-boundary surface of Ψ boundary cells one by one; use the same method in step 5 to obtain the non-boundary surface truncation point P _inr (x _m,1 ,y _m,1 ,z _m,1 ), stored in the truncation point data class P _inr (x _M,1 ,y _M,1 ,z _M,1 ), and counted, inr represents the non-boundary surface The label of the truncation point, the maximum count value of m in P _inr is recorded as M; Step 7. According to the Ψ boundary unit information obtained in step 5, retrieve the adjacent units of Ψ boundary units one by one, and take out the information of the faces, edges and points of the non-adjacent surfaces in the adjacent units. If multi-layer space curve is required The grid further loops to retrieve neighbors of neighbors; The adjacent units: two units sharing one face in the tetrahedron unit are defined as adjacent units, there are four adjacent units of non-boundary units, and there are less than four adjacent units of boundary units but at least one; Use the same method in step 5 to truncate each side to obtain the truncation point P _adj (x _k,1 ,y _k,1 ,z _k,1 ) of the adjacent unit and store it in the truncation point data class P _adj (x _K,1 , y _K,1 ,z _K,1 ), and counting, the maximum counting value of k in P _adj is K, and adj represents the label of the truncation point of the adjacent unit; Step 8. Traverse the boundary surface truncation point data class P _H (x _N,1 ,y _N,1 ,z _N,1 ) and projection point data class P _sub (x _N,1 ,y _N,1 ,z _N,1 ), calculate the boundary surface truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ) and its corresponding projection point P _sub (x _n,1 ,y _n,1 ,z _n,1 ) The cosine value of the angle between is recorded as cos(θ) _H-sub ;

Among them, x _{n, H} , y _{n, H} , z _{n, H} are the three-dimensional coordinates of the truncation point of the nth boundary surface, which means x _{n, sub} , y _{n, sub} , z _{n, sub} are the nth projection point 3D coordinate value; Step 9. Substitute all the truncation points obtained by truncation in steps 5, 6 and 7 into the radial basis function calculation formula, and the value of f(x, y, z) obtained by the radial basis function calculation formula is the truncation point shape variable; The calculation formula of radial basis function is:

where α is the coefficient,

is the weight, I _N is the total number of control points and I _N = N; />

is the radial basis function, />

Indicates the Euclidean distance from all truncation points to the corresponding control point, control point />

is the boundary surface truncation point P _H (x _N,1 ,y _N,1 ,z _N,1 ) in step 5, which is a subset of all boundary surface truncation points; First calculate the weight through the control point and projection point

When calculating weights />

X represents the truncation point of the boundary surface, so the total number is I _N = N; f(x,y,z) is a known quantity when calculating the weight, which represents the displacement of the control point, and f(x,y,z) takes the same boundary surface The minimum value of the Euclidean distance from the upper truncation point to the projected point; f(x,y,z)=min{║X _b -X _sub ║} In the case of the first-order truncation point, X _b , b∈{1,2,3} are the three truncation points corresponding to the three boundary edges on the same boundary surface, and X _sub is the coordinate value of the projection point corresponding to the truncation point; press First take the global numbering sequence of the outer loop of the boundary surface, then take the sequence of the local numbers of the inner loop of the boundary edge, traverse the truncation point of the boundary surface to calculate the displacement of the control point, and skip the repeated point; get the weight

Substituting into the truncation point to calculate f(x,y,z), at this time />

is the known quantity, f(x,y,z) is the quantity to be sought; when calculating the grid point displacement />

Among them, X=X _j , X _j represents all the truncation points obtained from the boundary surface truncation point in step 5, the non-boundary surface truncation point in step 6 and the adjacent unit truncation point in step 7, the total number of truncation points I _SUM =N+M +K; Step 10. Substitute the deformation variable f(x,y,z) calculated in step 9 into the coordinates of the truncation point (x _J,1 ,y _J,1 ,z _J,1 ) to calculate the spatial curvature grid coordinates (x _c , y _c , z _c ), the calculation formula is: (x _c ,y _c ,z _c )=(x _J,1 ,y _J,1 ,z _J,1 )+f(x,y,z),J∈{N},{M},{K} Step 11. The high-order grid points (x _c , y _c , z _c ) calculated in step 10 are used as new grid nodes, inserted into the grid data, and the grid topological relationship is updated; then the new spatial curved network The grid data is output in the corresponding grid data file format, and the output grid data file is imported into visualization. 2. complex model three-dimensional space curved grid division method as claimed in claim 1, is characterized in that, in described step 3: Firstly, the frontal advancing method is used to generate a hierarchical grid of W layers with the surface of the computational domain as the initial frontal surface, where W≥1, W∈{1,2,3,…,10}; Then, in the remaining computational domain, the initial linear tetrahedral element is generated by using the three-dimensional constrained Delaunay method with the current hierarchical grid surface as the initial constraint condition. 3. The complex model three-dimensional space curved mesh division method according to claim 1, characterized in that: in the step 5, λ=1, which is the midpoint of the side. 4. complex model three-dimensional space curved grid division method as claimed in claim 1, is characterized in that: In said step 5: If a high-precision curved mesh is to be constructed, new truncation points shall be continuously generated on the basis of the first truncation of each boundary element; with point P _H (x _n,1 ,y _n,1 ,z _n,1 ) and P ₁ (x ₁ ,y ₁ ,z ₁ ), and P _H (x _n,1 ,y _n,1 ,z _n,1 ) and P ₂ (x ₂ ,y ₂ ,z ₂ ) are endpoints respectively continue according to

Generate a truncation point, at this time ( _xi , y _i , z _i ) is a primary truncation point P _H (x _n,1 ,y _n,1 ,z _n,1 ), (x _j ,y _j ,z _j ) respectively is the value of vertices P ₁ (x ₁ ,y ₁ ,z ₁ ) and P ₂ (x ₂ ,y ₂ ,z ₂ ); Then project the obtained truncation points P _H (x _n,21 ,y _n,21 ,z _n,21 ) and P _H (x _n,22 ,y _n,22 ,z _n,22 ) to the computational domain boundary Get P _sub (x _n,21 ,y _n,21 ,z _n,21 ) and P _sub (x _n,22 ,y _n,22 ,z _n,22 ) and store them in the corresponding truncation point data class and projection point The data class, so loop to get more high-order truncation points and high-order projection points. 5. The complex model three-dimensional space curved grid division method according to claim 1, characterized in that: the grid data file format output in the step 11 adopts the CGNS format.

Images