CN114510775B - Method for dividing three-dimensional space curved grids of complex model - Google Patents
Method for dividing three-dimensional space curved grids of complex model Download PDFInfo
- Publication number
- CN114510775B CN114510775B CN202111637505.5A CN202111637505A CN114510775B CN 114510775 B CN114510775 B CN 114510775B CN 202111637505 A CN202111637505 A CN 202111637505A CN 114510775 B CN114510775 B CN 114510775B
- Authority
- CN
- China
- Prior art keywords
- boundary
- point
- points
- cut
- grid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Aviation & Aerospace Engineering (AREA)
- Automation & Control Theory (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Image Generation (AREA)
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
Technical Field
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.
Background
With the continuous development of computer technology, numerical simulation technology and numerical calculation technology are increasingly widely applied to related researches of aerospace vehicles. Compared with experimental methods such as wind tunnel test, the numerical simulation technology has the advantages of low experimental cost, strong repeatability, strong applicability and the like. Through numerical simulation of the aerospace vehicle, the flow characteristics, the thermal characteristics and the electromagnetic characteristics of the aerospace vehicle can be obtained, and then the design of relevant parameters of the aerospace vehicle is optimized according to the result of the numerical simulation. Various numerical calculation methods (such as finite element method, finite volume method, intermittent Galerkin method and the like) are widely applied to aircraft design in the numerical simulation process, and an efficient and accurate grid division process is a key step of the numerical calculation methods.
The grid division is the basis for carrying out high-precision numerical simulation and numerical calculation, and the quality of the grid directly influences the precision of the numerical simulation and the numerical calculation. In addition, in the numerical simulation process, the meshing process often occupies most of the workload of the whole numerical simulation process, so meshing, especially meshing of complex models, is an important research part in numerical simulation. There is a higher demand for meshing in complex models and numerical simulations under complex conditions. Complex models typically contain highly curved portions or locally tiny geometrically contiguous features. To ensure discrete accuracy and mesh quality, linear tetrahedral meshes in such geometric feature regions require ever shrinking mesh dimensions to match their geometric features, further leading to ever increasing mesh magnitudes, and thus to increased computational effort and computation time. For complex shapes, poor mesh quality can also lead to divergence in the computation process and thus to computation failure. Therefore, it is necessary to improve the quality of the grid by some new techniques under the complex conditions of supersonic speed or multi-scale, etc., and curved grid is a better way to solve the existing grid quality problem.
The existing curved meshing method basically comprises a direct method and an indirect method. The direct method uses less in numerical calculation because of the problems of complex calculation, poor adaptability, and the like. The indirect method can be regarded as grid deformation based on linear grid division, and is a method commonly used at present. In the indirect method, the mesh deformation mainly includes extrapolation, spring method, interpolation method, and the like. Extrapolation is difficult to apply in complex models; the spring theory is complex, and the calculation process is long and has poor stability. The radial basis function method is an interpolation method for constructing interpolation bases through the distances between interpolation points, and is widely applied to the fields of dynamic grids, grid deformation and the like in recent years due to good universality and good quality of deformed grids. However, radial basis function methods often require a large number of computations, and the solution process is time consuming and unstable.
Disclosure of Invention
Aiming at the problems or the defects, the invention provides a three-dimensional space curve meshing method of a complex model, which aims at solving the problems of poor universality and poor quality of meshing in the meshing process of the existing complex model. By providing a new control point and deformation point selection strategy, the selection of the control point and the deformation point is reduced while the grid quality is ensured, a large number of position calculations in the grid division process are reduced, the grid division strategy is optimized, the grid division time consumption is reduced, and the calculation force requirement is reduced.
The method for meshing the three-dimensional space curve of the complex model comprises the following steps:
and step 1, inputting a target three-dimensional physical model and grid parameters.
And 2, converting the target three-dimensional physical model into three-dimensional space geometric information, marking and storing calculation domain boundary information.
Step 3, generating based on the array plane propulsion method and the three-dimensional constrained Delaunay method in sequence according to the three-dimensional space geometric information in the step 2TInitial linear tetrahedral units, simultaneously according toTThe grid topology of the initial linear tetrahedral units constructs a hierarchy of tetra (volume) →face →edge →node, storing grid information.
Step 4, searching the grid information of the initial linear tetrahedral unit generated in the step 3 one by one and calculating the domain boundary information in the step 2TInitial linear tetrahedral units.
If the t-th cell is a boundary cell, the cell is attached with a boundary tag and is counted as ѱ -th boundary cell, t e {1,2,3, …,T}. Traversing step 3TA total of ψ boundary units, ѱ e {1,2,3, …, ψ }, can be obtained from the initial linear tetrahedral units.
And 5, searching the psi boundary units one by one according to the psi boundary units marked in the step 4, and taking out the boundary surface, boundary edge and boundary point information of the ѱ th boundary unit.
The boundary surface in the tetrahedron unit comprises three boundary edges and three vertexes, and the edges on all boundary surfaces in the tetrahedron unit are cut off to generate boundary surface cut-off points. And projecting the generated boundary surface cut-off points to the calculated domain boundary in the step 2 to obtain projection points, wherein the projection points and the boundary surface cut-off points have a one-to-one correspondence.
And 6, searching the surface, edge and point information of the non-boundary surfaces of the psi boundary units one by one according to the psi boundary unit information obtained in the step 5.
And 7, retrieving adjacent units of the psi boundary units one by one according to the psi boundary unit information obtained in the step 5, and taking out the information of the faces, edges and points of the non-adjacent faces in the adjacent units.
The adjacent unit: two cells sharing a face in a tetrahedral cell are defined as adjacent cells, four adjacent cells being non-boundary cells, less than four but at least one adjacent cells being boundary cells.
And 8, traversing all boundary surface cut-off points and corresponding projection points, and calculating an included angle cosine value between the boundary surface cut-off points and the corresponding projection points.
And 9, substituting all the truncated points obtained by the truncations of the steps 5, 6 and 7 into a radial basis function calculation formula, wherein the value obtained by the radial basis function calculation formula is the truncated point deformation.
The control points are boundary surface cut-off points in the step 5, namely a subset of all boundary surface cut-off points; firstly, calculating weights through control points and projection points; and then, according to the global numbering sequence of the outer layer circulation of the boundary surface, and then, the sequence of the local numbering of the inner layer circulation of the boundary edge is taken, the control point displacement is calculated by traversing the cut-off points of the boundary surface, and the point skipping is repeated.
Step 10, substituting the deformation calculated in the step 9 into the cut-off point coordinates to calculate the space curved grid coordinates (x) c ,y c ,z c ）。
Step 11, the space curved grid points calculated in the step 10 are used as newly added grid nodes, inserted into grid data and updated in grid topological relation. And then outputting the new space curved grid data in a corresponding grid data file format, and importing the output grid data file into a visualization.
Further, the specific process of the step 3 is as follows: firstly, a hierarchical grid of W layers is generated by adopting an array surface pushing method to push the calculation domain surface as an initial array surface, wherein W is more than or equal to 1, and W can be selected by a user according to actual conditions (the larger W is, the higher the spatial curved grid accuracy is, but the larger the calculated amount is, the W value is W epsilon {1,2,3, …,10 }).
And then generating initial linear tetrahedral units in the residual calculation domain by adopting a three-dimensional constrained Delaunay method by taking the current hierarchical grid surface as an initial constraint condition.
Further, in the step 5, when calculating the position of the break point, the ratio parameter λ takes 1, which is the midpoint of the edge.
Further, in step 5, if a curved grid with higher precision is to be constructed, a new cutoff point is continuously generated on the basis of the first cutoff of each boundary unit. And then projecting the new first-order cut-off points obtained respectively to the calculated domain boundary to obtain corresponding projection points and storing the corresponding cut-off point data class and the corresponding projection point data class, so that more high-order cut-off points and more high-order projection points are obtained in a circulating way.
Further, the basis function in the radial basis function in the step 9 is a tightly-supported radial basis function.
Further, the radial basis function in step 9 uses a quadratic polynomial.
Further, in the step 11, the grid data file format output adopts CGNS format.
Further, the step 6 adopts the same method as that in the step 5 to obtain the non-boundary surface cutting point.
Furthermore, in the step 7, each edge is truncated by the same method as in the step 5 to obtain the adjacent unit truncated points.
Further, in the step 7, the adjacent units of the adjacent units are further retrieved circularly for the multi-layer spatial curved grid.
The method for meshing the three-dimensional space curved grid of the complex model is realized by generating curved surface tetrahedron units by cutting off original linear tetrahedrons, can well meet the meshing requirement of the complex model under the condition of not lifting the grid magnitude, and is based on an array surface propulsion method and a three-dimensional constraint Delaunay method. In the generation process of the space curved grid points, only the boundary surface cut-off points in the step 5 are selected as control points, and the deformation points are also only selected as cut-off points of boundary units and adjacent units of the boundary units, so that the number of the control points and the deformation points is greatly reduced, and the deformation points, namely the calculation related to the space curved grid points, is greatly reduced. Compared with the prior art, the step 8 creatively provides a basis function selection mode based on an included angle cosine value; the calculation required for the deformation quantity can be greatly reduced by using the low-order polynomial basis function with simple form on the premise of ensuring the quality of the grid. The stability and the speed of grid division are greatly improved due to the reduction of deformation point calculation and the simplification of a basis function form.
In summary, the invention provides a new control point and deformation point selection strategy, so that the selection of the control point and the deformation point is reduced while the grid quality is ensured, a large amount of position calculation in the grid division process is reduced, the grid division strategy is optimized, the grid division time consumption is reduced, and the calculation force requirement is reduced. The grid division requirement of the complex model can be well met under the condition that the grid magnitude is not lifted.
Drawings
FIG. 1 is a schematic diagram of a three-dimensional curved grid of an exemplary object model;
FIG. 2 is a graph of a blunt cone (10) linear/curved grid comparison;
FIG. 3 is a schematic view of adjacent cells in cross section (two dimensions);
FIG. 4 is a schematic view of an X-51A linear grid;
FIG. 5 is a schematic view of an X-51A spatial curved grid;
fig. 6 is a flow chart of the present invention.
Detailed Description
The technical scheme of the invention is described in detail below with reference to the accompanying drawings and examples.
Referring to the attached drawings, the three-dimensional curved grid dividing method of the complex model comprises the following steps:
and step 1, inputting a target three-dimensional physical model and grid parameters.
And drawing a digital model based on an ACIS kernel according to a three-dimensional physical model (size parameters, formulas, drawings, real objects and the like), and reading grid parameter settings such as grid size, boundary type and the like through a program.
And 2, converting the target three-dimensional physical model into three-dimensional space geometric information, marking and storing calculation domain boundary information.
And (3) establishing Model classes, importing the three-dimensional space geometric information in the digital Model into the Model classes, and establishing boundary marks (BCtype which is a member of the Model classes in the program) according to the set boundary types.
Step 3, generating based on the array plane propulsion method and the three-dimensional constraint Delaunay method in sequence according to the three-dimensional space geometric information in the step 2TInitial linear tetrahedral units, simultaneously according toTThe grid topological relation of each initial linear tetrahedron unit constructs a hierarchy structure of a tetrad body, a face surface, edge edges and node points to store grid information.
The grid parameters and the three-dimensional space geometric information in Model classes are read, firstly, a multi-layer (two-layer in the ball embodiment) grid unit is generated by pushing the surface of a calculation area inwards by adopting an array plane pushing method, and then a grid of the residual area is generated by adopting a three-dimensional constraint Delaunay method and taking the current array plane as a constraint condition. And simultaneously constructing a tetraArray class, a faceArray class, an edgeArray class and a nodeArray class to store the unit information.
Step 4, searching the grid information of the initial linear tetrahedral unit generated in the step 3 one by one and calculating the domain boundary information in the step 2TInitial linear tetrahedral units.
If the t-th cell is a boundary cell, the cell is attached with a boundary tag and is counted as ѱ -th boundary cell, t e {1,2,3, …,T-a }; traversing step 3TA total of ψ boundary units, ѱ e {1,2,3, …, ψ }, can be obtained from the initial linear tetrahedral units.
In this embodiment: and inserting the boundary mark BCtype into the tetraArray class, and establishing a corresponding relation. And inserting BCtype (BCtype which is a member of the tetraArray class in the program) in the tetraArray class into faceArray class, and marking boundary surfaces, boundary edges and boundary points in units in the faceArray class and the nodeArray class.
Step 5, searching the psi boundary units one by one according to the psi boundary units marked in the step 4, and taking out boundary surface, boundary edge and boundary point information of the ѱ th boundary unit;
the boundary surface in the tetrahedral unit comprises three boundary edges and verticesP 1 （x 1 ,y 1 ,z 1 ）、P 2 （x 2 ,y 2 ,z 2 ) AndP 3 （x 3 ,y 3 ,z 3 ) X, y and z respectively represent coordinate values of the point in three dimensions under a three-dimensional Cartesian rectangular coordinate system; the first edge corresponds to the vertexP 1 AndP 2 the second edge corresponds toP 1 AndP 3 the third side corresponds toP 2 AndP 3 . Edges (tetrahedral unit faces contain three edges) on all boundary faces in the psi boundary units are cut off, and boundary face cut-off points are generated.
The position of the point of the truncated position of the first boundary edge (the second boundary edge and the third boundary edge are the same) passes through the vertexP 1 （x 1 ,y 1 ,z 1 ) AndP 2 （x 2 ,y 2 ,z 2 ) The calculation formula is as follows:
and intercept the generated boundary surfaceP H （x n,1 ,y n,1 ,z n,1 ）（HRepresenting the index of the boundary surface cut-off point) to the calculated domain boundary in step 2 to obtain a projection pointP sub （x n,1 ,y n,1 ,z n,1 ）（subRepresenting the projection point reference number, and the projection point and the boundary surface intercept point have one-to-one correspondence), lambda is a proportional parameter (in most cases, lambda=1 can be directly taken for omitting the calculation stepThe specific complex models such as the middle point and the high curvature of the edge can be adjusted based on the curvature of the physical model and the radial basis function, and the adjustment range is lambda epsilon [0.5,2]。
Each time a break point is generated, the break point will be brokenP H （x n,1 ,y n,1 ,z n,1 ) Logging of cut-off point data classesP H （x N,1 ,y N,1 ,z N,1 ) Middle and count%P H N represents the nth cut-off point, and the maximum count value isN). Simultaneously projecting all the projections onto the boundary surface to obtain projection pointsP sub （x n,1 ,y n,1 ,z n,1 ) Also stored in the proxel data classP sub （x N,1 ,y N,1 ,z N,1 ) Middle and count%P sub N represents the nth projection point, and the maximum count value is alsoN）。
Further, if a higher precision curved grid is to be constructed, new cutoff points continue to be generated on the basis of the first cutoff of each boundary cell. In dotsP H （x n,1 ,y n,1 ,z n,1 ) AndP 1 （x 1 ,y 1 ,z 1 ) And (b)P H （x n,1 ,y n,1 ,z n,1 ) AndP 2 （x 2 ,y 2 ,z 2 ) Respectively, continue according to the end points
A cut-off point is generated, at this time (x i ,y i ,z i ) Is a one-time cut-off pointP H （x n,1 ,y n,1 ,z n,1 ），（x j ,y j ,z j ) Are respectively vertexesP 1 （x 1 ,y 1 ,z 1 ) AndP 2 （x 2 ,y 2 ,z 2 ) Is a value of (2).
Then will be separately foundThe resulting cut-off pointP H （x n,21 ,y n,21 ,z n,21 ) AndP H （x n,22 ,y n,22 ,z n,22 ) Projection onto computational domain boundaries to obtainP sub （x n,21 ,y n,21 ,z n,21 ) AndP sub （x n,22 ,y n,22 ,z n,22 ) And stores the corresponding cut-off point data class and projection point data class. More high-order cutoff points and high-order projection points are circularly obtained, and in general, each edge is cut once to obtain one first-order cutoff point or cut twice to obtain three second-order cutoff points, so that the grid precision requirement can be met. Although more cut points can be circularly cut to improve the grid precision, the more the cut times are, the larger the calculation amount is, the slower the calculation speed is and the larger the calculation force requirement is. Therefore, calculation schemes above the second order cut-off point are not recommended in meshing.
The embodiment is as follows: and searching the boundary mark BCtype in the tetraArray class, and searching the edgeArray class and the nodeArray class corresponding to the unit when the unit is the boundary unit. And further searching the boundary mark BCtype of the boundary, and solving a cutoff point according to the point coordinates and formulas corresponding to the boundary when the boundary is searchedP H （x n,1 ,y n,1 ,z n,1 ) And a proxelP sub （x n,1 ,y n,1 ,z n,1 ). Obtaining a cutting pointP H （x n,1 ,y n,1 ,z n,1 ) And a proxelP sub （x n,1 ,y n,1 ,z n,1 ) Post-logging of cut-off point data classesP H （x N,1 ,y N,1 ,z N,1 ) AndP sub （x N,1 ,y N,1 ,z N,1 ) And count (nmax count value isN）。
Step 6, searching the surface, edge and point information of the non-boundary surface of the psi boundary units one by one according to the psi boundary unit information obtained in the step 5; the same method as in the step 5 is adopted for solvingTo obtain non-boundary surface cut-off pointP inr （x m,1 ,y m,1 ,z m,1 ) Storing cut-off point data classP inr （x M,1 ,y M,1 ,z M,1 ) And the number of the counting steps is counted,inrrepresenting the non-boundary surface discontinuity labels,P inr the M maximum count value of M is denoted as M.
And 7, retrieving adjacent units of the psi boundary units one by one according to the psi boundary unit information obtained in the step 5, taking out the information of the faces, edges and points of the non-adjacent faces in the adjacent units, and further circularly retrieving the adjacent units of the adjacent units if multi-layer space curved grids are needed.
The adjacent unit: two cells sharing a face in a tetrahedral cell are defined as adjacent cells, four adjacent cells being non-boundary cells, less than four but at least one adjacent cells being boundary cells.
Cutting off each edge by the same method in step 5 to obtain adjacent unit cut-off pointsP adj （x k,1 ,y k,1 ,z k,1 ) Logging of cut-off point data classesP adj （x K,1 ,y K,1 ,z K,1 ) In the process, and counting the number of the groups,P adj the K maximum count value of (c) is K,adjrepresenting adjacent cell cut-off point labels.
In this embodiment, the class members tetra [ a, B ] of the faceArray class are retrieved, and a and B are the global numbers of two units to which the face belongs. And searching adjacent units in the tetraArray class according to the global number A or B, and taking out the information of the sides and the points of the non-adjacent surfaces.
Step 8, traversing boundary surface cut-off point data classP H （x N,1 ,y N,1 ,z N,1 ) And proxel data classP sub （x N,1 ,y N,1 ,z N,1 ) All points in (1) calculating boundary surface cut-off pointsP H （x n,1 ,y n,1 ,z n,1 ) Corresponding to the projection pointP sub （x n,1 ,y n,1 ,z n,1 ) Cosine value of included angle between themIs recorded as cos @θ) H-sub And storing the corresponding position of each point in the nodeArray class.
Wherein x is n,H ,y n,H ,z n,H Is the firstnThree-dimensional coordinate values of the cutting points of the boundary surfaces represent x subn, ,y subn, ,z subn, Is the firstnThree-dimensional coordinate values of the respective projection points.
Step 9, substituting all the truncated points obtained in the steps 5, 6 and 7 into the radial basis function calculation formula, and obtaining the product through the radial basis function calculation formulaf(x,y,z) The value of (2) is the break point deformation; deposit temporary variable mesh-distance measure (total number)N+M+K) Is a kind of medium. The radial basis function is calculated as:
wherein the method comprises the steps ofαAs the coefficient of the light-emitting diode,
as the weight of the material to be weighed,I N is the total number of control points andI N =N；/>
as a function of the radial basis function,
representing the Euclidean distance from all the cut-off points to the corresponding control point, control point +.>
For the boundary surface intercept point in step 5P H （x N,1 ,y N,1 ,z N,1 ) I.e. a subset of all boundary surface intercept points;
in an embodiment, based on the maximum included angle cosine value of the intercept point (i.e., the spatial curved grid point), the radial basis function calculation formula used in the present invention is:
the radial basis function uses only quadratic polynomials in most cases, and the calculation required by deformation quantity can be greatly reduced by using the low-order polynomial basis function with simple form; and only when complex conditions such as high curvature and the like occur in a complex model, a high-order polynomial is adopted to meet the grid quality requirement.
Firstly, calculating weight by using control point (boundary surface cut-off point) and projection point
Weight is calculated +.>
In (a)XRepresents the boundary surface cut-off point (i.e., the boundary surface cut-off point in step 5)P H （x N,1 ,y N,1 ,z N,1 ) So the total isI N =NThe method comprises the steps of carrying out a first treatment on the surface of the When calculating the weightf (x,y,z) Is a known quantity, representing the displacement of the control point,f (x,y,z) Taking the minimum value of the Euclidean distance from the cutoff point to the projection point on the same boundary surface to avoid the occurrence of the condition of trivial solution;
f (x,y,z)=min{║X b -X sub ║}
in the case of a first-order cut-off pointX b B epsilon {1,2,3} is three cut-off points corresponding to three boundary edges on the same boundary surface,X sub coordinate values of projection points corresponding to the cut-off points; then according to the global numbering sequence of the outer layer circulation of the boundary surface, then the sequence of the local numbering of the inner layer circulation of the boundary edge is taken to traverse the cut-off point of the boundary surface to calculate the displacement of the control point, and the point skipping is repeated;
obtaining weight
Post substitution intercept point calculationf (x,y,z) At this time->
In order to be a known quantity,f (x,y,z) Is the quantity to be calculated; calculating the displacement of grid points +.>
In (a)X=X j ，X j Representing all the cut-off points obtained by the cut-off points of the boundary surface in the step 5, the cut-off points of the non-boundary surface in the step 6 and the cut-off points of the adjacent units in the step 7, and the total number of the cut-off pointsI SUM =N+M+K；
Step 10, the deformation calculated in the step 9 is calculatedf (x,y,z) Substituting the cut-off point coordinates (x J,1 ,y J,1 ,z J,1 ) To calculate the space curved grid coordinates (x c ,y c ,z c ) The calculation formula is as follows:
（x c ,y c ,z c ）=（x J,1 ,y J,1 ,z J,1 ）+f (x,y,z),J∈{N},{M},{K}
traversing the temporary variable mesh-discrete found to obtain the deformation corresponding to each cut-off point, and adding the deformation to the initial position of the cut-off point to obtain the deformed cut-off point (the added cut-off point position is the space curved grid point), namely the space curved grid point coordinate. And then the obtained space curved grid point coordinates are stored into the cut-off point class to replace the original cut-off point coordinates.
Step 11. Higher order grid points (x) c ,y c ,z c ) As newly added grid nodes, namely, the space curved grid point coordinates obtained by calculation in the step 10 are read from the cut-off point class, as newly added grid nodes node5 to node10, the newly added grid nodes are inserted into grid data, and the grid topological relation is updated; the new space curved grid data is then formatted with a corresponding grid data file(e.g., in CGNS format) and import the output mesh data file into the visualization.
This embodiment uses 3 target objects as examples to verify the validity of the present invention; the final meshing spatial curved mesh results of the embodiments are shown in the right-hand diagrams of fig. 1 and 2, and in fig. 5.
Fig. 1 is a global surface view of a space curved grid with a sphere as a target object, a middle view is a cross-sectional view of the space curved grid, and a right view is a partial enlarged view.
Fig. 2 is a graph of the spatial curved grid results for a blunt cone (10 °). The left graph is the initial linear tetrahedral mesh obtained in the step three, and the right graph is a schematic diagram of the completed space curved mesh.
Fig. 4 and 5 are graphs of the results of an X-51A aircraft. Fig. 4 is a schematic diagram of the initial linear tetrahedral mesh obtained in step three, and fig. 5 is a schematic diagram of the completed spatial curved mesh.
From the comparison of the left graph and the right graph in fig. 2, the comparison of fig. 4 and fig. 5 can show that the space curve grid division method provided by the invention realizes the space curve grid with high precision under the condition of ensuring the quality of the grid and the unchanged basic topological structure. The existing curved grid technology basically only realizes the curved grid of the object plane, although the curved grid of the object plane improves the calculation precision at the boundary, the numerical calculation precision is often limited to the precision of the internal grid. The invention adopts a brand new strategy for selecting control points and deformation points (namely space curved grid points), expands curved grids to space inner grid to realize high-precision space curved grid division of complex model. In conclusion, the invention is practically effective.
In summary, the invention provides a new strategy for selecting the control points and the deformation points, and the grid division requirement of the complex model can be well met under the condition of not lifting the grid magnitude. Compared with the prior art: 1. the invention realizes the space curved grid division, and the prior art is basically an object plane curved grid; 2. adopting a new control point and deformation point selection strategy and a calculation method, and well meeting the grid division requirement of the complex model under the condition of not improving the grid magnitude; 3. the radial basis function is selected according to the cosine value of the included angle. The invention reduces the selection of control points and deformation points while ensuring the quality of grids, reduces a large number of position calculations in the grid dividing process, optimizes the grid dividing strategy, reduces the time consumption of grid division and reduces the calculation force requirement. The stability and the speed of grid division are greatly improved due to the reduction of deformation point calculation and the simplification of a basis function form.
Claims (5)
1. The method for meshing the three-dimensional space curve of the complex model is characterized by comprising the following steps of:
step 1, inputting a target three-dimensional physical model and grid parameters;
step 2, converting the target three-dimensional physical model into three-dimensional space geometric information, marking and storing calculation domain boundary information;
step 3, generating T initial linear tetrahedron units sequentially based on an array surface propulsion method and a three-dimensional constraint Delaunay method according to the three-dimensional space geometric information in the step 2, and simultaneously constructing a hierarchical structure of a tetrahedron, a face, edge and node points according to grid topological relations of the T initial linear tetrahedron units to store grid information;
step 4, according to the grid information of the initial linear tetrahedral units generated in the step 3 and the calculated domain boundary information in the step 2, searching the T initial linear tetrahedral units generated in the step 3 one by one;
if the T-th cell is a boundary cell, attaching a boundary mark to the cell and calculating as a T-th boundary cell, wherein T is {1,2,3, …, T }; traversing the T initial linear tetrahedral units in the step 3 to obtain total psi boundary units, psi E {1,2,3, …, psi };
step 5, searching the psi boundary units one by one according to the psi boundary units marked in the step 4, and taking out boundary surface, boundary edge and boundary point information of the psi boundary units;
the boundary surface comprises three boundary edges and a vertex P 1 (x 1 ,y 1 ,z 1 )、P 2 (x 2 ,y 2 ,z 2 ) And P 3 (x 3 ,y 3 ,z 3 ) X, y and z represent three dimensions of the point in a three-dimensional Cartesian coordinate system respectivelyCoordinate values; the first edge corresponds to the vertex P 1 And P 2 The second edge corresponds to P 1 And P 3 The third edge corresponds to P 2 And P 3 The method comprises the steps of carrying out a first treatment on the surface of the Cutting off edges on all boundary surfaces in the psi boundary units to generate boundary surface cutting-off points, wherein the tetrahedral unit surfaces comprise three edges;
the position of the point of the first boundary edge cut-off position passes through the vertex P 1 (x 1 ,y 1 ,z 1 ) And P 2 (x 2 ,y 2 ,z 2 ) And calculating the second boundary edge and the third boundary edge, wherein the calculation formula is as follows:
and the generated boundary surface intercept point P H (x n,1 ,y n,1 ,z n,1 ) Projecting to the calculated domain boundary in the step 2 to obtain a projection point P sub (x n,1 ,y n,1 ,z n,1 ) Wherein H represents the boundary surface cut-off point mark, sub represents the projection point mark, and the projection point and the boundary surface cut-off point have a one-to-one correspondence; lambda is a proportional parameter, and is adjusted based on a curvature and radial basis function of a physical model, and the adjustment range is lambda epsilon [0.5,2]；
Each time a break point is generated, the break point P will be broken H (x n,1 ,y n,1 ,z n,1 ) Logging of cut-off point data class P H (x N,1 ,y N,1 ,z N,1 ) In (1), count, P H N represents the nth cut-off point, and the maximum count value is N; all the projections are projected to the boundary surface to obtain projection points P sub (x n,1 ,y n,1 ,z n,1 ) Also stored in the proxel data class P sub (x N,1 ,y N,1 ,z N,1 ) In (1), count, P sub N represents the nth projection point, and the maximum count value is also N;
step 6, searching the surface, edge and point information of the non-boundary surface of the psi boundary units one by one according to the psi boundary unit information obtained in the step 5; obtaining non-boundary surface cutoff by the same method in step 5Point P inr (x m,1 ,y m,1 ,z m,1 ) Storing the cut-off point data class P inr (x M,1 ,y M,1 ,z M,1 ) And counting, inr represents the non-boundary surface cut-off point mark, P inr M maximum count value is marked as M;
step 7, searching adjacent units of the psi boundary units one by one according to the psi boundary unit information obtained in the step 5, taking out the information of the surfaces, edges and points of the non-adjacent surfaces in the adjacent units, and further circularly searching the adjacent units of the adjacent units if multi-layer space curved grids are needed;
the adjacent unit: two units sharing one plane in the tetrahedral units are defined as adjacent units, four adjacent units of non-boundary units are provided, and less than four but at least one adjacent units of boundary units are provided;
cutting off each edge by the same method in step 5 to obtain adjacent unit cut-off points P adj (x k,1 ,y k,1 ,z k,1 ) Logging of cut-off point data class P adj (x K,1 ,y K,1 ,z K,1 ) In (1), count, P adj K maximum count value is K, and adj represents the mark of the adjacent unit cut-off point;
step 8, traversing boundary surface cut-off point data class P H (x N,1 ,y N,1 ,z N,1 ) And proxel data class P sub (x N,1 ,y N,1 ,z N,1 ) All points in (1) calculating boundary surface cut-off point P H (x n,1 ,y n,1 ,z n,1 ) Corresponding to the projection point P sub (x n,1 ,y n,1 ,z n,1 ) The cosine value of the included angle between them is recorded as cos (theta) H-sub ；
Wherein x is n,H ,y n,H ,z n,H Three-dimensional coordinate value representing x for the nth boundary surface intercept point n,sub ,y n,sub ,z n,sub Three-dimensional coordinate values of the nth projection point;
step 9, substituting all the truncated points obtained in the steps 5, 6 and 7 into a radial basis function calculation formula, wherein the f (x, y, z) value obtained by the radial basis function calculation formula is the truncated point deformation;
the radial basis function is calculated as:
where alpha is a coefficient of which,
as weight, I N Is the total number of control points and I N ＝N；/>
Is radial basis function>
Representing the Euclidean distance from all the cut-off points to the corresponding control point, control point +.>
For boundary surface intercept point P in step 5 H (x N,1 ,y N,1 ,z N,1 ) I.e. a subset of all boundary surface intercept points;
firstly, calculating weights through control points and projection points
Weight calculation +.>
Wherein X represents the boundary surface cut-off point, so the total is I N =n; f (x, y, z) is a known quantity when the weight is calculated, and represents the displacement quantity of the control point, and f (x, y, z) takes the minimum value of the Euclidean distance from the cut-off point to the projection point on the same boundary surface;
f(x,y,z)＝min{║X b -X sub ║}
x in the case of first-order cut-off points b B epsilon {1,2,3} is three cut-off points corresponding to three boundary edges on the same boundary surface, X sub Coordinate values of projection points corresponding to the cut-off points; according to the global numbering sequence of the outer layer circulation of the boundary surface, then the sequence of the local numbering of the inner layer circulation of the boundary edge is taken to traverse the cut-off point of the boundary surface to calculate the displacement of the control point, and the point is skipped;
obtaining weight
Then substituting the cut-off points to calculate f (x, y, z), at this time +.>
F (x, y, z) is the quantity to be solved for; calculating the displacement of grid points +.>
Wherein x=x j ，X j Representing all the cut-off points obtained by the cut-off points of the boundary surface in the step 5, the cut-off points of the non-boundary surface in the step 6 and the cut-off points of the adjacent units in the step 7, and the total number of the cut-off points I SUM ＝N+M+K；
Step 10, substituting the deformation f (x, y, z) calculated in the step 9 into the coordinate (x) of the cutting point J,1 ,y J,1 ,z J,1 ) To calculate the space curved grid coordinates (x c ,y c ,z c ) The calculation formula is as follows:
(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. Higher order grid points (x) c ,y c ,z c ) As a newly added grid node, inserting the newly added grid node into grid data, and updating a grid topological relation; and then outputting the new space curved grid data in a corresponding grid data file format, and importing the output grid data file into a visualization.
2. The method for meshing the three-dimensional space curve of the complex model according to claim 1, wherein in the step 3:
firstly, adopting an array surface propelling method to compute a domain surface as an initial array surface to inwards propel and generate a hierarchical grid of a W layer, wherein W is more than or equal to 1, W is {1,2,3, …,10};
then, the initial linear tetrahedral unit is generated in the residual calculation domain by adopting the Delaunay method in the three-dimensional constraint Delaunay method by taking the current hierarchical grid surface as an initial constraint condition.
3. The method for meshing a three-dimensional space curve of a complex model according to claim 1, wherein: in the step 5, λ=1 is taken, which is the midpoint of the edge.
4. The method for meshing a three-dimensional space curve of a complex model according to claim 1, wherein:
in the step 5:
if a high-precision curved grid is to be constructed, continuing to generate a new cutoff point on the basis of the first cutoff of each boundary unit; at point P H (x n,1 ,y n,1 ,z n,1 ) And P 1 (x 1 ,y 1 ,z 1 ) And P H (x n,1 ,y n,1 ,z n,1 ) And P 2 (x 2 ,y 2 ,z 2 ) Respectively, continue according to the end points
A cut-off point is generated, at this time (x i ,y i ,z i ) For the primary cut-off point P H (x n,1 ,y n,1 ,z n,1 )，(x j ,y j ,z j ) Respectively is the vertex P 1 (x 1 ,y 1 ,z 1 ) And P 2 (x 2 ,y 2 ,z 2 ) Is a value of (2);
then the separately obtained cut-off points P H (x n,21 ,y n,21 ,z n,21 ) And P H (x n,22 ,y n,22 ,z n,22 ) Projection onto computational domain boundary to obtain P sub (x n,21 ,y n,21 ,z n,21 ) And P sub (x n,22 ,y n,22 ,z n,22 ) And storing corresponding cut-off point data class and projection point data class, thus circularly obtaining a plurality of high-order cut-off points and high-order projection points.
5. The method for meshing a three-dimensional space curve of a complex model according to claim 1, wherein: and in the step 11, the grid data file format output adopts a CGNS format.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111637505.5A CN114510775B (en) | 2021-12-30 | 2021-12-30 | Method for dividing three-dimensional space curved grids of complex model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111637505.5A CN114510775B (en) | 2021-12-30 | 2021-12-30 | Method for dividing three-dimensional space curved grids of complex model |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114510775A CN114510775A (en) | 2022-05-17 |
| CN114510775B true CN114510775B (en) | 2023-06-27 |
Family
ID=81548505
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111637505.5A Active CN114510775B (en) | 2021-12-30 | 2021-12-30 | Method for dividing three-dimensional space curved grids of complex model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114510775B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115186074B (en) * | 2022-06-24 | 2023-07-21 | 深圳市规划和自然资源数据管理中心 | Method for simulating soil pH value spatial distribution pattern based on Meta analysis |
| CN116401916B (en) * | 2023-03-20 | 2024-01-26 | 北京云境智仿信息技术有限公司 | Method, device, medium and equipment for generating high-quality three-dimensional grid |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109636913A (en) * | 2018-12-04 | 2019-04-16 | 山东理工大学 | Triangle gridding increment topology joining method based on Delaunay subdivision |
| CN112015735A (en) * | 2020-08-20 | 2020-12-01 | 西安数峰信息科技有限责任公司 | Data storage structure and data storage method of unstructured grid |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10424112B2 (en) * | 2012-05-14 | 2019-09-24 | Autodesk, Inc. | Mesh boundary smoothing |
-
2021
- 2021-12-30 CN CN202111637505.5A patent/CN114510775B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109636913A (en) * | 2018-12-04 | 2019-04-16 | 山东理工大学 | Triangle gridding increment topology joining method based on Delaunay subdivision |
| CN112015735A (en) * | 2020-08-20 | 2020-12-01 | 西安数峰信息科技有限责任公司 | Data storage structure and data storage method of unstructured grid |
Non-Patent Citations (1)
| Title |
|---|
| 用逐点插入法生成Delaunay四面体自适应网络;骆冠勇 等;计算力学学报;第24卷(第6期);917-922 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114510775A (en) | 2022-05-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114510775B (en) | Method for dividing three-dimensional space curved grids of complex model | |
| CN100468418C (en) | Method and program for generating volume data from boundary representation data | |
| CN102509339B (en) | Method for simplifying vertex clustering of three-dimensional models with texture constraint | |
| JPH077426B2 (en) | Shape design method or shape prediction method | |
| Shu | High order numerical methods for time dependent Hamilton-Jacobi equations | |
| CN111243094B (en) | Three-dimensional model accurate voxelization method based on lighting method | |
| CN114492250A (en) | Curved surface mesh generation method and system based on recursive decomposition and computer equipment | |
| CN110069854A (en) | Multiple resolution TWENO format is to the analogy method that can press flow field problems | |
| CN100454341C (en) | Process type ground fast drawing method based on fractal hierarchical tree | |
| Natarajan et al. | Simplification of three-dimensional density maps | |
| Seshadri et al. | Supporting multi-point fan design with dimension reduction | |
| Keyser et al. | Efficient and accurate B-rep generation of low degree sculptured solids using exact arithmetic: I—representations | |
| Jekeli | Spline representations of functions on a sphere for geopotential modeling | |
| Marinić-Kragić et al. | Superimposed RBF and B-spline parametric surface for reverse engineering applications | |
| Erath et al. | Integrating a scalable and effcient semi-Lagrangian multi-tracer transport scheme in HOMME | |
| Weimer et al. | Subdivision schemes for thin plate splines | |
| Zhang et al. | A nonoscillatory discontinuous Galerkin transport scheme on the cubed sphere | |
| CN113342999A (en) | Variable-resolution-ratio point cloud simplification method based on multi-layer skip sequence tree structure | |
| Langbein et al. | An efficient point location method for visualization in large unstructured grids. | |
| CN114510677B (en) | Neutron transport equation processing method based on intermittent finite element and computer program product | |
| Crampton et al. | Detecting and approximating fault lines from randomly scattered data | |
| CN113111612A (en) | Discrete point cloud repeated point fast searching method based on self-adaptive space subdivision | |
| Brandts et al. | A geometric toolbox for tetrahedral finite element partitions | |
| CN111666689A (en) | Characteristic line tracking method, reactor core neutron physical calculation method and device | |
| Mittal et al. | Mixed-Order Meshes through rp-adaptivity for Surface Fitting to Implicit Geometries |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |