黏性边界层网格自动生成1)

甘洋科 刘剑飞

(北京大学工学院力学与工程科学系，北京100871)

黏性边界层网格自动生成1)

甘洋科 刘剑飞2)

(北京大学工学院力学与工程科学系，北京100871)

高雷诺数黏性流动在壁面附近存在边界层，在计算模拟中自动生成可靠且有效的计算网格仍然是计算流体力学存在的瓶颈.三棱柱/四面体混合网格技术在一定程度上缓解了这个困难.然而，对于复杂外形的情况，在边界层内自动高效生成高质量的三棱柱单元仍然十分困难.常用的层推进法在凹凸区域及角点处生成的边界层网格单元质量较差，边界层网格最外层尺寸不均匀.针对这些问题，发展了一种黏性边界层网格自动生成方法，通过顶点周围边的二面角识别物面网格特征确定多生长方向，预估并调整生长高度处理相交情况.同时提出一种三维前沿尺寸调节方式，提高了边界层网格单元的正交性，保证了边界层网格与远场网格尺寸的光滑过渡.通过ONERA M6翼型以及带发动机短舱的DLR-F6翼身组合体等外形的网格生成实例及绕流数值模拟，将计算值与标准实验值进行对比，结果表明：该方法能够自动高效地生成满足数值计算需求的混合网格.

边界层网格，生长法向，特征边，相交，尺寸调节

引言

近年来，计算流体力学成为学术界研究的热点，在航空航天、汽车、热能动力等工程领域都有着重要的地位，广泛用于设计、分析和优化.在计算流体力学模拟中，需要借助计算机数值求解物理量的控制方程.当前主要采用的计算方法，包括有限差分法、有限体积法、有限元法等，都需要对计算区域进行网格划分，进而在网格单元节点上离散求解偏微分方程.30多年来，尽管计算流体力学取得了很大的成就，网格生成依然是需要解决的关键问题之一.大批商业网格生成软件，包括Gambit，Gridgen以及ANSYS ICEM-CFD等，虽然具有功能全面、使用方便、技术服务好等优点，但在生成含边界层网格时的自动性及可靠性较差，对于复杂外形需要大量手工操作，有时甚至无法生成.高雷诺数黏性流动中不能自动生成可靠且有效的网格仍然是计算流体力学中的瓶颈问题.

边界层是流体运动中贴近壁面的薄层，其厚度与物体特征尺寸的比值为一个小量，通常小于千分之一.如果采用各向同性网格，则生成的边界层网格单元数量巨大，大大提高了计算的代价，对于大规模的复杂问题，其计算代价是目前工程上不能接受的.而如果计算域都采用大长宽比的各向异性单元则会得到质量很差的数值解[1].因此，边界层网格生成的难点在于如何用尽可能少的网格单元达到捕捉层内流动特性的目的.

三维黏性混合网格策略，即在边界层内生成高质量的半结构化三棱柱单元而在其他远场区域生成各向同性的四面体网格单元，是目前最常用的方法，这种方法兼具结构网格和非结构网格的优点.Pirzadeh[2]提出的层推进法(advancing layer method)是其中最具影响的一种.层推进方法从需要生长边界层的表面三角网格出发，沿通过表面顶点的某一个方向布置各向异性网格的节点，然后根据表面网格的拓扑关系连接这些布置的节点来生成一层一层的三棱柱单元，有时由于求解器的原因，三棱柱单元可以转化为四面体，如图1所示.

图1 三棱柱单元生长示意图Fig.1 Prism elements growing example

围绕层推进法，很多学者做了不同程度的改进，Kallingderis等[3]改进了推进向量的求法，提出了可见区域的概念，利用该方法求得的向量可以保证生成的单元都有正的体积.L¨ohner[4]提出了连接模型边界上生成的半结构化网格与其他区域生成的各向同性网格的方法，在半结构化边界层网格生成过程中，检测形状差,反转、相交、尺寸不适合的单元，将相应的单元删除.Aubry等[5]提出了一种求节点法向的方法，并将其用于表面各向异性网格的生成[67]，最近又利用Voronoi图来确定非流形角点处的生长法向，进而更好地生成边界层网格，并集成到软件中[8].Connell等[9]实现了全四面体黏性网格的生成，在确保相容性的前提下将三棱柱单元全部转化为四面体单元，通过删除单元的方式处理相交问题.Garimella等[10]引进多生长方向，提高了非流形角点处的单元正交性及单元质量，提出采用收缩、平滑、裁剪等操作处理层推进过程中的单元自交及前沿相交情况.Athanasiadis等[1112]则提出了构造虚拟点来实现多法向生长及法向融合操作，提高了凹凸区域黏性网格的正交性及结构化特性.Sharov等[13]通过生成各向异性的表面网格来适应表面凹凸区域的半结构化黏性网格生成.王刚等[1415]在这方面做了一些工作，将层推进方法应用到黏性网格生成.陈建军等[16]也在已有的研究基础上，集成各种处理方法，发展了一套混合网格生成程序.

