WO2022126348A1 - 一种基于各向异性体调和场的附面层网格生成方法 - Google Patents

一种基于各向异性体调和场的附面层网格生成方法 Download PDF

Info

Publication number
WO2022126348A1
WO2022126348A1 PCT/CN2020/136329 CN2020136329W WO2022126348A1 WO 2022126348 A1 WO2022126348 A1 WO 2022126348A1 CN 2020136329 W CN2020136329 W CN 2020136329W WO 2022126348 A1 WO2022126348 A1 WO 2022126348A1
Authority
WO
WIPO (PCT)
Prior art keywords
mesh
boundary layer
tetrahedral
volume
harmonic field
Prior art date
Application number
PCT/CN2020/136329
Other languages
English (en)
French (fr)
Inventor
王胜法
朱一鸣
郑晓朋
雷娜
罗钟铉
陈富卫
王永杰
张帆
Original Assignee
大连理工大学
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 大连理工大学 filed Critical 大连理工大学
Priority to PCT/CN2020/136329 priority Critical patent/WO2022126348A1/zh
Priority to US17/605,659 priority patent/US20220414282A1/en
Publication of WO2022126348A1 publication Critical patent/WO2022126348A1/zh

Links

Images

Classifications

    • 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/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • 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
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Definitions

  • the invention belongs to the technical fields of computational fluid dynamics, numerical simulation, computer-aided design and manufacturing, and relates to a method for generating boundary layer grids based on anisotropic volume harmonic fields, which is suitable for generating boundary layer grids of complex curved surfaces .
  • the control of the volume harmonic field by the tensor makes the mesh of the boundary layer induced by it more flexible and controllable.
  • the boundary layer is a thin flowing layer in which the viscous force close to the object surface cannot be ignored in the flow around the high Reynolds number, and the quality of the boundary layer mesh directly determines the quality of the numerical simulation.
  • a layered anisotropic prismatic mesh perpendicular to the object must be used to capture the boundary layers near viscous walls.
  • Prismatic mesh generation around viscous walls has always been the focus of research in computational fluid dynamics, numerical simulation, computer-aided design and manufacturing technology.
  • the existing methods generally do not consider the problem of anisotropy, that is, the thickness of each layer is basically the same. When faced with some special numerical simulation requirements, it is difficult for isotropy to capture some tiny physical features. How to make the generation of prismatic meshes more flexible and controllable is the focus of the research on boundary layer meshes.
  • the present invention proposes a method for generating boundary layer meshes based on anisotropic volume harmonic fields.
  • the method belongs to the generation of boundary layer grids based on solving partial differential equations, and includes three invention contents:
  • boundary surface mesh construction based on the Minkowski sum, and the tetrahedral background mesh (discrete computational domain) generation of the boundary layer space.
  • a boundary layer mesh (prismatic mesh) generation strategy based on the advancing distance and advancing direction calculated by the anisotropic volume harmonic field.
  • a method for generating boundary layer meshes based on anisotropic volume harmonic field the steps are as follows:
  • the boundary layer space is defined as: the space between the boundary surface mesh of the calculated Minkowski sum and the original surface mesh (object surface).
  • TetGen Use the software TetGen to perform tetrahedral mesh division in the boundary layer space, and detect whether there are four points of tetrahedral elements in the tetrahedral background mesh that are simultaneously located on the boundary surface mesh of the boundary layer space. Subdivide if it exists until none of the above occurs.
  • the volume harmonic energy of the fixed original surface mesh (object surface) is a constant value a (usually 1)
  • the volume harmony energy of the fixed outer surface mesh is constant b (usually 0)
  • the edge The weights are set to the classic Cotangent Weight, and the volume harmonic field on the tetrahedral background mesh is calculated; then, a breadth-first search (BFS) is performed starting from the tetrahedral elements close to the original surface mesh (object surface), Find the tetrahedral unit that satisfies the difference between the energy values on each edge and is smaller than the threshold T slit (usually the value is [0.01, 0.1]), denoted as the set R slit ; Subdivision (the number of subdivisions generally takes the value [2,10]). Delaunay processing is used to optimize the locally subdivided tetrahedral background mesh, and then Laplace smoothing is used to
  • ⁇ i is the vertex on the tetrahedral background grid
  • [x 1 , x 2 , x 3 ] is the standard three-dimensional orthogonal frame
  • ⁇ 1 , ⁇ 2 , ⁇ 3 are three standard orthogonal frames respectively
  • the scaling factor in the direction can be regarded as placing an ellipsoid on each vertex, as shown in Figure 2, the length of the major and minor axes of the ellipse indicates the degree of control of the body harmonic field by the tensor along the major and minor axes. .
  • H is the vector composed of the values of the body harmonic field acting on each vertex
  • L is the weight matrix
  • ⁇ i is the vertex on the tetrahedral background mesh
  • e ij is the edge connecting ⁇ i and ⁇ j on the tetrahedral background mesh
  • e ij is the edge connecting ⁇ i and ⁇ j on the tetrahedral background mesh
  • N( ⁇ i ) is the set of vertices adjacent to ⁇ i
  • W(e ij ) is the edge weight for solving the Laplace equation:
  • ⁇ >0 is a control factor (generally the value is [0.001, 100])
  • the smaller the ⁇ the more the edge weight W(e ij ) is affected by the tensors T( vi ), T(v j ) Larger, on the contrary, the edge weight W(e ij ) is less affected by the tensors T(v i ) and T(v j ).
  • E is the set of edges on the tetrahedral background mesh and H( vi ) is the anisotropic volume harmonic energy at vi .
  • edge weights can be simplified as:
  • the volume harmony energy of the fixed original surface mesh (object surface) is a constant value a (usually 1), and the volume harmony energy of the fixed outer surface mesh is constant b (usually 0); according to actual needs , automatically add local anisotropy tensor control; set the maximum number of iterations T iter (usually set to 2000); set the cut-off threshold T energy (usually set to 1.0*10 -8 ) for optimizing the volume harmonic energy; through the formula ( 7), update H(v i ) iteratively; after each update of the body harmonic energy on all vertices, calculate the body harmonic energy K(H) once; iterative process until the maximum number of iterations T iter is satisfied, or K(H) The cutoff threshold T energy is reached.
  • the gradient change rate along the direction d will be limited; the boundary layer network generated by the body harmonic field controlled by the formula (8) Lattice (prismatic mesh), along the direction d, intuitively, the overall thickness of the boundary layer mesh decreases significantly.
  • the classical Cotangent Weight calculates the volume harmonic field on the tetrahedral background mesh;
  • the difference between the energy values on the strip edges is smaller than the threshold T slit (generally valued at [0.01, 0.1]) tetrahedral unit, denoted as the set R slit ;
  • the tensor is calculated, expressed as:
  • the gradient change rate of the volume harmonic field at the slit between the multi-connected branches will be limited; the volume harmonic field controlled by the formula (10) will be limited.
  • the boundary layer mesh (prismatic mesh) generated by the field, intuitively, the distortion of the boundary layer mesh is significantly reduced at the slit.
  • the advancing distance of the leading edge node is controlled by the gap between the anisotropic volume and field isosurfaces.
  • the present invention converts the expected grid thickness input by the user into sampling energy, and calculates the position of each layer node through the sampling energy.
  • the thickness of each boundary layer grid is calculated. Then, set the vertex on the object surface as the front edge node, select a front edge node with a curvature close to 0, trace the outer surrounding surface mesh along the gradient line of the volume harmonic field, and calculate the boundary layer of each layer according to the calculation.
  • each tetrahedral element is a linear space, and the forward direction of the front node is Under the guidance of , the position of the leading edge node after advancing can be easily determined by sampling the energy.
  • the forward direction of the front node is obtained by performing weighted Laplace smoothing on the gradient direction of the volume harmonic field.
  • the front node Under the guidance of the advancing distance and the advancing direction, the front node can calculate a new family of advancing positions; the boundary layer mesh is obtained by the directional connection of the advancing positions of all the front nodes according to the original surface mesh (object surface) topology.
  • the method for generating boundary layer meshes based on anisotropic volume harmonic fields proposed by the present invention includes three beneficial effects:
  • the traditional boundary layer mesh generation method based on solving partial differential equations only considers global information, and lacks local control and flexibility.
  • the present invention realizes the control of the volume harmonic field by automatically constructing the local anisotropy tensor. On the one hand, it enables it to perceive the local geometric information, and strengthens the control and flexibility of the local grid generation (mainly for the attached grids).
  • the generated boundary layer mesh can be made dense in one or several directions, so that subtle physical features can be captured.
  • the present invention uses the distance between the isosurfaces of the body harmonic field as a guide to control the advancing distance of the front nodes of each layer, that is,
  • Each layer of nodes in the boundary layer grid has the same volume harmonic energy value in the tetrahedral background grid (discrete computational domain), and they are all on the same layer volume harmonic field isosurface.
  • the coupling relationship between the volume harmonic field and the thickness of the boundary layer mesh is established, so that the thickness of the boundary layer mesh is more flexible and controllable, so that the local control of the boundary layer can be achieved by locally controlling the volume harmonic field.
  • the purpose of layer mesh thickness is to meet complex practical needs.
  • boundary layer meshes Prior meshes
  • partial differential equations generally directly uses the gradient direction as the advancing direction, but it is easy to introduce a large number of prismatic elements with negative volume or zero volume at concave edges and grooves.
  • weighted Laplacian smoothing to the gradient direction to get a smoother forward direction.
  • the focus is on the design and selection of the weights, which will directly affect the result of the boundary layer mesh generation.
  • the strategy adopted by the present invention uses the gradient direction to calculate the initial advance position, and directly calculates the quality of the generated prismatic unit in combination with the current position information, and uses the quality of the current prismatic unit as the basis for weight setting, which effectively avoids the Generation of negative volume prismatic elements.
  • Fig. 1 is the algorithm flow chart of the present invention
  • Fig. 2 is the schematic diagram of the tensor acting on the volume harmonic field on the vertex
  • Figure 3 is a schematic diagram of a Laplacian smooth weight design for optimizing the forward direction
  • Figure 4 is a schematic diagram of the plane model generating the boundary layer mesh based on the anisotropic volume harmonic field, (a) the original surface (object surface) mesh of the plane model; (b) the surface mesh surrounding the plane model; (c) the plane The original (object surface) mesh of the model; (d) the cross section of the standard volume harmonic field isosurface of the aircraft model; (e) the cross section of the anisotropic volume harmonic field isosurface of the aircraft model; (f) the aircraft model based on each The boundary layer network of anisotropic body harmonic fields.
  • the algorithm flow of the present invention is shown in Figure 1, which generally includes five steps: constructing the boundary surface mesh of the Minkowski sum; generating a tetrahedral background mesh (discrete computational domain) for the boundary layer space; constructing an isotropic Anisotropic tensor; computes anisotropic volume harmonic fields; generates boundary layer meshes (prismatic meshes) based on anisotropic volume harmonic fields.
  • the input of the algorithm of the present invention includes 1 original surface mesh (object surface) and 3 parameters.
  • the original surface mesh can be a triangular mesh or a quadrilateral mesh; the three parameters are the mesh thickness L 1 of the first layer of the boundary layer, the thickness growth factor ⁇ of the boundary layer, and the boundary layer The number of layers n of the grid.
  • the specific implementation of generating a boundary layer mesh based on an aircraft model based on anisotropic volume harmonic field is an example of the present invention, as shown in FIG. 4 , and the specific steps are as follows:
  • Input aircraft model triangular mesh
  • TetGen Use the commercial software TetGen to perform tetrahedral mesh division for the boundary layer space between the original surface mesh (object surface) and the surrounding surface mesh to obtain a tetrahedral background mesh (discrete computational domain); detection Whether there are four points of the tetrahedral element in the tetrahedral background mesh are located on the boundary surface mesh of the boundary layer space at the same time, if so, subdivide until the above situation does not occur;
  • the volume harmonic energy on the original surface mesh (object surface) is fixed to 1, and the volume harmonic energy on the fixed outer surface mesh is 0; based on the anisotropy tensor, each edge is calculated by formula (6).
  • the control factor ⁇ is set to 0.05; the maximum number of iterations is set to 2000, and the energy cutoff threshold is set to 1.0*10 -8 ;
  • the body harmonic energy value on the vertex is iteratively updated by formula (7), each iteration is 50 times , calculate the new body harmonic energy K(H) according to formula (5); if the difference between the current body harmonic energy value and the previous body harmonic energy value is less than the cut-off threshold, the iteration ends, otherwise, continue to iterate until the maximum number of iterations is reached ;

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Complex Calculations (AREA)
  • Image Generation (AREA)

