三维空间非结构网格生成方法研究

三维空间非结构网格生成方法研究
郭丹;肖素梅;庞宇飞;戴峰
【摘要】With the rapid development of Computational Fluid Dynamics, the problems to be solved with CFD are becom-ing more and more complex. However, the traditional unified body-fitted grids structure has been unable to meet the needs of development of CFD. Unstructured mesh appears and attracts the CFD researchers'attention with new features by breaking the habit of structured thinking. And the unstructured grid with boundary layer is always a difficult problem in unstructured grid generation. The paper studies the grid generating method for unstructured tetrahedral, at the same time, combined with the technological advantages of Spider software platform structure grid parameters of boundary layer pro-pulsion, based on it, unstructured mesh generation module UGCS development is done. The large number of examples of grid generation results, mainly grid quality and numerical calculation, fully verify the reliability and the robustness of algorithm.%随着计算流体力学领域待解决问题复杂程度的不断提高,传统的统一贴体结构网格已不能很好地满足针对复杂外形的高精度网格生成需求,而非结构网格以其独特优势受到CFD工作者的普遍关注.带有附面层的非结构网格是非结构网格生成的难点.进行了非结构四面体网格的生成方法研究,同时结合Spider软件平台中结构网格参数化附面层推进的技术优势,进行了基于Spider软件平台中非结构网格生成模块"UGCS"的开发.通过分析大量网格生成实例中网格质量和数值计算结果,验证了算法的可靠性与鲁棒性.
【期刊名称】《计算机工程与应用》
【年(卷),期】2017(053)022
【总页数】6页(P45-49,54)
【关键词】复杂外形;非结构网格;附面层网格;网格生成自动化
【作者】郭丹;肖素梅;庞宇飞;戴峰
【作者单位】西南科技大学制造科学与工程学院,四川绵阳 621010;西南科技大学制造科学与工程学院,四川绵阳 621010;中国空气动力研究与发展中心,四川绵阳621000;西南科技大学制造科学与工程学院,四川绵阳 621010
【正文语种】中文
【中图分类】TK05;TH6
由于非结构网格没有对应于ξ,η,ς的坐标线，这就使得我们能够充分利用网格的灵活性，来控制网格，使其放在合适的位置。

某种程度上来说，构造一个非结构网格可被视为一种艺术，因为可以随意构造网格的形状并且将他们放在物理平面中的任意位置上。

必须利用计算机程序来实现网格生成的自动化[1]。

尽管非结构网格已经在结构力学的有限元计算中应用很多年了，但是在CFD领域中它们还相对较新[2]。

基于Spider软件平台，利用结构网格的参数化附面层推进功能[3]来获取非结构的附面层三棱柱网格，同时利用文中的三维空间四面体网格生成算法，对外场非结构网格进行了填充。

给出了某翼型的非结构网格生成实例，通过生成网格的数值模拟和数据比对证明了算法的鲁棒性好，通过做一套可计算网格耗时比对证明自动化程度高。

文章将结构网格转化为非结构时的方法进行曲面三角化处理[4]。

采用这种方法，
是基于下面三方面考虑：首先，基于软件开发的刚性需求，为了能够快速地看到非结构网格生成效果，验证三维四面体网格的生成算法鲁棒性高，且生成网格能够用于数值计算；其次，表面结构网格剖分为非结构网格，简单又直接，而且效率很高；最后，也是非常重要的一点，采用这种方法可以充分利用结构网格中附面层推进的特有优势，为后期生成带有附面层的非结构网格建立基础。

将结构网格中的Domain数据按照u,v指标进行小四边形网格单元的对角连接，
即可完成结构网格向非结构网格的转化，但有一点需要注意，就是每个小四边形分解成的两个三角形，其角点必须按顺序连接[5]。

2.2.1 算法说明
具体算法如下：
（1）首先定义一个能够覆盖所有将要插入点的凸壳，原始凸壳一般是由8个顶点构成的立方体。