以上所述的方法，在出现相交情况的时候采用删除单元或停止推进的方式进行处理，这样会导致各向同性网格生成器的初始表面网格尺寸不均匀，存在长宽比很大的狭长单元，不利于各向同性网格的生成;同时边界层网格的层状结构也会破坏.另外，该方法在凹凸脊及角点等几何特征处需要特殊处理，现有的处理方法都会产生一定的问题，尚未合理地解决通用性的问题.

除了层推进方法，一些学者也提出了其他方法.Dyedov等[17]将偏移面(face o ff setting)方法[18]应用到黏性混合网格生成中，考虑面追踪，并直接推进表面三角形单元，事先确定推进距离的上限，避免相交检测，并提出采用变分优化的方法对三棱柱网格进行质量改善，该方法通过求解拉格朗日演化方程的方式来生成三棱柱边界层网格，并将该混合网格生成方法应用到生物力学计算中，取得了一定效果.

与层推进方法不同，Ito等[19]以及Alauzet等[20]发展了一种黏性混合网格的生成方法，先生成全域的非结构各向同性网格，然后采用挤压内部网格的方式生长边界层半结构化网格.前者不改变各向同性网格的拓扑结构，后者结合拓扑可变网格变形方法.该方法可以避免网格生长过程中的相交检测，但受限于网格变形技术，由于网格变形能力有限，在凹角点区域可以出现单元反转的现象，同时由于挤压，内部非结构网格的质量下降，影响计算的精度.

张来平等[2123]在生成壁面三棱柱网格方面发展了一种聚合的方法，将物面附近各向异性的四面体单元聚合生成三棱柱单元，从而得到半结构化的边界层网格，并将生成的混合网格应用到非定常数值计算中[24-25].

本文在前人研究成果的基础上，基于前沿层推进方法，发展了一套黏性边界层网格自动生成程序，其中的细节处理方法，包括前沿多法向的确定以及生长高度的调节等，很好地解决了凹凸区域及角点附近的边界层网格生成稳定性问题，能够自动高效地生成可靠且有效的边界层网格.同时发展了一种三维尺寸调节方式，可以为各向同性网格生成器提供均匀尺寸的初始表面网格，进而得到高质量的远场各向同性网格，有利于提高流体计算解的质量.

1 边界层网格生成流程

本文算法步骤如下:

输入：计算区域物面三角形网格，边界层首层厚度hf，厚度增长比例hr，层数nl.

输出：边界层网格，及最外层三角形网格前沿.

Step1读入物面三角形网格数据，建立网格数据结构及拓扑关系.

Step2调整各节点总生长高度，此时生长高度为用户给定总生长高度

Step 2.1提取前沿特征信息，包括特征边及角点等信息，确定节点生长线数量及方向.

Step 2.2根据特征信息生长对应的单元，采用收缩高度的方式处理单元自交(反转)问题，得到各点的合法生长高度.

Step 2.3对新前沿进行全局相交检测，收缩高度避免全局相交，得到各点新的生长高度.采用拉普拉斯方法进行高度光滑处理得到最终的总生长高度.

Step3根据调整后的总生长厚度，得到每层生长高度，依次重复Step2.1—Step2.2，直到达到边界层数nl，对最外层前沿进行全局相交检测，采用收缩高度的方式处理.

2 前沿节点生长法向确定

2.1 定义

边界层网格生成过程中的前沿为表面三角形网格.物面初始三角形网格剖分为初始前沿.前沿上节点 V，其周围相邻节点记为 Pi(i=1,2,···)，对应的邻边记为 ei(i=1,2,···).

V的邻边ei两侧的三角面片之间的有向二面角记为 θi，θi较大 (大于 240°)的边定义为特征边. 邻边中存在三条及以上特征边的节点称为角点.如图2标出了边e3及对应的二面角θ3.