Abstract

本发明公开一种基于各向异性体调和场的附面层网格生成方法，属于计算流体力学、数值模拟仿真、计算机辅助设计、与制造技术领域。首先采用闵科夫斯基和的边界曲面网格构造求解体调和场时所需的四面体背景网格，然后根据需求自动化地添加各向异性张量，并在张量的控制下计算各向异性体调和场，最后结合特殊的加权拉普拉斯平滑计算生成附面层网格所需的前进方向。本发明基于闵科夫斯基和边界曲面网格构造四面体背景网格的策略降低了计算耗时和内存浪费，通过自动化添加各向异性张量能够可控地局部调节附面层网格厚度，结合特殊的加权拉普拉斯光滑对前进方向进行优化显著地提高了附面层网格的生成质量。

Description

一种基于各向异性体调和场的附面层网格生成方法 技术领域
本发明属于计算流体力学、数值模拟仿真、计算机辅助设计，以及制造技术领域，涉及一种基于各向异性体调和场生成附面层网格的方法，适用于复杂曲面的附面层网格生成。通过张量对体调和场的控制，使其诱导生成的附面层网格更具有灵活性和可控性。
背景技术
附面层是高雷诺数绕流中紧贴物面的粘性力不可忽略的流动薄层，而附面层网格的质量直接决定着数值模拟效果的好坏。在高雷诺数流动的气动模拟中，必须采用垂直于物体的层状各向异性棱柱网格来捕捉粘性壁面附近的附面层。粘性壁面周围的棱柱网格生成一直是计算流体力学、数值模拟仿真、计算机辅助设计与制造技术领域的研究重点。棱柱网格生成的主流方法有两种：前沿结点推进方法，基于求解偏微分方程方法。但是已有的方法一般不考虑各向异性的问题，即每一层增长的厚度基本都一样，面对一些特殊的数值模拟需求时，各向同性很难捕捉一些微小的物理特征。如何使得生成棱柱网格更具有灵活性和可控性，是附面层网格研究的重心。
发明内容
基于上述的问题，本发明提出一种基于各向异性体调和场生成附面层网格的方法。本方法属于基于求解偏微分方程生成附面层网格，包含3个发明内容：
1.基于闵可夫斯基和的边界曲面网格构造，以及附面层空间的四面体背景网格(离散计算域)生成。
2.基于局部张量控制的各向异性体调和场的计算。
3.基于各向异性体调和场计算的前进距离和前进方向的附面层网格(棱柱 网格)生成策略。
本发明的技术方案：
一种基于各向异性体调和场的附面层网格生成方法，步骤如下：
(1)基于闵可夫斯基和的边界曲面网格构造，以及附面层空间的四面体背景网格(离散计算域)生成，具体步骤如下：
a)输入原始曲面网格(一般为三角形网格或者四边形网格)和半径为r(一般取值为能够包裹原始曲面网格的最小长方体的对角线长度乘以一个系数c，c一般取值在[0.05,0.3])的球网格，计算闵可夫斯基和的边界曲面网格。
b)对初步得到的闵可夫斯基和的边界曲面网格执行网格优化处理，包含非流形消除和自相交消除，最终得到二维流形的边界曲面网格。
c)定义附面层空间为：计算得到的闵可夫斯基和的边界曲面网格与原始曲面网格(物面)之间的空间。
d)对附面层空间使用软件TetGen进行四面体网格剖分，检测四面体背景网格中是否存在四面体单元的四个点同时位于附面层空间的边界曲面网格上。如果存在则进行细分，直到无上述情况出现。
e)对四面体背景网格的狭缝处进行局部细分。首先，固定原始曲面网格(物面)的体调和能量为常值a(一般取值为1)，固定外包围曲面网格的体调和能量为常值b(一般取值为0)，边权重设置为经典的余切权重(Cotangent Weight)，计算四面体背景网格上的体调和场；然后，从靠近原始曲面网格(物面)的四面体单元开始做广度优先搜索(BFS)，寻找满足每条边上的能量值之差同时小于阈值T slit(一般取值在[0.01,0.1])的四面体单元，记作集合R slit；最后，对集合R slit内的四面体单元进行细分(细分次数一般取值在[2,10])。对局部细分后的四面体背景网格使用delaunay处理进行优化，再使用拉普拉斯平滑辅助优化， 最终得到高质量的四面体背景网格。
(2)基于局部张量控制的各向异性体调和场的计算，具体内容如下：
2.1)四面体背景网格顶点上张量的定义
T(υ i)＝γ 1x 1x 1 Τ+γ 2x 2x 2 Τ+γ 3x 3x 3 Τ；(1)
其中，υ i是四面体背景网格上的顶点，[x 1,x 2,x 3]是标准的三维正交标架，γ 1,γ 2,γ 3分别作为标准正交标架三个方向上的缩放因子。直观上，张量对体调和场的控制可以看作在每一个顶点上放置一个椭球，如图2所示，椭圆的长短轴边长表示着张量沿着长短轴方向对体调和场的控制程度。
2.2)各向异性体调和场在四面体背景网格顶点上的定义
LH＝0；(2)
其中，H是体调和场作用在每一个顶点上的值组合而成的向量；L是权重矩阵，表达式为：
Figure PCTCN2020136329-appb-000001
其中，υ i是四面体背景网格上的顶点，e ij是四面体背景网格上连接υ i和υ j的边，e ij是四面体背景网格上连接υ i和υ j的边，N(υ i)是与υ i相邻顶点的集合，W(e ij)是求解拉普拉斯方程的边权重：
Figure PCTCN2020136329-appb-000002
其中，δ>0，是一个控制因子(一般取值在[0.001，100])，δ越小则边权重W(e ij)受到张量T(v i)，T(v j)的影响越大，反之，则边权重W(e ij)受到张量T(v i)，T(v j)的影响越小。
2.3)四面体背景网格中的各向异性体调和场的计算
定义各向异性体调和能量K(H)：
Figure PCTCN2020136329-appb-000003
其中，E是四面体背景网格上的边集合，H(v i)是在v i处的各向异性体调和能量。
在基于迭代法计算的各向异性体调和场的框架下，边权重可以简化为：
W(e ij)＝exp(T(e ij)/δ)；(6)
在基于迭代法计算的各向异性体调和场的框架下，对于顶点上体调和能量的优化，表示为：
Figure PCTCN2020136329-appb-000004
基于迭代法的各向异性体调和场的计算的算法流程如下：
固定原始曲面网格(物面)的体调和能量为常值a(一般取值为1)，固定外包围曲面网格的体调和能量为常值b(一般取值为0)；根据实际需求，自动化地添加局部各向异性张量控制；设置最大的迭代次数T iter(一般设置为2000)；设置优化体调和能量的截断阈值T energy(一般设置为1.0*10 -8)；通过公式(7)，迭代地更新H(v i)；每次更新完所有顶点上的体调和能量，计算一次体调和能量K(H)；迭代过程直到满足最大的迭代次数T iter，或者K(H)达到截断阈值T energy。
2.4)自动化的各向异性张量控制的构造
2.4.1)一般情况下，限制体调和场沿着特定方向d的梯度变化速率，表示为：
Figure PCTCN2020136329-appb-000005
将公式(8)构造的张量作用于体调和场的计算中，则沿着方向d上的梯度变化速率将受到限制；基于公式(8)张量控制的体调和场生成的附面层网格(棱柱网格)，沿着方向d，直观上，附面层网格的总体厚度明显减小。
2.4.2)限制凹边和凹槽处的体调和场梯度变化速率，为了使得计算得到的体调和场等值面更贴合于物面，避免诱导生成的附面层网格在凹边和凹槽处产生较大畸变，自动构造局部张量方式如下：
首先，设置原始曲面网格(物面)的能量为常量a(一般取值为1)，外包围曲面网格的能量常量b(一般取值为0)，要求a>b，边权重设置为经典的余切权重(Cotangent Weight)，计算四面体背景网格上的体调和场；然后，从靠近原始曲面网格(物面)的四面体单元开始做广度优先搜索(BFS)，寻找满足每条边上的能量值之差同时小于阈值T slit(一般取值在[0.01,0.1])的四面体单元，记作集合R slit；最后，计算张量，表示为：
Figure PCTCN2020136329-appb-000006
将公式(9)构造的张量作用于体调和场的计算中，则在凹边和凹槽处的体调和场梯度变化速率将受到限制；基于公式(9)张量控制的体调和场生成的附面层网格(棱柱网格)，在凹边和凹槽处，直观上，附面层网格的畸变明显减小。
2.4.3)限制多连通分支之间狭缝处的体调和场梯度变化速率，延缓体调和场在狭缝处鞍点的产生，提高诱导生成附面层网格在多连通分支之间狭缝处的质量。自动构造局部张量方式如下，设两个相互靠近的体模型P，Q：首先，分别设H 1(P)＝a，H 1(Q)＝b和H 2(P)＝b，H 2(Q)＝a作为狄利克雷边界条件(a一般取值为1；b一般取值为0)，边权重设置为经典的余切权重(Cotangent Weight)，计算两个标准体调和场H 1，H 2；然后，计算张量，表示为：
Figure PCTCN2020136329-appb-000007
将公式(10)构造的张量作用于体调和场的计算中，则在多连通分支之间狭缝处的体调和场梯度变化速率将受到限制；基于公式(10)张量控制的体调和场 生成的附面层网格(棱柱网格)，直观上，在狭缝处附面层网格的畸变明显减小。
(3)基于各向异性体调和场计算的前进距离和前进方向的附面层网格(棱柱网格)生成策略，具体内容如下：
3.1)前沿结点前进距离的计算
前沿结点的前进距离，通过各向异性体调和场等值面之间的间隙来控制。本发明将用户输入的期望网格厚度转化为采样能量，通过采样能量计算每一层结点的位置。具体实现方式：
首先，根据用户输入的第一层附面层厚度L 1，附面层厚度增长速度因子α，附面层层数n，计算出每一层附面层网格的厚度。然后，将物面上的顶点设置为前沿结点，选择一个曲率接近0的前沿结点沿着体调和场的梯度线追溯到外包围曲面网格，并根据计算出的每一层附面层网格厚度，在体调和场中提取n个采样能量；最后，对于离散化在四面体背景网格中的体调和场，每一个四面体单元内都是一个线性空间，在前沿结点前进方向的引导下，前沿结点前进后的位置能够很容易通过采样能量确定。
3.2)前沿结点前进方向的计算
前沿结点的前进方向，通过对体调和场的梯度方向进行加权拉普拉斯光滑化后得到。具体实现方式：
首先，计算前沿结点当前位置的梯度方向，一般为所在四面体单元内等值面的法向量方向；然后，记当前位置为p i，在梯度方向和下一个采样能量的引导下，计算出下一个位置
Figure PCTCN2020136329-appb-000008
如图3所示；最后，将拉普拉斯光滑的权重设置为
Figure PCTCN2020136329-appb-000009
其中，p(一般取值为4)和q(一般取值为2)是两个控制参数，另外
Figure PCTCN2020136329-appb-000010
在这种权重表示下，能够直接检测出在前进方向和前进距离引导下生成的附面层网格是否存在负体积单元，有效保障了生成附面层网格的质量；一般拉普拉斯光滑次数设置为100次。
3.3)附面层网格(棱柱网格)的生成
前沿结点在前进距离和前进方向的引导下，计算得到一族新的前进位置；附面层网格由所有前沿结点的前进位置根据原始曲面网格(物面)拓扑的有向连接得到。
本发明的有益效果：
基于上述的发明内容，本发明提出的基于各向异性体调和场生成附面层网格的方法包含3个有益效果：
(1)传统构造四面体背景网格(离散计算域)的方法采用长方体或者球体作为外包围曲面网格，这样的策略引入大量冗余的四面体单元，带来额外的计算量消耗。基于闵可夫斯基和边界曲面网格的四面体背景网格(离散计算域)生成能够有效地消除冗余四面体单元，从而提高内存的利用率和算法的执行效率。
(2)传统基于求解偏微分方程的附面层网格生成方法只考虑全局信息，欠缺局部上的控制力和灵活性。本发明通过自动化地构造局部各向异性张量来实现对体调和场的控制，一方面，使其能够感知局部的几何信息，加强对局部生成网格的控制力和灵活性(主要是针对附面层网格厚度的控制力)；另一方面，可以根据实际需求，使得生成的附面层网格沿着某个或者某几个方向变得稠密，以至于能够捕捉细微的物理特征。
(3)本发明中对于前沿结点的前进距离和前进方向的计算策略具有以下好 处：
a)为了让体调和场的构造与附面层网格的生成更好地结合，本发明将体调和场等值面之间的距离作为引导，控制每一层前沿结点的前进距离，即附面层网格中的每一层结点，在四面体背景网格(离散计算域)中的体调和能量值是相等的，它们都处于同一层体调和场等值面上。基于这种策略，建立体调和场与附面层网格厚度之间的耦合关系，使附面层网格厚度更灵活、更可控，以致于能够通过局部控制体调和场达到局部控制附面层网格厚度的目的，以满足复杂的实际需求。
传统基于偏微分方程生成附面层网格(棱柱网格)的方法，一般直接使用梯度方向作为前进方向，然而容易在凹边和凹槽处引入大量负体积或者零体积的棱柱体单元。除此之外，存在许多工作将加权拉普拉斯平滑应用于梯度方向上，得到更光滑的前进方向。对于加权拉普拉斯平滑策略，重点是权重的设计和选择，这将直接影响附面层网格生成的结果。本发明所采取的策略，利用梯度方向计算初始前进位置，并结合当前位置信息直接计算所生成棱柱体单元的质量，将当前棱柱体单元的质量作为权重设置的根据，一定程度上，有效地避免了负体积棱柱体单元的生成。
附图说明
图1为本发明的算法流程图；
图2为张量在顶点上作用于体调和场的示意图；
图3为优化前进方向的拉普拉斯光滑的权重设计的示意图；
图4为飞机模型基于各向异性体调和场生成附面层网格的示意图，(a)飞机模型原始曲面(物面)网格；(b)飞机模型外包围曲面网格；(c)飞机模型原始(物面)网格；(d)飞机模型标准体调和场等值面的横截面；(e)飞 机模型各向异性体调和场等值面的横截面；(f)飞机模型基于各向异性体调和场的附面层网络。
具体实施方式
以下结合附图和技术方案，进一步详述本发明的具体实施方式。
本发明的算法流程如图1所示，总体包含5个步骤：构造闵科夫斯基和的边界曲面网格；对附面层空间生成四面体背景网格(离散计算域)；构造各向异性张量；计算各向异性体调和场；基于各向异性体调和场生成附面层网格(棱柱网格)。本发明的算法的输入包含1份原始曲面网格(物面)和3个参数。其中，原始曲面网格(物面)可以为三角形网格或者四边形网格；3个参数分别为第一层附面层网格厚度L 1，附面层网格厚度增长因子α以及附面层网格的层数n。
本实施例以飞机模型基于各向异性体调和场生成附面层网格的具体实施为本发明的示例，如图4所示，具体步骤如下：
1.输入飞机模型(三角形网格)；输入参数L 1＝1.0*10 -2，α＝1.15，n＝60；根据输入参数L 1，α，n，计算每一层附面层网格的期望厚度{L 1,L 2,...,L n}，其中，L i+1＝L i*α；计算附面层网格的期望总体厚度
Figure PCTCN2020136329-appb-000011
2.设置小球网格半径为L total*1.5，计算原始曲面网格(物面)的闵科夫斯基和的边界曲面网格作为外包围曲面网格；消除外包围曲面网格的非流型区域和自相交区域；
3.为原始曲面网格(物面)和外包围曲面网格之间的附面层空间使用商业化软件TetGen进行四面体网格剖分，得到四面体背景网格(离散计算域)；检测四面体背景网格中是否存在四面体单元的四个点同时位于附面层空间的边界曲面网格上，如果存在则进行细分，直到无上述情况出现；
4.固定原始曲面网格(物面)上的体调和能量为1，固定外包围曲面网格上的体调和能量为0，边权重设置为经典的余切权重(Cotangent Weight)，计算四面体背景网格上的体调和场；以依附着原始曲面网格(物面)的四面体单元为起始位置执行广度优先搜索(BFS)，寻找满足每条边上的能量值之差同时小于阈值0.05的四面体单元，记作集合R slit；根据公式(9)计算局部张量；
5.固定原始曲面网格(物面)上的体调和能量为1，固定外包围曲面网格上的体调和能量为0；基于各向异性张量，通过公式(6)计算出每一条边上的权重，其中，控制因子δ设置为0.05；设置最大迭代次数为2000，设置能量截断阈值为1.0*10 -8；通过公式(7)迭代更新顶点上的体调和能量值，每迭代50次，根据公式根据公式(5)计算新的体调和能量K(H)；若当前体调和能量值与上一次体调和能量值之差小于截断阈值则迭代结束，否则继续迭代，直到达到最大迭代次数；
6.将原始曲面网格(物面)上的顶点设置为前沿结点，选择一个曲率最接近于0的前沿结点沿着体调和场的梯度线追溯到外包围曲面网格，根据计算出的每一层附面层网格厚度{L 1,L 2,...,L n}，在追溯的轨迹上提取n个采样能量{K 1,K 2,...,K n}；
7.将前沿结点当前所在位置的梯度方向作为初始前进方向，以公式(11)作为加权拉普拉斯光滑的权重对前进方向进行优化，平滑次数设置为100次；
8.基于优化后的前进方向和采样能量K i，计算出前沿结点推进后的位置；每推进前沿结点一次，重新计算前进方向；根据原始曲面网格(物面)的拓扑对所有前沿结点所分别对应的n个推进位置进行简单的有向连接，得到附面层网格(棱柱网格)。