以立方体的6个面构成6个四面体，建立原始Delaunay数据结构。

（2）在内部引入一个新的点。

（3）依据Delaunay原则，给要删除的四面体做标记。

（4）建立空洞边界。

（5）在剩余四面体中查找被标记的相邻四面体。

（6）找出新四面体的顶点。

（7）建立现有四面体的关系。

（8）更新数据结构，处理新旧四面体及对应关系。

（9）再引入下一个点，重复同样的操作即可。

2.2.2 数据结构
该算法数据结构的设计采用C++编程语言，根据面向对象思想，设计了一个四面
体派生类UGTetCell，它继承了非结构中的基类UGCell。

整个UGTetCell的数据结构包括一个整型数组，即四面体的四个顶点int ver[4]；一个list容器，用来存放当前四面体的四个邻居list＜UGTetCell*＞tet_neighbor；一个四面体外接球心的点坐标指针UGPoint tet_center及其半径double tet_r。

整个四面体集合的定义：list＜UGTetCell*＞ug_tets。

2.2.3 算法执行
为了更好地理解，现对算法做出详细的说明：
（1）首先需要一个把物面网格包含进去的大凸壳，本文选择立方体的8个顶点作为凸壳的初始顶点，这8个点的坐标确定了背景网格区域的大小，通常选择至少10倍模型尺寸。

（2）假设物面网格点是已知的，能够按照一定顺序输入且不考虑流场内的插点问题。

（3）根据外接球准则，在每个四面体的外接球内部，一旦有新点插入，该四面体将被标记以便后续删除。

即：
其中x(tet0,i)表示tet0的第i个邻居的球心坐标，r(tet0,i)表示tet0的第i个邻居的外接球半径。

若满足条件则继续查找下一个，若都不满足则停止查找。

这样大大提高了查找效率。

（4）定义一个储存被标记四面体的临时容器变量，这些被标记的四面体的顶点构成了包围新点的凸壳。

（5）查找被标记四面体的邻居，修改它们的邻居变量tet_neighbor，将被标记的四面体单元删除掉，同时该位置用新四面体元素代替。

（6）一个新四面体由新顶点和被标记四面体的某三个顶点组成。

而这三个点的正确组合可以在计算第5步的同时得到。

（7）首先确定新四面体之间的关系，然后建立新四面体与在第5步已标记的邻居
之间的关系。

（8）对每个四面体的顶点坐标ver，球心坐标tet_center和半径tet_r，每个四面体的邻居单元tet_neighbor进行更新修改。

2.2.4 正确性检测
由实际编程经验可知，上述算法第3步中外接球准则的判断是最棘手的部分。

本
文采用双精度数据类型处理三维问题，因此在尺寸极小的情况下容易出现误差，如果仅把式（2）作为判断条件，在多次插入点后可能会出现错误，将导致严重后果。

本文在程序运行时采取了三项措施以解决此问题。

一，按照特定顺序储存四面体的四个顶点，以使体积结果均为正。

二，检测生成四面体时拓扑的有效性。

三，检测每个单元的邻居，显然，两个互为邻居的四面体在彼此的tet_neighbor中都保存了指向对方的指针。

在程序设计时，外接球准则应采用下式判断：
如果满足上面的条件，一个坐标为xnew的新点就落在了一个半径为r，球心坐标为xε的外接球之内。

一个最慢的方法就是全局搜索，即逐个搜索判断，但若采用
前文建立的list容器数据结构，查找效率就得到了显著的提高。

假设找到一个外接球包含新点的四面体tet0，则以tet0为起点进行查找。

将式（1）应用于tet0的邻居，即：
其中，ε是一个极小的值，根据不同情况取值。

经验证明在一定范围内ε的取值与四面体的质量呈正比关系。

当上式介于0和ε之间时，此时先不插入新点，而是
先将其放到新点列表的末尾，再插下个新点，因为四面体拓扑结构发生了改变，所以几乎不会出现临界情况。