2.2 生长法向计算

计算每个节点周围的有向二面角，确定其邻边中的特征边，根据特征边将节点周围的三角面片分组，如图 2，V点周围存在三条特征边 VP1,VP3,VP4，则将周围的三角面片分为3组.每一组确定一条生长方向ni(i=1,2,3)，由该组三角面片的外法向加权得到，Ni为第i组三角面片数，mj为第 j个三角片外法向，权值wj可简单取为1或者采用文献[26]中提出的角度权值

图2 生长法向确定Fig.2 Calculate growth directions

以上的法向计算策略能够很好地满足推进单元正交性及可视化条件，而且三条法向就足够了.因此本文设定节点的生长方向最多有三条.当节点周围存在三条以上特征边时，则从中挑出有向二面角最大的三条特征边作为最后的分组标准.

节点生长法向确定后，周围三角片在该节点处有一个对应的生成法向；特征边在该点处有一到两个生长方向，可由边两侧的三角片的对应情况确定.当节点存在三个生长法向时，该节点会单独生长出一个四面体单元.

3 生长高度确定

初始生长高度为用户给定，即程序的输入中包含的第一层网格高度hf，生长高度增长比率hr及层数nl.由于节点的生长法向不一定都垂直于物面，因此有必要对初始生长高度进行适当的调整.按照调整后的高度信息逐层推进生长边界层，在推进过程中，由于边界形状的复杂性，可能出现单元自交及前沿全局相交的问题.因此也需要对当前的生长高度进行调节.

3.1 初始生长高度

在用户给定初始生长高度等信息作为输入后，由于上一节确定的节点生长方向不一定垂直于节点周围所有的三角面片，因此为了尽量保证推进前沿与当前前沿的距离符合用户给定的信息，需要对初始生长高度进行适当的调整.

对于二维情况，由于生长法向是节点外角的N等分线 (N为生长法向数目)，因此只需要考虑法向与物面垂线夹角即可.图3为二维示意图，前沿推进高度为h，节点V1有一个法向，则生长高度为h/sinθ，点V2处θ=90°，所以生长高度为h.

图3 二维初始高度调节Fig.3 2D initial height setting

对于三维情况，由于采用多生长方向，同时法向的确定为分组三角面片法向的加权平均，因此为了保证推进前沿与当前前沿尽量平行，则需要采用加权法调整每个节点的初始生长高度.如图4，节点V周围的一个三角片Ti，其单元法向为nT，它对应的分组节点生长法向为ni(节点V可能有多个法向，图4只显示了对应于Ti的法向ni)，则初始高度调整后为

遍历节点V周围所有三角片加权平均得到节点V调整后的初始节点生长高度hV，记T(V)为节点V周围的三角片集合，NT(V)为集合元素个数，则有

图4 三维初始高度调节Fig.4 3D initial height setting

3.2 自交

自交即单元出现反转，如图5(a)所示.当出现自交时，收缩自交单元节点的生长高度，二维单元则收缩至交点处，即法向融合，这样可以避免前沿小尺寸出现，同时也可以减少之后处理自交的频次.

图5 自交处理二维示意图Fig.5 Dealing with 2D self-intersection

对于三维情况，当三棱柱单元出现反转时，如图6，记三角形ABC推进到三角形A′B′C′，面积由S变为S′，则采取收缩高度的方式处理，收缩因子γ为

自交处理完以后，需要采用拉普拉斯方法对前沿各点的生长高度进行光滑处理，这样可以避免三棱柱单元邻边突变的情形.

图6 自交处理三维示意图Fig.6 Dealing with 3D self-intersection

3.3 全局相交

对于复杂的计算域，前沿层推进的过程中，不同物面推进的前沿之间可能存在相交的情况，对于前沿规模非常大的问题，全局相交检测是一个耗时的过程.为了减少相交检测的范围，可以建立八叉树(二维四叉树)结构使得全局检测局部化.尽管八叉树结构提高了全局相交检测的效率，但以往的层推进方法每一层都需要检测，使效率下降.为了提高边界层网格生成的效率，本文的方法先对全局预测一个总生长高度，即先采用总生长高度生长一层边界层，在单元生长的过程中进行自相交处理，得到每个表面节点的生长高度，然后进行拉普拉斯光滑处理，最后进行一次全局相交检测，采用收缩的方式得到每个节点的总生长高度，进而得到合适的每层生长高度.采用预测得到的生长高度进行边界层网格推进则可以减少全局相交检测的次数，提高效率.