Claims (1)

  1. 一种基于各向异性体调和场的附面层网格生成方法，其特征在于，步骤如下：
    (1)基于闵可夫斯基和的边界曲面网格构造，以及附面层空间的四面体背景网格生成
    a)输入原始曲面网格和半径为r的球网格，计算闵可夫斯基和的边界曲面网格；
    b)对初步得到的闵可夫斯基和的边界曲面网格执行网格优化处理，包含非流形消除和自相交消除，最终得到二维流形的边界曲面网格；
    c)定义附面层空间为：计算得到的闵可夫斯基和的边界曲面网格与原始曲面网格之间的空间；
    d)对附面层空间进行四面体网格剖分，检测四面体背景网格中是否存在四面体单元的四个点同时位于附面层空间的边界曲面网格上；如果存在则进行细分，直到无上述情况出现；
    e)对四面体背景网格的狭缝处进行局部细分，并对局部加密后的四面体背景网格进行优化，再使用拉普拉斯平滑辅助优化，最终得到高质量的四面体背景网格；
    (2)基于局部张量控制的各向异性体调和场的计算
    1)张量在四面体背景网格顶点上的定义：
    T(υ i)＝γ 1x 1x 1 Τ+γ 2x 2x 2 Τ+γ 3x 3x 3 Τ；(1)
    其中，υ i是四面体背景网格上的顶点，[x 1,x 2,x 3]是标准的三维正交标架，γ 1,γ 2,γ 3分别作为标准正交标架三个方向上的缩放因子；
    2)各向异性体调和场在四面体背景网格顶点上的定义：
    LH＝0；(2)
    其中，H是各向异性体调和场作用在每一个顶点上的值组合而成的向量；L是 权重矩阵，表达式为：
    Figure PCTCN2020136329-appb-100001
    其中，υ i是四面体背景网格上的顶点，e ij是四面体背景网格上连接υ i和υ j的边，N(υ i)是与υ i相邻顶点的集合，W(e ij)是求解拉普拉斯方程的边权重：
    Figure PCTCN2020136329-appb-100002
    其中，δ>0，是一个控制因子，δ越小则边权重W(e ij)受到张量T(v i)，T(v j)的影响越大，反之，则边权重W(e ij)受到张量T(v i)，T(v j)的影响越小；
    3)四面体背景网格中的各向异性体调和场的计算：
    定义各向异性体调和能量K(H)：
    Figure PCTCN2020136329-appb-100003
    其中，E是四面体背景网格上的边集合，H(v i)是在v i处的各向异性体调和能量；
    在基于迭代法计算的各向异性体调和场的框架下，边权重简化为：
    W(e ij)＝exp(T(e ij)/δ)；(6)
    在基于迭代法计算的各向异性体调和场的框架下，对于顶点上体调和能量的优化，表示为：
    Figure PCTCN2020136329-appb-100004
    基于迭代法的各向异性体调和场的计算的算法流程如下：
    固定原始曲面网格的体调和能量为常值a，固定外包围曲面网格的体调和能量为常值b；根据实际需求，自动化地添加局部各向异性张量控制；设置最大的 迭代次数T iter；设置优化体调和能量的截断阈值T energy；通过公式(7)，迭代地更新H(v i)；每次更新完所有顶点上的体调和能量，计算一次体调和能量K(H)；迭代过程直到满足最大的迭代次数T iter，或者K(H)达到截断阈值T energy；
    4)自动化的各向异性张量控制的构造：
    4.1)限制体调和场沿着特定方向d的梯度变化速率，表示为：
    Figure PCTCN2020136329-appb-100005
    将公式(8)构造的张量作用于体调和场的计算中，则沿着方向d上的梯度变化速率将受到限制；基于公式(8)张量控制的体调和场生成的附面层网格，沿着方向d，直观上，附面层网格的总体厚度明显减小；
    4.2)限制凹边和凹槽处的体调和场梯度变化速率，为了使得计算得到的体调和场等值面更贴合于物面，避免诱导生成的附面层网格在凹边和凹槽处产生较大畸变，自动构造局部张量方式如下：
    首先，设置原始曲面网格的体调和能量为常量a，外包围曲面网格的能量常量b，要求a>b，计算四面体背景网格上的体调和场；然后，从靠近原始曲面网格的四面体单元开始做广度优先搜索，寻找满足每条边上的能量值之差同时小于阈值T slit的四面体单元，记作集合R slit；最后，计算张量，表示为：
    Figure PCTCN2020136329-appb-100006
    将公式(9)构造的张量作用于体调和场的计算中，则在凹边和凹槽处的体调和场梯度变化速率将受到限制；基于公式(9)张量控制的体调和场生成的附面层网格，在凹边和凹槽处，直观上，附面层网格的畸变明显减小；
    4.3)限制多连通分支之间狭缝处的体调和场梯度变化速率，延缓体调和场在狭缝处鞍点的产生，提高诱导生成附面层网格在多连通分支之间狭缝处的质 量；自动构造局部张量方式如下，设两个相互靠近的体模型P，Q：首先，分别设H 1(P)＝a，H 1(Q)＝b和H 2(P)＝b，H 2(Q)＝a作为狄利克雷边界条件，计算两个标准体调和场H 1，H 2；然后，计算张量，表示为：
    Figure PCTCN2020136329-appb-100007
    将公式(10)构造的张量作用于体调和场的计算中，则在多连通分支之间狭缝处的体调和场梯度变化速率将受到限制；基于公式(10)张量控制的体调和场生成的附面层网格，直观上，在狭缝处附面层网格的畸变明显减小；
    (3)基于各向异性体调和场计算的前进距离和前进方向的附面层网格生成策略
    前沿结点前进距离的计算：
    前沿结点的前进距离，通过各向异性体调和场等值面之间的间隙来控制；将用户输入的期望网格厚度转化为采样能量，通过采样能量计算每一层结点的位置；具体实现方式：
    首先，根据用户输入的第一层附面层厚度L 1，附面层厚度增长速度因子α，附面层层数n，计算出每一层附面层网格的厚度；然后，将物面上的顶点设置为前沿结点，选择一个曲率接近0的前沿结点沿着体调和场的梯度线追溯到外包围曲面网格，并根据计算出的每一层附面层网格厚度，在体调和场中提取n个采样能量；最后，对于离散化在四面体背景网格中的体调和场，每一个四面体单元内都是一个线性空间，在前沿结点前进方向的引导下，前沿结点前进后的位置能够很容易通过采样能量确定；
    a)前沿结点前进方向的计算：
    前沿结点的前进方向，通过对体调和场的梯度方向进行加权拉普拉斯光滑化后得到；具体实现方式：
    首先，计算前沿结点当前位置的梯度方向，为所在四面体单元内等值面的法向量方向；然后，记当前位置为p i，在梯度方向和下一个采样能量的引导下，计算出下一个位置
    Figure PCTCN2020136329-appb-100008
    最后，将拉普拉斯光滑的权重设置为
    Figure PCTCN2020136329-appb-100009
    其中，p和q是两个控制参数，另外
    Figure PCTCN2020136329-appb-100010
    在这种权重表示下，能直接检测出在前进方向和前进距离引导下生成的附面层网格是否存在负体积单元，有效保障了生成附面层网格的质量；
    b)附面层网格的生成：
    前沿结点在前进距离和前进方向的引导下，计算得到一族新的前进位置；附面层网格由所有前沿结点的前进位置根据原始曲面网格拓扑的有向连接得到。
