CN114510775B - 一种复杂模型三维空间曲网格划分方法 - Google Patents

一种复杂模型三维空间曲网格划分方法 Download PDF

Info

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
Application number
CN202111637505.5A
Other languages
English (en)
Other versions
CN114510775A (zh
Inventor
黄澍程
徐立
黄桃
李斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202111637505.5A priority Critical patent/CN114510775B/zh
Publication of CN114510775A publication Critical patent/CN114510775A/zh
Application granted granted Critical
Publication of CN114510775B publication Critical patent/CN114510775B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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

本发明属于飞行器设计、飞行器数值模拟和数值计算技术、网格划分技术领域，具体涉及一种复杂模型三维空间曲网格划分方法。本发明提出一种新的控制点和空间曲网格点选取策略，在不提升网格量级的情况下很好的满足复杂模型的网格划分需求。在保证网格质量的同时减少了控制点和形变点的选取，减少了网格划分过程中大量的位置计算，优化了网格划分策略。并提出了基于夹角余弦值的基函数选择方式，在保证网格质量的前提下，使用形式简单的低次多项式基函数的能大量减少形变量所需的计算，减少了网格划分耗时以及降低了算力需求。有效解决了现有复杂模型网格划分过程中存在网格划分普适性较差，网格质量不佳的问题。

Description

一种复杂模型三维空间曲网格划分方法
技术领域
本发明属于飞行器设计、飞行器数值模拟和数值计算技术、网格划分技术领域，具体涉及一种复杂模型三维空间曲网格划分方法。
背景技术
随着计算机技术的不断发展，数值模拟技术和数值计算技术在航空航天飞行器的相关研究中的应用越来越广泛。相比于风洞试验等实验方法，数值模拟技术有着实验成本低，可重复性强，适用性强等优势。通过航空航天飞行器的数值模拟，可以得到航空航天飞行器的流动特性、热学特性以及电磁特性，再进一步根据数值模拟的结果，优化航空航天飞行器相关参数的设计。在数值模拟过程中各种数值计算方法（如有限元法、有限体积法、间断伽辽金法等）广泛应用于飞行器设计中，而高效精准的网格划分过程是这些数值计算方法的一个关键步骤。
网格划分是进行高精度数值模拟和数值计算的基础，网格质量的好坏直接影响着数值模拟和数值计算的精度。并且在数值模拟过程中，网格划分过程往往占据了整个数值模拟过程的大部分工作量，因而网格划分尤其是复杂模型的网格划分是数值模拟中十分重要的研究部分。在复杂模型及复杂条件下的数值模拟中对网格划分有着更高的要求。复杂模型中通常包含高曲率的部分或者局部微小的几何邻近特征。为保证离散精度和网格质量，在此类几何特征区域中线性四面体网格需要不断缩小网格尺寸来匹配其几何特征，进一步导致网格量级的不断上升，从而导致计算量的增大以及计算时间的增加。对于复杂外形，较差的网格质量还有可能导致计算过程的发散，进而导致计算失败。因此在超声速或多尺度等复杂条件下通过一些新的技术来提升网格质量是十分必要的，而曲网格是解决现有网格质量问题的一种较好方式。
现有曲网格划分方法基本包含直接法和间接法。直接法因计算复杂且适应性差等问题在数值计算中使用较少。间接法可以看作是基于线性网格划分的网格变形，是目前常用的方法。在间接法中，网格变形主要包含外推法、弹簧法和插值法等。外推法难以在复杂模型中取得应用；弹簧法理论较为复杂，计算过程耗时较长且稳定性较差。径向基函数法是一种通过插值点之间距离构建插值基的插值法，因其通用性好、变形网格质量好，近些年被广泛应用于动态网格以及网格变形等领域。但是径向基函数法往往需要进行大量的计算，求解过程非常耗时且不稳定。
发明内容
针对上述存在的问题或不足，为解决现有复杂模型网格划分过程中存在网格划分普适性较差，网格质量不佳的问题，本发明提供了一种复杂模型三维空间曲网格划分方法。通过提出一种新的控制点和形变点选取策略，在保证网格质量的同时减少了控制点和形变点的选取，减少了网格划分过程中大量的位置计算，优化了网格划分策略，减少了网格划分耗时以及降低了算力需求。
具体方案如下，一种复杂模型三维空间曲网格划分方法，包括以下步骤：
步骤1. 输入目标三维物理模型以及网格参数。
步骤2．将目标三维物理模型转换为三维空间几何信息，标记并存储计算域边界信息。
步骤3．依据步骤2中的三维空间几何信息，依次基于阵面推进法和三维约束德洛内（Delaunay）方法生成T个初始线性四面体单元，同时依据T个初始线性四面体单元的网格拓扑关系构造tetra（体）→face（面）→edge（边）→node（点）的层次结构存储网格信息。
步骤4．根据步骤3中生成的初始线性四面体单元的网格信息以及步骤2中的计算域边界信息，逐个检索步骤3生成的T个初始线性四面体单元。
如果第t个单元为边界单元则给该单元附上边界标记并计为第ѱ个边界单元，t∈｛1,2,3,…,T｝。遍历步骤3中T个初始线性四面体单元，可得到共计Ψ个边界单元，ѱ∈｛1,2,3,…,Ψ｝。
步骤5．根据步骤4中标记的Ψ个边界单元，逐个检索Ψ个边界单元，取出第ѱ个边界单元的边界面、边界边和边界点信息。
四面体单元中的边界面包含三条边界边以及三个顶点，将Ψ个边界单元中所有边界面上的边截断，生成边界面截断点。并将生成的边界面截断点投影至步骤2中的计算域边界得到投影点，且投影点和边界面截断点存在一一对应关系。
步骤6．根据步骤5中得到的Ψ个边界单元信息，逐个检索Ψ个边界单元的非边界面的面、边和点信息。
步骤7．根据步骤5中得到的Ψ个边界单元信息，逐个检索Ψ个边界单元的相邻单元，取出相邻单元中非相邻面的面、边和点的信息。
所述相邻单元：四面体单元中共用一个面的两个单元定义为相邻单元，非边界单元的相邻单元有四个，边界单元的相邻单元少于四个但至少有一个。
步骤8．遍历所有边界面截断点及其对应的投影点，计算边界面截断点与其对应投影点之间的夹角余弦值。
步骤9．将步骤5、6和7截断得到的全部截断点代入径向基函数计算式中，通过径向基函数计算式求得的值即为截断点形变量。
控制点为步骤5中边界面截断点，即全部边界面截断点的子集；先通过控制点和投影点计算权重；然后按先取边界面外层循环的全局编号顺序，再取边界边内层循环的局部编号的先后顺序遍历边界面截断点计算控制点位移量，重复点跳过。
步骤10．将步骤9计算出来的形变量代入截断点坐标以计算空间曲网格坐标（xc,yc,zc）。
步骤11．步骤10计算所得的空间曲网格点作为新增网格节点，***网格数据中，并更新网格拓扑关系。然后将新的空间曲网格数据用相应的网格数据文件格式输出，并将输出的网格数据文件导入可视化。
进一步的，所述步骤3的具体过程为：首先采用阵面推进法以计算域表面为初始阵面向内部推进生成W层的层次网格，其中W≥1，W可由用户根据实际情况选取（W越大空间曲网格精度越高，但计算量越大，W取值为W∈｛1,2,3,…,10｝）。
然后在剩余计算域中以当前层次网格表面为初始约束条件，采用三维约束德洛内（Delaunay）法生成初始线性四面体单元。
进一步的，所述步骤5中计算截断点位置时，比例参数λ取1，即为边的中点。
进一步的，所述步骤5中若要构造更高精度的曲网格，则在每个边界单元第一次截断的基础之上继续生成新的截断点。然后将分别求得的新一阶截断点投影到计算域边界得到相应的投影点并存入截断点数据类和投影点数据类，如此循环得到更多的高阶截断点和高阶投影点。
进一步的，所述步骤9中径向基函数中基函数为紧支型径向基函数。
进一步的，所述步骤9中径向基函数使用二次多项式。
进一步的，所述步骤11中网格数据文件格式输出采用CGNS格式。
进一步的，所述步骤6采用步骤5步中同样的方法求得非边界面截断点。
进一步的，所述步骤7采用步骤5步中同样的方法将每条边截断求得相邻单元截断点。
进一步的，所述步骤7中，对于多层空间曲网格则进一步循环检索相邻单元的相邻单元。
本发明提出的复杂模型三维空间曲网格划分方法通过对原有线性四面体截断生成曲面四面体单元来实现，能在不提升网格量级的情况下很好的满足复杂模型网格划分需求，初始线性网格基于阵面推进法和三维约束德洛内（Delaunay）法。空间曲网格点生成过程中仅选取步骤5中边界面截断点作为控制点，形变点也仅选取边界单元和边界单元相邻单元的截断点，大量减少了控制点和形变点的数量，从而大量减少了形变点即空间曲网格点相关的计算。相比于现有技术，步骤8创新性的提出了基于夹角余弦值的基函数选择方式；在保证网格质量的前提下，使用形式简单的低次多项式基函数的能大量减少形变量所需的计算。形变点计算减少与基函数形式简化使得网格划分的稳定性和速度得到了较大的提升。
综上所述，本发明通过提出一种新的控制点和形变点选取策略，在保证网格质量的同时减少了控制点和形变点的选取，减少了网格划分过程中大量的位置计算，优化了网格划分策略，减少了网格划分耗时以及降低了算力需求。在不提升网格量级的情况下很好的满足复杂模型的网格划分需求。
附图说明
图1为实施例目标模型三维空间曲网格示意图；
图2为钝锥（10°）线性/曲网格对比图；
图3为相邻单元剖面（二维）示意图；
图4为X-51A线性网格示意图；
图5为X-51A空间曲网格示意图；
图6为本发明流程图。
具体实施方式
下面结合附图和实施例来详细说明本发明的技术方案。
参照附图，一种复杂模型三维曲网格划分方法，包括以下步骤：
步骤1．输入目标三维物理模型以及网格参数。
根据三维物理模型（尺寸参数、公式、图纸、实物等）绘制基于ACIS内核的数字模型，通过程序读取网格尺寸，边界类型等网格参数设置。
步骤2．将目标三维物理模型转换为三维空间几何信息，标记并存储计算域边界信息。
建立Model类，将数字模型中三维空间几何信息导入模型Model类中，根据设置的边界类型建立边界标记（程序中为Model类的成员BCtype）。
步骤3．依据步骤2中的三维空间几何信息，依次基于阵面推进法和三维约束德洛内Delaunay方法生成T个初始线性四面体单元，同时依据T个初始线性四面体单元的网格拓扑关系构造tetra体→face面→edge边→node点的层次结构存储网格信息。
读取网格参数和Model类中的三维空间几何信息，首先采用阵面推进法由计算区域表面向内部推进生成多（球实施例中为两层）层网格单元，然后利用三维约束德洛内（Delaunay）法以当前阵面为约束条件生成剩余区域的网格。同时构造tetraArray类、faceArray类、edgeArray类、nodeArray类来存储单元信息。
步骤4．根据步骤3中生成的初始线性四面体单元的网格信息以及步骤2中的计算域边界信息，逐个检索步骤3生成的T个初始线性四面体单元。
如果第t个单元为边界单元则给该单元附上边界标记并计为第ѱ个边界单元，t∈｛1,2,3,…,T｝；遍历步骤3中T个初始线性四面体单元，可得到共计Ψ个边界单元，ѱ∈｛1,2,3,…,Ψ｝。
本实施例中：将边界标记BCtype***tetraArray类，并建立对应关系。再将tetraArray类中BCtype（程序中为tetraArray类的成员BCtype）的***faceArray类，edgeArray类，nodeArray类中，标记单元中的边界面、边界边和边界点。
步骤5．根据步骤4中标记的Ψ个边界单元，逐个检索Ψ个边界单元，取出第ѱ个边界单元的边界面、边界边和边界点信息；
四面体单元中的边界面包含三条边界边以及顶点P 1（x1,y1,z1）、P 2（x2,y2,z2）和P 3（x3,y3,z3）,x、y和z分别代表点在三维笛卡尔直角坐标系下三个维度的坐标值；第一条边对应顶点P 1和P 2，第二条边对应P 1和P 3，第三条边对应P 2和P 3。将Ψ个边界单元中所有边界面上的边（四面体单元面包含三条边）截断，生成边界面截断点。
第一条边界边（第二和第三条边界边同理）截断位置点的位置通过顶点P 1（x1,y1,z1）和P 2（x2,y2,z2）计算，其计算公式为：
Figure 532136DEST_PATH_IMAGE001
并将生成的边界面截断点P H （xn,1,yn,1,zn,1）（H代表边界面截断点标号）投影至步骤2中的计算域边界得到投影点P sub （xn,1,yn,1,zn,1）（sub代表投影点标号，且投影点和边界面截断点存在一一对应关系），λ为比例参数（大部分情况下为省略计算步骤可直接取λ=1即为边的中点，高曲率等特别的复杂模型中可以基于物理模型曲率、径向基函数来调整，调整范围为λ∈[0.5,2]。
每生成一个截断点，就将截断点P H （xn,1,yn,1,zn,1）存入截断点数据类P H （xN,1,yN,1,zN,1）中，并计数（P H 中n代表第n个截断点，最大计数值为N）。同时将所有投影至边界表面得到投影点P sub （xn,1,yn,1,zn,1）也存入投影点数据类P sub （xN,1,yN,1,zN,1）中，并计数（P sub 中n代表第n个投影点，最大计数值也为N）。
进一步的，如果要构造更高精度的曲网格，则在每个边界单元第一次截断的基础之上继续生成新的截断点。以点P H （xn,1,yn,1,zn,1）和P 1（x1,y1,z1），以及P H （xn,1,yn,1,zn,1）和P 2（x2,y2,z2）分别为端点继续根据
Figure 744943DEST_PATH_IMAGE002
生成截断点，此时（xi,yi,zi）为一次截断点P H （xn,1,yn,1,zn,1），（xj,yj,zj）为分别为顶点P 1（x1,y1,z1）和P 2（x2,y2,z2）的值。
然后将分别求得的截断点P H （xn,21,yn,21,zn,21）和P H （xn,22,yn,22,zn,22）投影到计算域边界得到P sub （xn,21,yn,21,zn,21）和P sub （xn,22,yn,22,zn,22）并存入相应的截断点数据类和投影点数据类。如此循环得到更多的高阶截断点和高阶投影点，一般情况下每条边截断一次得到一个一阶截断点或者截断两次得到三个二阶截断点就可以满足网格精度需求。虽然可以循环截取更多截断点来提高网格精度，但是截断次数越多带来的计算量越大，计算速度越慢且算力需求越大。因此在网格划分中不建议采用高于二阶截断点的计算方案。
本实施例为：在tetraArray类检索边界标记BCtype，当该单元为边界单元时检索单元对应的edgeArray类，nodeArray类。再进一步检索边的边界标记BCtype，检索到边界边时，通过该边界边对应的点坐标和公式求截断点P H （xn,1,yn,1,zn,1）和投影点P sub （xn,1,yn,1,zn,1）。求得截断点P H （xn,1,yn,1,zn,1）和投影点P sub （xn,1,yn,1,zn,1）后存入截断点数据类P H （x N,1,y N,1,z N,1）和P sub （x N,1,y N,1,z N,1）中，并计数（n最大计数值为N）。
步骤6．根据步骤5中得到的Ψ个边界单元信息，逐个检索Ψ个边界单元的非边界面的面、边和点信息；采用步骤5步中同样的方法求得非边界面截断点P inr （xm,1,ym,1,zm,1），存入截断点数据类P inr （xM,1,yM,1,zM,1），并计数，inr代表非边界面截断点标号，P inr 中m最大计数值记为M。
步骤7．根据步骤5中得到的Ψ个边界单元信息，逐个检索Ψ个边界单元的相邻单元，取出相邻单元中非相邻面的面、边和点的信息，如果需要多层空间曲网格则进一步循环检索相邻单元的相邻单元。
所述相邻单元：四面体单元中共用一个面的两个单元定义为相邻单元，非边界单元的相邻单元有四个，边界单元的相邻单元少于四个但至少有一个。
采用步骤5步中同样的方法将每条边截断求得相邻单元截断点P adj （xk,1,yk,1,zk,1）存入截断点数据类P adj （xK,1,yK,1,zK,1）中，并计数，P adj 中k最大计数值为K，adj代表相邻单元截断点标号。
本实施例中，检索faceArray类的类成员tetra[A,B],A和B即为面所属的两个单元全局编号。再根据全局编号A或者B，在tetraArray类中检索相邻单元，取出非相邻面的边和点的信息。
步骤8．遍历边界面截断点数据类P H （xN,1,yN,1,zN,1）和投影点数据类P sub （xN,1,yN,1,zN,1）中所有点，计算边界面截断点P H （xn,1,yn,1,zn,1）与其对应投影点P sub （xn,1,yn,1,zn,1）之间的夹角余弦值，记为cos(θ) H-sub ，存入nodeArray类中每个点对应的位置。
Figure 811250DEST_PATH_IMAGE003
其中xn,H,yn,H,zn,H为第n个边界面截断点的三维坐标值，代表xn,sub ,yn,sub ,zn,sub 为第n个投影点的三维坐标值。
步骤9．将步骤5、6和7截断得到的全部截断点代入径向基函数计算式中，通过径向基函数计算式求得的f(x,y,z)的值即为截断点形变量；存入临时变量mesh-distancecurved（总数N+M+K）中。径向基函数计算式为：
Figure 411996DEST_PATH_IMAGE004
其中α为系数，
Figure 547442DEST_PATH_IMAGE005
为权重，I N 为控制点总数且I N =N；/>
Figure 157415DEST_PATH_IMAGE006
为径向基函数，
Figure 920971DEST_PATH_IMAGE007
表示全部截断点到对应控制点的欧式距离，控制点/>
Figure 82831DEST_PATH_IMAGE008
为步骤5中边界面截断点P H （xN,1,yN,1,zN,1），即全部边界面截断点的子集；
在实施例中，基于截断点（即空间曲网格点）的最大夹角余弦值，本发明所使用的径向基函数计算式为：
Figure 299049DEST_PATH_IMAGE009
径向基函数在大部分情况下仅使用二次多项式，使用形式简单的低次多项式基函数能大量减少形变量所需的计算；仅在复杂模型出现高曲率等复杂情况则采用高次多项式来满足网格质量要求。
先通过控制点（边界面截断点）和投影点计算权重
Figure 712713DEST_PATH_IMAGE005
，计算权重时/>
Figure 206142DEST_PATH_IMAGE007
中X代表边界面截断点（即为步骤5中边界面截断点P H （xN,1,yN,1,zN,1），所以总数为I N =N；计算权重时f (x,y,z)为已知量，代表控制点位移量，f (x,y,z)取同一边界面上截断点到投影点的欧式距离的最小值，以避免有平凡解的情况出现；
f (x,y,z)=min{║X b -X sub ║}
其中一阶截断点情况下X b ,b∈｛1,2,3｝为同一边界面上的三条边界边对应的三个截断点，X sub 为截断点相对应的投影点的坐标值；然后按先取边界面外层循环的全局编号顺序，再取边界边内层循环的局部编号的先后顺序遍历边界面截断点计算控制点位移量，重复点跳过；
求得权重
Figure 945428DEST_PATH_IMAGE005
后代入截断点计算f (x,y,z)，此时/>
Figure 74708DEST_PATH_IMAGE005
为已知量，f (x,y,z)为待求量；计算网格点位移时/>
Figure 636270DEST_PATH_IMAGE010
中X=X j ，X j 代表步骤5中边界面截断点、步骤6中非边界面截断点和步骤7中相邻单元截断点求得的全部截断点，截断点总数I SUM =N+M+K；
步骤10．将步骤9计算出来的形变量f (x,y,z)代入截断点坐标（xJ,1,yJ,1,zJ,1）以计算空间曲网格坐标（xc,yc,zc），其计算式为：
（xc,yc,zc）=（xJ,1,yJ,1,zJ,1）+f (x,y,z),J∈{N},{M},{K}
遍历临时变量mesh-distancecurved取出每个截断点对应的形变量，与截断点初始位置相加求得形变后的截断点（相加后的截断点位置即为空间曲网格点）位置，即为空间曲网格点坐标。再将求得空间曲网格点坐标存入截断点类，代替原有截断点坐标。
步骤11．步骤10计算所得的高阶网格点（xc,yc,zc）作为新增网格节点，即从截断点类读取步骤10计算所得的空间曲网格点坐标，作为新增网格节点node5至node10，***网格数据中，并更新网格拓扑关系；然后将新的空间曲网格数据用相应的网格数据文件格式（如CGNS格式）输出，并将输出的网格数据文件导入可视化。
本实施例采用3个目标物体做出示例以验证本发明的有效性；各实施例最终网格划分的空间曲网格结果如图1和图2的右图，以及图5所示。
图1是以球体为目标物体，左图空间曲网格全局表面视图，中间图为空间曲网格剖面图，右图为局部放大图。
图2是钝锥（10°）空间曲网格结果图。左图为步骤三所得初始线性四面体网格，右图为完成后的空间曲网格示意图。
图4和5是X-51A飞行器的结果图。图4为步骤三所得初始线性四面体网格，图5为完成后的空间曲网格示意图。
从图2中左图和右图的对比，图4和图5的对比可以看出，在保证网格质量和基本拓扑结构不变的情况下，本发明所提出的空间曲网格划分方法实现了高精度的空格曲网格。现有曲网格技术基本上只实现了物面曲网格，虽然物面曲网格提升了边界处的计算精度，数值计算精度往往受限制于内部网格的精度。本发明采取了全新的控制点和形变点（即空间曲网格点）选取策略，将曲网格拓展至空间内层网格中实现复杂模型的高精度空间曲网格划分。综上可见本发明切实有效。
综上可见，本发明通过提出一种新的控制点和形变点选取策略，在不提升网格量级的情况下很好的满足复杂模型的网格划分需求。相比现有技术：1. 本发明实现了空间曲网格划分，现有技术基本为物面曲面网格；2. 采用新的控制点和形变点选取策略以及计算方法，在不提升网格量级的情况下很好的满足复杂模型的网格划分需求；3. 径向基函数根据夹角余弦值选取。本发明在保证网格质量的同时减少了控制点和形变点的选取，减少了网格划分过程中大量的位置计算，优化了网格划分策略，减少了网格划分耗时以及降低了算力需求。形变点计算减少与基函数形式简化使得网格划分的稳定性和速度得到了较大的提升。

Claims (5)

1.一种复杂模型三维空间曲网格划分方法，其特征在于，包括以下步骤：
步骤1.输入目标三维物理模型以及网格参数；
步骤2.将目标三维物理模型转换为三维空间几何信息，标记并存储计算域边界信息；
步骤3.依据步骤2中的三维空间几何信息，依次基于阵面推进法和三维约束德洛内Delaunay方法生成T个初始线性四面体单元，同时依据T个初始线性四面体单元的网格拓扑关系构造tetra体→face面→edge边→node点的层次结构存储网格信息；
步骤4.根据步骤3中生成的初始线性四面体单元的网格信息以及步骤2中的计算域边界信息，逐个检索步骤3生成的T个初始线性四面体单元；
如果第t个单元为边界单元则给该单元附上边界标记并计为第ψ个边界单元，t∈{1,2,3,…,T}；遍历步骤3中T个初始线性四面体单元，可得到共计Ψ个边界单元，ψ∈{1,2,3,…,Ψ}；
步骤5.根据步骤4中标记的Ψ个边界单元，逐个检索Ψ个边界单元，取出第ψ个边界单元的边界面、边界边和边界点信息；
边界面包含三条边界边以及顶点P1(x1,y1,z1)、P2(x2,y2,z2)和P3(x3,y3,z3)，x、y和z分别代表点在三维笛卡尔直角坐标系下三个维度的坐标值；第一条边对应顶点P1和P2，第二条边对应P1和P3，第三条边对应P2和P3；将Ψ个边界单元中所有边界面上的边截断，生成边界面截断点，四面体单元面包含三条边；
第一条边界边截断位置点的位置通过顶点P1(x1,y1,z1)和P2(x2,y2,z2)计算，第二和第三条边界边同理，其计算公式为：
Figure QLYQS_1
并将生成的边界面截断点PH(xn,1,yn,1,zn,1)投影至步骤2中的计算域边界得到投影点Psub(xn,1,yn,1,zn,1)，其中H代表边界面截断点标号，sub代表投影点标号，且投影点和边界面截断点存在一一对应关系；λ为比例参数，基于物理模型曲率、径向基函数来调整，调整范围为λ∈[0.5,2]；
每生成一个截断点，就将截断点PH(xn,1,yn,1,zn,1)存入截断点数据类PH(xN,1,yN,1,zN,1)中，并计数，PH中n代表第n个截断点，最大计数值为N；同时将所有投影至边界表面得到投影点Psub(xn,1,yn,1,zn,1)也存入投影点数据类Psub(xN,1,yN,1,zN,1)中，并计数，Psub中n代表第n个投影点，最大计数值也为N；
步骤6.根据步骤5中得到的Ψ个边界单元信息，逐个检索Ψ个边界单元的非边界面的面、边和点信息；采用步骤5步中同样的方法求得非边界面截断点Pinr(xm,1,ym,1,zm,1)，存入截断点数据类Pinr(xM,1,yM,1,zM,1)，并计数，inr代表非边界面截断点标号，Pinr中m最大计数值记为M；
步骤7.根据步骤5中得到的Ψ个边界单元信息，逐个检索Ψ个边界单元的相邻单元，取出相邻单元中非相邻面的面、边和点的信息，如果需要多层空间曲网格则进一步循环检索相邻单元的相邻单元；
所述相邻单元：四面体单元中共用一个面的两个单元定义为相邻单元，非边界单元的相邻单元有四个，边界单元的相邻单元少于四个但至少有一个；
采用步骤5中同样的方法将每条边截断求得相邻单元截断点Padj(xk,1,yk,1,zk,1)存入截断点数据类Padj(xK,1,yK,1,zK,1)中，并计数，Padj中k最大计数值为K，adj代表相邻单元截断点标号；
步骤8.遍历边界面截断点数据类PH(xN,1,yN,1,zN,1)和投影点数据类Psub(xN,1,yN,1,zN,1)中所有点，计算边界面截断点PH(xn,1,yn,1,zn,1)与其对应投影点Psub(xn,1,yn,1,zn,1)之间的夹角余弦值，记为cos(θ)H-sub；
Figure QLYQS_2
其中xn,H,yn,H,zn,H为第n个边界面截断点的三维坐标值，代表xn,sub,yn,sub,zn,sub为第n个投影点的三维坐标值；
步骤9.将步骤5、6和7截断得到的全部截断点代入径向基函数计算式中，通过径向基函数计算式求得的f(x,y,z)的值即为截断点形变量；
径向基函数计算式为：
Figure QLYQS_3
其中α为系数，
Figure QLYQS_4
为权重，IN为控制点总数且IN＝N；/>
Figure QLYQS_5
为径向基函数，/>
Figure QLYQS_6
表示全部截断点到对应控制点的欧式距离，控制点/>
Figure QLYQS_7
为步骤5中边界面截断点PH(xN,1,yN,1,zN,1)，即全部边界面截断点的子集；
先通过控制点和投影点计算权重
Figure QLYQS_8
计算权重时/>
Figure QLYQS_9
中X代表边界面截断点，所以总数为IN＝N；计算权重时f(x,y,z)为已知量，代表控制点位移量，f(x,y,z)取同一边界面上截断点到投影点的欧式距离的最小值；
f(x,y,z)＝min{║Xb-Xsub║}
其中一阶截断点情况下Xb,b∈{1,2,3}为同一边界面上的三条边界边对应的三个截断点，Xsub为截断点相对应的投影点的坐标值；按先取边界面外层循环的全局编号顺序，再取边界边内层循环的局部编号的先后顺序遍历边界面截断点计算控制点位移量，重复点跳过；
求得权重
Figure QLYQS_10
后代入截断点计算f(x,y,z)，此时/>
Figure QLYQS_11
为已知量，f(x,y,z)为待求量；计算网格点位移时/>
Figure QLYQS_12
中X＝Xj，Xj代表步骤5中边界面截断点、步骤6中非边界面截断点和步骤7中相邻单元截断点求得的全部截断点，截断点总数ISUM＝N+M+K；
步骤10.将步骤9计算出来的形变量f(x,y,z)代入截断点坐标(xJ,1,yJ,1,zJ,1)以计算空间曲网格坐标(xc,yc,zc)，其计算式为：
(xc,yc,zc)＝(xJ,1,yJ,1,zJ,1)+f(x,y,z),J∈{N},{M},{K}
步骤11.步骤10计算所得的高阶网格点(xc,yc,zc)作为新增网格节点，***网格数据中，并更新网格拓扑关系；然后将新的空间曲网格数据用相应的网格数据文件格式输出，并将输出的网格数据文件导入可视化。
2.如权利要求1所述复杂模型三维空间曲网格划分方法，其特征在于，所述步骤3中：
首先，采用阵面推进法以计算域表面为初始阵面向内部推进生成W层的层次网格，其中W≥1，W∈{1,2,3,…,10}；
然后，在剩余计算域中以当前层次网格表面为初始约束条件，采用三维约束德洛内Delaunay法生成初始线性四面体单元。
3.如权利要求1所述复杂模型三维空间曲网格划分方法，其特征在于：所述步骤5中取λ＝1，即为边的中点。
4.如权利要求1所述复杂模型三维空间曲网格划分方法，其特征在于：
所述步骤5中：
如果要构造高精度的曲网格，则在每个边界单元第一次截断的基础之上继续生成新的截断点；以点PH(xn,1,yn,1,zn,1)和P1(x1,y1,z1)，以及PH(xn,1,yn,1,zn,1)和P2(x2,y2,z2)分别为端点继续根据
Figure QLYQS_13
生成截断点，此时(xi,yi,zi)为一次截断点PH(xn,1,yn,1,zn,1)，(xj,yj,zj)分别为顶点P1(x1,y1,z1)和P2(x2,y2,z2)的值；
然后将分别求得的截断点PH(xn,21,yn,21,zn,21)和PH(xn,22,yn,22,zn,22)投影到计算域边界得到Psub(xn,21,yn,21,zn,21)和Psub(xn,22,yn,22,zn,22)并存入相应的截断点数据类和投影点数据类，如此循环得到多的高阶截断点和高阶投影点。
5.如权利要求1所述复杂模型三维空间曲网格划分方法，其特征在于：所述步骤11中网格数据文件格式输出采用CGNS格式。
CN202111637505.5A 2021-12-30 2021-12-30 一种复杂模型三维空间曲网格划分方法 Active CN114510775B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111637505.5A CN114510775B (zh) 2021-12-30 2021-12-30 一种复杂模型三维空间曲网格划分方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111637505.5A CN114510775B (zh) 2021-12-30 2021-12-30 一种复杂模型三维空间曲网格划分方法

Publications (2)

Publication Number Publication Date
CN114510775A CN114510775A (zh) 2022-05-17
CN114510775B true CN114510775B (zh) 2023-06-27

Family

ID=81548505

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111637505.5A Active CN114510775B (zh) 2021-12-30 2021-12-30 一种复杂模型三维空间曲网格划分方法

Country Status (1)

Country Link
CN (1) CN114510775B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115186074B (zh) * 2022-06-24 2023-07-21 深圳市规划和自然资源数据管理中心 基于Meta分析的模拟土壤pH值空间分布格局的方法
CN116401916B (zh) * 2023-03-20 2024-01-26 北京云境智仿信息技术有限公司 高质量三维体网格的生成方法、装置、介质及设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109636913A (zh) * 2018-12-04 2019-04-16 山东理工大学 基于Delaunay剖分的三角网格增量拓扑拼接方法
CN112015735A (zh) * 2020-08-20 2020-12-01 西安数峰信息科技有限责任公司 一种非结构化网格的数据存储结构及数据存储方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10186079B2 (en) * 2012-05-14 2019-01-22 Autodesk, Inc. Adaptively joining meshes

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109636913A (zh) * 2018-12-04 2019-04-16 山东理工大学 基于Delaunay剖分的三角网格增量拓扑拼接方法
CN112015735A (zh) * 2020-08-20 2020-12-01 西安数峰信息科技有限责任公司 一种非结构化网格的数据存储结构及数据存储方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
用逐点***法生成Delaunay四面体自适应网络;骆冠勇 等;计算力学学报;第24卷(第6期);917-922 *

Also Published As

Publication number Publication date
CN114510775A (zh) 2022-05-17

Similar Documents

Publication Publication Date Title
CN114510775B (zh) 一种复杂模型三维空间曲网格划分方法
CN102509339B (zh) 一种带纹理约束的三维模型顶点聚类简化方法
WO2021203711A1 (zh) 一种基于几何重建模型的等几何分析方法
Shu High order numerical methods for time dependent Hamilton-Jacobi equations
CN100454341C (zh) 一种基于分形层次树的过程式地形快速绘制方法
CN114492250A (zh) 基于递归分解的曲面网格生成方法及***、计算机设备
Pan et al. Subdivision-based isogeometric analysis for second order partial differential equations on surfaces
Natarajan et al. Simplification of three-dimensional density maps
Jekeli Spline representations of functions on a sphere for geopotential modeling
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
Seshadri et al. Supporting multi-point fan design with dimension reduction
CN113342999A (zh) 一种基于多层跳序树结构的变分辨率点云简化方法
CN110765298B (zh) 矢量数据几何属性解耦的瓦片编码方法
Langbein et al. An efficient point location method for visualization in large unstructured grids.
CN114510677B (zh) 基于间断有限元的中子输运方程处理方法、计算机程序产品
Crampton et al. Detecting and approximating fault lines from randomly scattered data
WO2021249374A1 (zh) 特征线追踪方法、堆芯中子物理计算方法和装置
CN113111612A (zh) 一种基于自适应空间剖分的离散点云重复点快速查找方法
Hou et al. Knot optimization for biharmonic b-splines on manifold triangle meshes
Brandts et al. A geometric toolbox for tetrahedral finite element partitions
Mittal et al. Mixed-Order Meshes through rp-adaptivity for Surface Fitting to Implicit Geometries
Wang et al. An automated mesh generation algorithm for curved surfaces of ship longitudinal structures
CN114004175B (zh) 一种快速查找全域壁面距离和无量纲壁面距离的方法

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