图7所示的是一个二维全局收缩的示意图，收缩后的总生长高度能够较好地避免内部边界层网格前沿的全局相交，因此，在很大程度上能够避免多次的全局相交检测及收缩操作.

图7 全局相交处理二维示意图Fig.7 Dealing with global intersection(2D example)

4 单元生长及前沿尺寸调节

4.1 单元生长及细分

根据第2节可知，前沿存在三种特征信息:角点，特征边以及三角片.因此，本文方法的单元生长方式共有如图8四种情况，粗点及粗线在当前前沿上，箭头表示生长方向:角点生长出一个四面体单元，三角片生长出一个三棱柱单元，特征边则可以生长出一个三棱柱单元或者一个退化的三棱柱单元.

图8 单元生长示意图Fig.8 Element growing types

由于大多数流体求解器局限于纯四面体网格[7]，为了适应大多数流体求解程序的要求，本文实现了将得到的边界层网格转化全四面体网格的过程:连接三棱柱或退化三棱柱的侧面对角线，将三棱柱细分为三个各向异性的四面体单元，如图9.

为了保持单元之间的相容性，即单元相邻当且仅当单元具有共同的节点、边或面.然而，不是所有的连接方式都能对三棱柱单元进行有效的剖分，如图9(a)的S型连接方式则不能对三棱柱单元进行剖分得到三个四面体单元.因此在连接三棱柱侧边四边形对角线的时候需要进行特殊的处理.

排除以上的S型连接方式，可以得到合法连接的充要条件是存在两条有公共端点的对角线.如图9(b)连接方式将三棱柱剖分为ABCF,ABFE,AEFD三个四面体.

图9 三棱柱侧面对角线连接方式Fig.9 Connecting types of diagonals for lateral faces of a prism

因此本文处理连接对角线的方式为：对角线统一由编号小的点连向编号大的生长点，如图10编号i1＞i0＞i2，则对角线由Pi0和Pi2连向编号大的点Pi1的生长点，Pi2连向Pi0的生长点.对于全部的三棱柱和退化的三棱柱的侧面都进行如此连接处理，则可以得到全局相容的四面体网格.

图10 节点拓扑连接关系示意图Fig.10 Diagonal choice for prism tetrahedronization

4.2 前沿尺寸调节

边界层网格生成结束之后，边界层最外层前沿将作为各向同性网格生成器的表面网格，因此最外层前沿三角网格的质量直接关系到生成的各向同性网格的质量.所以，很有必要在生成边界层网格的同时进行前沿尺寸调节，为之后的各向同性网格生成提供高质量的三角形网格输入.

在边界层生长推进的过程中，其前沿的尺寸可能不断发生变化.二维情况下，在生成边界层网格时，表面线段在推进后的尺寸可能会变大(变小的情况按自交情况处理)，如图11中的BC推进到AD，AD的长度变大，记

文献[27]在生成二维混合网格时对这种情况进行了处理，提出当α大于给定值(本文取1.5)时，则取AD中点E，将四边形单元ABCD细分为三个三角形.

本文将其推广应用到三维尺寸控制中.由于三维边界前沿在推进过程中只出现四种情形，角点产生四面体，特征边和三角面片产生三棱柱，不管哪种情形，其表面只存在两种形状:三角形或四边形.首先标记尺寸大的边，三角形的处理方式为标记边的中点与对应的顶点相连，然后消除标记，依次进行，直至所有标记边处理结束;而对于四边形，则采用图11的细分方式，同时也加入三角形的操作方式;参照图12，红色粗线为标记的尺寸较大的边，蓝色线为新加入线.对应于三维中的调节操作可参见附录A,另外有两种情况需要进行一次边置换操作，如图13.

图12 三角形及四边形处理方式Fig.12 Subdivision of triangles and quads

图13 边置换操作(蓝色线)Fig.13 Edge swap operation(blue line)

5 边界层网格生成实例及数值验证

程序在输入物面三角形网格之后，能够根据给定的边界层网格第一层高度、层数及总高度，自动生成边界层网格.以下实例的远场网格均采用前沿推进法[28]生成.同时，为了验证本文生成的混合网格的计算有效性，采用本文方法生成的混合网格对二维NACA0012翼型、三维ONERA M6翼型、带发动机短舱F6翼身组合体进行黏性绕流计算.