PCT/CN2020/136329 2020-12-15 2020-12-15 一种基于各向异性体调和场的附面层网格生成方法 WO2022126348A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/CN2020/136329 WO2022126348A1 (zh) 2020-12-15 2020-12-15 一种基于各向异性体调和场的附面层网格生成方法
US17/605,659 US20220414282A1 (en) 2020-12-15 2020-12-15 Boundary layer mesh generation method based on anisotropic volume harmonic field

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2020/136329 WO2022126348A1 (zh) 2020-12-15 2020-12-15 一种基于各向异性体调和场的附面层网格生成方法

Publications (1)

Publication Number Publication Date
WO2022126348A1 true WO2022126348A1 (zh) 2022-06-23

Family

ID=82059811

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/136329 WO2022126348A1 (zh) 2020-12-15 2020-12-15 一种基于各向异性体调和场的附面层网格生成方法

Country Status (2)

Country Link
US (1) US20220414282A1 (zh)
WO (1) WO2022126348A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114820991A (zh) * 2022-06-30 2022-07-29 中国空气动力研究与发展中心计算空气动力研究所 一种非结构附面层网格交叉处理方法及装置

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117131832B (zh) * 2023-10-23 2024-02-02 巨霖科技(上海)有限公司 一种仿真元器件的构建方法、装置及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105550465A (zh) * 2016-01-06 2016-05-04 五邑大学 一种基于橫风非定常理论的车厢横断面优化方法
CN106844994A (zh) * 2017-02-09 2017-06-13 苏州大学 本构模型与有限元结合的脉络膜新生血管生长预测方法
CN107798730A (zh) * 2017-10-27 2018-03-13 中国空气动力研究与发展中心计算空气动力研究所 一种结构网格附面层网格自动生成方法
US20190347854A1 (en) * 2018-05-11 2019-11-14 Canon Kabushiki Kaisha System and method for processing a graphic object

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105550465A (zh) * 2016-01-06 2016-05-04 五邑大学 一种基于橫风非定常理论的车厢横断面优化方法
CN106844994A (zh) * 2017-02-09 2017-06-13 苏州大学 本构模型与有限元结合的脉络膜新生血管生长预测方法
CN107798730A (zh) * 2017-10-27 2018-03-13 中国空气动力研究与发展中心计算空气动力研究所 一种结构网格附面层网格自动生成方法
US20190347854A1 (en) * 2018-05-11 2019-11-14 Canon Kabushiki Kaisha System and method for processing a graphic object

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114820991A (zh) * 2022-06-30 2022-07-29 中国空气动力研究与发展中心计算空气动力研究所 一种非结构附面层网格交叉处理方法及装置
CN114820991B (zh) * 2022-06-30 2022-09-16 中国空气动力研究与发展中心计算空气动力研究所 一种非结构附面层网格交叉处理方法及装置