实践证明在ε值选取适当时，可依次插入全部物面节点。

本文取插点前物面节点间最小距离的1/10作为全局误差容忍值tolerance，即
ε[6]。

采用上述检测方法确保了本方法的鲁棒性。

实践测试证明，通过上述检测手段，进
行了复杂外形进行外场填充[7]，并能输出通用的CGNS文件进行定性计算。

当前的Spider软件开发中还未加入针对非结构网格检查的可视化显示，故将数据保存为CGNS数据文件在Pointwise软件中进行了与Spider相同位置的切面检查。

如图1所示，以带有吊舱的F6模型进行的非结构外场填充，可以看出网格质量非常好。

目前，非结构的附面层网格生成技术一直是研究的难点[8]，除了Pointwise的T-Rex技术，将各向异性四面体融合来获得附面层的三棱柱单元比较成熟外，还有一种层推进技术[9]，这种方法通过数学方法进行判断处理较为频繁，而且容易出现
推进不成功的情况，导致程序崩溃。

文中结合Spider中框架线推进的技术优势，将推进得到的结构网格块，按照逐“层”的方式，沿着法向量推进的方向（l方向），将每层中的上下两面l和
l+1(l≥0)，按照j和k指标进行面剖分，方法与2.1节类似，将得到的角点标号按照右手定则排列，得到三棱柱的六个角点标号。

另外一种得到附面层单元的方式，是将附面层剖分为六面体单元以及与之紧相贴的金字塔单元[10]，这种体单元的剖分方式与前面所述的三棱柱剖分方式大同小异，主要是金字塔单元的构造需要计算法向量方向以及推进高度，由于本文主要以三棱柱单元的附面层作为计算用例，因此此种方式的网格单元剖分在此不作详述。

对于推进的附面层网格，首层推进高度可以达到0.001 mm级别，这种方法既能
保证程序稳定性，又能保证附面层内单元的网格质量[11]。

而且，结构化网格无论采用哪种生成方法，都可以使用此种转化[12]。

使用这种方法，既能保证网格质量，又能快速完成整个非结构网格的生成。

以747的头部为例，图2给出了带有附面
层的完整网格示例。

3.1.1 控制方程
三维N-S方程可以写成如下的无量纲形式：
其中，V为单元体积，n为单元表面外法向，N为单元的表面总数(k=1,2,…,N)，Sk为k的表面面积，H 和 Hv分别为表面法向无黏和黏性通量。

对于无黏通量进行分裂：
其中：
其中，U为守恒变量；E,F,G为无黏通量；Ev,Fv,Gv为黏性通量。

3.1.2 空间离散
对式（4）在控制单元内进行积分，可得：
其中∇U为守恒变量的梯度，而ϕi和ϕj为左右两侧的限制系数（显然0＜ϕ＜1）。

限制器采用了改进后的Barth-Jesperen限制器。

对于黏性通量采用中心格式，即：
对于湍流流动，引入了多种湍流模型，如一方程SA模型，两方程k-ω模型等，
湍流模型方程的计算与RANS方程相似，不同之处为对流项采用了一阶迎风离散。

3.1.3 隐式时间离散
对于式（5）进行隐式处理，可得：
式（9）可以进一步写为：
其中，A为Jacobian矩阵，并进一步整理为：
采用通量分裂法，可以将左侧的质量矩阵退化为一个对角阵。

为了改进收敛效率，这里采用了块LU-SGS（BLU-SGS）方法。

该算例在Windows 7旗舰版64位系统，16 GB内存，4核处理器上进行并行计算，计算条件是马赫数（Ma）为0.839 5，雷诺数（Re）为1.172E+7，攻角为3.06°。

采用Roe格式，LU-SGS时间推进，分别采用了SA和SST两种湍流模型，利用所生成的网格对两种湍流模型的差异进行了比较研究[13]。

在进行非结构网格生成时，将机翼前缘网格沿流向和前缘翼稍处进行了加密，生成三棱柱/四面体混
合网格，节点总数为848 002，单元面总数为5 473 552，单元总数为2 367 698，
附面层推进的第一层高度为5.0E－6。