图14是自交处理和尺寸调节结合的二维实例图.从图中可以看出，调整生长高度融合生长法向之后，单元自交得到很好的处理(凹角区域).

图14 自交处理和尺寸调节二维实例Fig.14 2D example dealing with self-intersection and size controlling

图15和图16为一个凹槽黏性网格图，展示了凹角区域收缩生长高度的效果.

图15 凹槽黏性网格Fig.15 Viscous mesh example of a cavity

图16 凹槽黏性网格内部放大图Fig.16 Enlarged inside view of the cavity mesh

5.1 NACA0012翼型

图17是一个二维NACA0012翼型混合网格生成实例，生成20层边界层网格，含8394个四边形，102个三角形;远场包括11302个三角形.图18为翼型尾部放大图.从二维NACA0012实例可以看出，采用尺寸调节后，最外层前沿线段长度比较均匀，因此边界层网格与远场网格衔接较好.

图17 NACA0012二维混合网格Fig.17 2D hybrid mesh for NACA0012

图18 NACA0012二维混合网格尾部放大图Fig.18 Enlarged view of the 2D hybrid mesh of NACA0012

图19是一个三维NACA0012混合网格生成实例，生成10层边界层网格，包括36819个四面体；远场包含22787个四面体.图20为内部视图，可以看出多生长方向提高了边界层网格的正交性.

图19 NACA0012三维网格Fig.19 3D hybrid mesh of NACA0012

图20 NACA0012三维网格内部视图Fig.20 Inside view of the 3D hybrid mesh of NACA0012

二维NACA0012翼型绕流的入口边界马赫数为0.8，攻角为0°，采用Spalart-Allmaras湍流模型，得到的压力分布云图如图21所示，其上表面壁面附近的速度矢量图如图22所示，从结果可以看出本文采用的混合网格能够捕捉壁面附近及流场内部的信息，网格是有效的.

图21 二维NACA0012压力分布云图Fig.21 Pressure contours of NACA0012 wing

图22 机翼上表面附近速度矢量图Fig.22 Velocity vectors near upper wing surface

5.2 三维ONERA M6翼型

采用本文方法生成M6的黏性混合网格如图23所示，图24为内部视图.

图24 ONERA M6混合网格内部视图Fig.24 Inside view of hybrid mesh of M6 wing

为了与现有实验数据进行对比，三维M6翼型绕流采用标准实验工况，Ma=0.8395，攻角为3.06°，采用k-ω湍流模型，机翼不同站位上的压力系数分布计算值与实验数据[29]对比如图25～图29.选取的站位y/b分别为0.2，0.44，0.65，0.8，0.9.

图25 y/b=0.2压力系数分布Fig.25 Pressure coeffients on the wing surface at y/b=0.2

图26 y/b=0.44压力系数分布Fig.26 Pressure coeffients on the wing surface at y/b=0.44

图27 y/b=0.65压力系数分布Fig.27 Pressure coeffients on the wing surface at y/b=0.65

图28 y/b=0.8压力系数分布Fig.28 Pressure coeffients on the wing surface at y/b=0.8

图29 y/b=0.9压力系数分布Fig.29 Pressure coeffients on the wing surface at y/b=0.9

从图中对比可以看出，压力系数分布总体上吻合得较好，各站位上的激波位置捕捉较为理想.总之，从对比结果可以看出本文生成的黏性混合网格能够有效应用到黏性绕流计算中.

5.3 带发动机短舱的F6翼身组合体

采用本文方法对F6翼身组合体进行黏性混合网格生成，图30为内部网格剖分图.

图30 DLR-F6混和网格内部视图Fig.30 Inside view of hybrid mesh of DLR-F6

计算中，为了与实验数据进行对比，采用Ma=0.75，Re=3×106，α=1.738的计算条件，文献[30]表明计算结果与湍流模型的选取有关，本文的湍流模型采用k-ω模型.将得到的机翼不同站位上的压力系数与实验值[31]进行对比.图31～图38展示了计算结果与所有8个实验站位上的对比图.从图中可以看出，计算结果与实验得到的物面压力系数分布较为符合，同时激波位置的捕捉也较为准确，因此本文生成的混合网格运用到F6翼身组合体的绕流计算是有效的.

图31 y/b=0.15压力系数分布Fig.31 Pressure coeffients on the wing surface at y/b=0.15

图32 y/b=0.239压力系数分布Fig.32 Pressure coeffients on the wing surface at y/b=0.239