Also Published As

Publication number Publication date
US20220414282A1 (en) 2022-12-29

Similar Documents

Publication Publication Date Title
CN112613206B (zh) 一种基于各向异性体调和场的附面层网格生成方法
Sun et al. Structural shape optimization by IGABEM and particle swarm optimization algorithm
WO2022126348A1 (zh) 一种基于各向异性体调和场的附面层网格生成方法
Ito Challenges in unstructured mesh generation for practical and efficient computational fluid dynamics simulations
Ling et al. Spectral quadrangulation with feature curve alignment and element size control
Ceze et al. Drag prediction using adaptive discontinuous finite elements
Loppi et al. Locally adaptive pseudo-time stepping for high-order flux reconstruction
Lu et al. Parallel mesh adaptation for high-order finite element methods with curved element geometry
Ebrahimi et al. Shape modeling based on specifying the initial B-spline curve and scaled BFGS optimization method
Chen et al. MGNet: a novel differential mesh generation method based on unsupervised neural networks
Du et al. Super resolution generative adversarial networks for multi-fidelity pressure distribution prediction
CN110689620A (zh) 一种多层次优化的网格曲面离散样条曲线设计方法
Peng et al. An automatic isotropic/anisotropic hybrid grid generation technique for viscous flow simulations based on an artificial neural network
Wang et al. A projective transformation-based topology optimization using moving morphable components
Toal et al. Geometric filtration using POD for aerodynamic design optimization
Yang et al. Isogeometric double-objective shape optimization of free-form surface structures with Kirchhoff–Love shell theory
Shen et al. Hexahedral mesh adaptation based on posterior-error estimation
Ali et al. Optimal multi-block mesh generation for CFD
Loseille et al. Developments on the P^ 2 cavity operator and Bézier Jacobian correction using the simplex algorithm.
Chang Development of Physics-Based Transition Models for Unstructured-Mesh CFD Codes Using Deep Learning Models
Chauhan et al. Wing shape optimization using FFD and twist parameterization
Xiao et al. Smooth gradation of anisotropic meshes using log–euclidean metrics
Du et al. A Fully Automated Adaptive Sampling Strategy for Reduced-Order Modeling of Flow Fields
Li et al. A fully differentiable GNN-based PDE Solver: With Applications to Poisson and Navier-Stokes Equations
CN118070621B (zh) 固壁边界的处理方法、装置、终端设备及存储介质

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20965366

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20965366

Country of ref document: EP

Kind code of ref document: A1