图3为计算用的网格。

3.3.1 计算结果分析
计算步数为50 000步，从图4中可知，在迭代20 000步后，两种湍流模型的计算残差基本不再下降。

3.3.2 流场结果分析
分别采用两种湍流模型对文中的研究方法生成的网格进行了表面压力和密度等值线的绘制，从图5中（a）和（b）可看出，二者分布形态基本相同，均可捕捉到明
显的“λ”激波，并且通过比较图5中（c）和（d）显示的密度等值线可知，SA
模型计算的机翼梢部的分离区域更大。

3.3.3 物面压力分析
各站位处x起始值和终止值如表1所示。

从图6中不同站位处的截面压力分布可知，两种湍流模型结果与实验值相符，整
体来说，SA的计算结果更接近实验值，除了不同湍流模型预测的激波位置有略微
偏差，其余位置的压力系数几乎重合。

对于同一网格，采用两种不同模型进行分析计算，且计算结果基本都与实验值相吻合，由此也说明了通过本文方法生成的网格质量很好[14]。

文章采用所提出的空间四面体网格生成方法对747头部模型以及某型号机翼进行
了远场四面体网格填充，同时利用了Spider的参数化框架线推进功能，完成非结构的附面层网格生成[15]。

通过将Spider生成的非结构网格导入至Pointwise中，检查网格质量，效果很好；对网格质量的分析和典型位置的气动力数值模拟计算，得出了下述结论：
（1）文章对Delaunay算法重新进行了数据结构设计，使得远场的四面体网格生成耗时更少，效率更高。

（2）采用两种不同的湍流模型，对同一翼型网格进行了数值计算分析，结果显示，
两种模型都能很好地与实验值相吻合。

由此也说明该算法有较好的健壮性，生成的网格质量令人满意。

【相关文献】
[1]Lohner R，Aubry R.Recent development in tet-meshing[C]//Proceedings of the 16th International Meshing Roundtable.Seattl：Spring Publishing，2008.
[2]John D，Anderson J putational fluid dynamics：the basics with
applications[M].Beijing：Tsinghua University Press，2002.
[3]蔡云龙.基于参数化的气动计算网格自动生成技术研究[D].四川绵阳：西南科技大学，2014.
[4]张来平，杨永键，呙超，等.三维复杂外形的非结构网格自动生成技术与应用[J].计算物理，1999（5）：552-558.
[5]朱培烨.空间曲面的非结构网格生成[J].航空学报，2001，22（2）：151-154.
[6]段旭鹏.复杂外形动态混合网格生成技术及非定常计算方法研究[D].绵阳：中国空气动力研究与发展中心，2008.
[7]孙少鹏.非结构网格生成技术的研究[J].空气动力学学报，1996，14（1）：19-25.
[8]陈建军.非结构化网格生成及其并行化的若干问题研究[D].杭州：浙江大学，2006.
[9]赵钟.复杂外形高雷诺数数值模拟的混合网格生成与多重网格法研究[D].绵阳：中国空气动力研究与发展中心，2011.
[10]刘学强，伍贻兆.三维复杂外形非结构网格的生成及应用[J].南京航空航天大学学报，2000，32（5）：533-538.
[11]Smith R J，Johnston L J.Automatic grid generation and flow solution for complex geometries[J].AIAA Journal，1996，34（6）：1120-1124.
[12]刘晶.结构与非结构网格的生成、转化及应用[D].南京：南京理工大学，2006.
[13]肖中云.旋翼流场数值模拟方法研究[D].绵阳：中国空气动力研究与发展中心，2007.
[14]牟斌.流动控制数值模拟研究[D].绵阳：中国空气动力研究与发展中心，2006.
[15]Pirzadeh S.Three-dimensional unstructured viscous grids by the advancing-layers method[J].AIAA Journal，1996，34（1）：43-49.。