图33 y/b=0.331压力系数分布Fig.33 Pressure coeffients on the wing surface at y/b=0.331

图34 y/b=0.377压力系数分布Fig.34 Pressure coeffients on the wing surface at y/b=0.377

图35 y/b=0.411压力系数分布Fig.35 Pressure coeffients on the wing surface at y/b=0.411

图36 y/b=0.514压力系数分布Fig.36 Pressure coeffients on the wing surface at y/b=0.514

图37 y/b=0.638压力系数分布Fig.37 Pressure coeffients on the wing surface at y/b=0.638

图38 y/b=0.847压力系数分布Fig.38 Pressure coeffients on the wing surface at y/b=0.847

6 结论

本文提出的黏性边界层网格自动生成策略可以很好地解决物面凹凸区域及角点处边界层网格的生长问题，避免产生非法单元，同时减少全局相交检测的次数，提高了网格生成的效率;另外还发展了一种三维尺寸调节方法，保证了边界层网格最外层三角形网格的尺寸统一，有利于生成高质量的各向同性远场网格.最后将本文生成的二维及三维的黏性混合网格应用到高雷诺数黏性绕流计算中，得到的计算结果表明本方法生成的黏性混合网格能够满足实际计算的要求，得到合理的计算结果.结合网格变形方法，如弹簧类方法，代数插值类方法[32]等,本文提出的混合网格生成策略也将很快应用到非定常流场计算[3334]中.

本文在处理自交或全局相交时只采用收缩因子的形式，在容易产生自交的凹角区域，其法向推进距离小于其他区域，前沿尺寸也相应减小，如图39所示为凹角区域收缩效果，在该区域的未推进空间过小，生成各向同性网格会有一定困难.目前的处理方式是增加该区域表面网格的密度，减小表面网格尺寸，即生成的初始表面网格在凹凸脊及角点区域的单元会加密.同样，在一些狭缝区域，两侧壁面同时推进的情况下，由于空间狭小，会造成未推进区域过小，给各向同性网格生成造成困难.本文实例中尚未涉及此种特别狭窄的区域，全局的前沿合并策略是解决这类情况的一种好的方式，将在今后的工作中进一步考虑和完善.

图39 凹角网格收缩效果Fig.39 Mesh shrinking at concave ridges

另外，本文的程序将表面网格直接作为输入，计算域边界的表面三角剖分对于边界层网格的顺利进行具有重要意义[6,35]，因此，在生成边界层网格之前，如何获得适用于边界层推进的高质量的表面网格值得进一步研究.

1 Sahni O,Jansen KE,Shephard MS,et al.Adaptive boundary layer meshing for viscous fl ow simulations.Engineering with Computers,2008,24(3):267-285

2 Pirzadeh S.Unstructured viscous grid generation by the advancinglayers method.AIAA Journal,1994,32(8):1735-1737

3 Kallinderis Y,Ward S.Prismatic grid generation for threedimensional complex geometries.AIAA Journal,1993,31(10):1850-1856

4 Lohner R.Matching semi-structured and unstructured grids for Navier–Stokes calculations.AIAA Paper 1993-3348,Orlando,FL,U.S.A.1993

5 Aubry R,L¨ohner R.On the‘most normal’normal.Communications in Numerical Methods in Engineering,2008,24(12):1641-1652

6 Aubry R,Dey S,Mestreau E,et al.An entropy satisfying boundary layer surface mesh generation.SIAM Journal on Scienti fi c Computing,2015,37(4):A1957-A1974

7 Aubry R,Dey S,Karamete K,et al.Smooth anisotropic sources with application to three-dimensional surface mesh generation.Engineering with Computers,2016,32(2):313-330

8 Dey S,Aubry R,Karamete B,et al.Capstone:A geometry-centric platform to enable physics-based simulation and design of systems.Computing in Science&Engineering,2016,18(1):32-39

9 Connell SD,Braaten ME.Semistructured mesh generation for threedimensional Navier-Stokes calculations.AIAA Journal,1995,33(6):1017-1024

10 Garimella RV,Shephard MS.Boundary layer mesh generation for viscous fl ow simulations. International Journal for Numerical Methods in Engineering,2000,49(1):193-218

11 Athanasiadis AN,Deconinck H.Object-oriented three-dimensional hybrid grid generation.International Journal for Numerical Methods in Engineering,2003,58(2):301-318

12 Athanasiadis AN,Deconinck H.A folding/unfolding algorithm for the construction of semi-structured layers in hybrid grid generation.Computer Methods in Applied Mechanics and Engineering,2005,194(48):5051-5067

13 Sharov D,Luo H,Baum JD,et al.Unstructured Navier-Stokes grid generation at corners and ridges.International Journal for Numerical Methods in Fluids,2003,43(6-7):717-728

14 王刚,叶正寅,陈迎春.三维非结构黏性网格生成方法.计算物理,2001,18(5):402-406(Wang Gang,Ye Zhengyin,Chen Yingchun.Generation of three-dimensional unstructured viscous grids.Chinese Journal of Computational Physics,2001,18(5):402-406(in Chinese))

15 王刚,叶正寅.三维非结构混合网格生成与N-S方程求解.航空学报,2003,24(5):385-390(Wang Gang,Ye Zhengyin.Generation of three dimensional mixed and unstructured grids and its application in solving Navier-Stokes equations.Acta Aeronautica et Astronautica Sinica,2003,24(5):385-390(in Chinese))

16 陈建军,曹建,徐彦等.适应复杂外形的三维黏性混合网格生成算法.计算力学学报,2014,31(3):363-370(Chen Jianjun,Cao Jian,Xu Yan,et al.Hybrid mesh generation algorithm for viscous computations of complex aerodynamics con fi gurations.Chinese Journal of Computational Mechanics,2014,31(3):363-370(in Chinese))

17 Dyedov V,Einstein DR,Jiao X,et al.Variational generation of prismatic boundary-layer meshes for biomedical computing.International Journal for Numerical Methods in Engineering,2009,79(8):907-945

18 Jiao X.Face o ff setting:A uni fi ed approach for explicit moving interfaces.Journal of Computational Physics,2007,220(2):612-625

19 Ito Y,Nakahashi K.Improvements in the reliability and quality of unstructured hybrid mesh generation.International Journal for Numerical Methods in Fluids,2004,45(1):79-108

20 Alauzet F,Marcum D.A closed advancing-layer method with connectivity optimization-based mesh movement for viscous mesh generation.Engineering with Computers,2015,31(3):545-560

21 Zhang L,Zhao Z,Chang X,et al.A 3D hybrid grid generation technique and a multigrid/parallel algorithm based on anisotropic agglomeration approach.Chinese Journal of Aeronautics,2013,26(1):47-62

22 Zhong Z,Zhang L P,Xin H E.Hybrid grid generation technique for complex geometries based on agglomeration of anisotropic tetrahedrons.Acta Aerodynamica Sinica,2013,31(1):34-39

23 张来平,张涵信.复杂无粘流场数值模拟的矩形/三角形混合网格技术.力学学报,1998,30(1):104-108(Zhang Laiping,Zhang Hanxin.A Cartesian/triangular hybrid grid solver for complex inviscid fl ow fi elds.Acta Mechanica Sinica,1998,30(1):104-108(in Chinese))

24 张来平,王志坚,张涵信.动态混合网格生成及隐式非定常计算方法.力学学报,2004,36(6):664-672(Zhang Laiping,Wang Zhijian,Zhang Hanxin.Hybrid moving grid generation and implicit algorithm for unsteady fl ows.Acta Mechanica Sinica,2004,36(6):664-672(in Chinese))

25 张来平,常兴华,段旭鹏,等.基于动态混合网格的不可压非定常流计算方法.力学学报,2007,39(5):577-586(Zhang Laiping,Chang Xinghua,Duan Xupeng,et al.Implicit algorithm based on dynamic hybrid mesh for incompressible unsteady fl ows.Acta Mechanica Sinica,2007,39(5):577-586(in Chinese))

26 Th¨urrner G,W¨uthrich C A.Computing vertex normals from polygonal facets.Journal of Graphics Tools,1998,3(1):43-46

27 Dussin D,Fossati M,Guardone A,et al.Hybrid grid generation for two-dimensional high-Reynolds fl ows.Computers&Fluids,2009,38(10):1863-1875

28 Geuzaine C,Remacle JF.Gmsh:A 3-D fi nite element mesh generator with built-in pre-and post-processing facilities.International JournalforNumericalMethodsinEngineering,2009,79(11):1309-1331

29 Schmitt V,Charpin F.Pressure distributions on the ONERA-M6-Wing at transonic Mach numbers,experimental data base for computer program assessment.Report of the Fluid Dynamics Panel Working Group 04,AGARD AR-138,1979.

30 石磊,杨云军,周伟江等.两种湍流模型在高速旋转翼身组合弹箭中的对比研究.力学学报,2017,49(1):84-92(Shi Lei,Yang Yunjun,Zhou Weijiang.A comparative study of two turbulence models for Magnus e ff ect in spinning projectile.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1):84-92(in Chinese))

31 La fl in KR,Brodersen O.Summary of data from the Second AIAA CFD Drag Prediction Workshop.AIAA2004-0555,2004

32 刘中玉,张明锋,聂雪媛等.一种基于径向基函数的两步法网格变形策略.力学学报,2015,47(3):534-538(Liu Zhongyu,Zhang Mingfeng,Wei Xueyuan,et al.A two-step mesh deformation strategy based on radial basis function.Chinese Journal of Theoretical and Applied Mechanics,2015,47(3):534-538(in Chinese))

33 常兴华,马戎,张来平等.基于计算流体力学的“虚拟飞行”技术及初步应用.力学学报,2015,47(4):596-604(Chang Xinghua,Ma Rong,Zhang Laiping,et al.Study on CFD-based numerical virtual fl ight technology and preliminary application.Chinese Journal of Theoretical and Applied Mechanics,2015,47(4):596-604(in Chinese))

34 童福林,李新亮,唐志共.激波与转捩边界层干扰非定常特性数值分析.力学学报,2017,49(1):93-104(Tong Fulin,Li Xinliang,Tang Zhigong.Numerical analysis of unsteady motion in shock wave/transitional boundary layer interaction.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1):93-104(in Chinese))

35 Aubry R.Ensuring a smooth transition from semi-structured surface boundary layer mesh to fully unstructured anisotropic surface mesh//53rd AIAA Aerospace Sciences Meeting and Exhibit,AIAA,Reston,VA,2015,AIAA-2015-1507

AUTOMATIC VISCOUS BOUNDARY LAYER MESH GENERATION1)

Gan Yangke Liu Jianfei2)
(Department of Mechanics and Engineering Science，College of Engineering，Peking University，Beijing 100871,China)

High Reynolds number fl uid fl ows have boundary layers at the wall.To automatically generate robust and valid boundary layer mesh for the simulations is still the bottleneck problem of computational fl uid dynamics.Prisms/Tetrahedra hybrid mesh leads to signi fi cant savings in mesh size and solution costs.However,it’s still difficult to generate prismatic elements of high quality within boundary layers of complex models.Previous advancing layers techniques sometimes lead to invalid meshes and poor quality elements at concave/convex ridges and sharp corners.To improve these situations,we present a strategy for automatically generating viscous boundary layer mesh.In this method,multiple growth directions at ridges and corners are well de fi ned by the dihedral angles around the vertices and the growth heights are adjusted appropriately.Therefore,boundary layer mesh grows well at sharp corners,convex and concave ridges of the domain.We also decrease the number of global intersection checks by prede fi ning the total growth heights before generating elements through one global check,which improves the efficiency of mesh generation.At the same time,we develop a 3D strategy of mesh size control to get a size uniform triangular mesh of the outer boundary of the boundary layer mesh,which is bene fi cial to generate far- fi eld isotropic mesh of high quality.Finally,mesh examples and the viscous fl ow simulations including 2D and 3D are presented.In 3D,the hybrid mesh over the standard ONERA M6 and DLR-F6 con fi gurations are generated with the present method.The numerical results agree very well with experiment data which indicates that the hybrid viscous meshes generated by the proposed method are e ff ective and efficient.

boundary layer mesh，growth normal，feature edge，intersection，size control

O242.21

A

10.6052/0459-1879-17-154

2017–05–04收稿，2017–07–05 录用,2017–07–07 网络版发表.

1)国家自然科学基金(11172005)和国家重点基础研究发展计划(2010CB832701)资助项目.

2)刘剑飞，副教授，主要研究方向：网格生成，计算几何，计算机图形学.E-mail：liujianfei@pku.edu.cn

甘洋科,刘剑飞.黏性边界层网格自动生成.力学学报,2017,49(5):1029-1041

Gan Yangke,Liu Jianfei.Automatic viscous boundary layer mesh generation.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1029-1041

附录A三维边界层网格生成中单元尺寸调节操作

角点生成的四面体Tetraheron generated from one corner

三角片生成的三棱柱Prism generated from one triangle

特征边生成的三棱柱Prism generated from one feature edge