面向降落伞稳态CFD计算的网格生成方法研究

靳宏宇 吴壮志 王奇 贾贺 荣伟,3

面向降落伞稳态CFD计算的网格生成方法研究

靳宏宇1吴壮志1王奇2贾贺2荣伟2,3

（1 北京航空航天大学计算机学院，北京 100191）（2 北京空间机电研究所，北京 100094）（3 北京市航空智能遥感装备工程技术研究中心，北京 100094）

降落伞稳态流体动力学（Computational Fluid Dynamics，CFD）计算的网格生成是一个带边界层的三维复杂域网格生成问题，目前用于边界层计算的混合网格生成方法，往往存在计算繁琐、应用范围窄、自动化差、以及难以适应复杂外形的缺点。文章提出了一种结合约束德洛内（Delaunay）网格生成和网格前沿推进技术的降落伞稳态CFD计算的流场网格生成方法，实现了降落伞网格和包括边界层区域在内的流场网格一体化、全自动、高品质的生成。该方法采用网格前沿推进法来生成边界层区域的网格节点；算法整体上采用三维约束Delaunay网格生成技术，边界层层节点集合在网格生成过程中作为约束Delaunay三角化的约束条件，避免了复杂的网格求交计算和拓扑处理。网格生成实例表明，该方法能全自动生成降落伞稳态CFD网格，生成的网格品质、网格规模等满足降落伞稳态CFD计算的需求，同时该方法具有一定的通用性。

约束德洛内三角化 网格生成 边界层 降落伞 航天返回

0 引言

网格生成是计算流体动力学（Computational Fluid Dynamics，CFD）数值计算的基础和前提，占据整个计算周期人力时间的60%左右，而且网格品质的好坏直接关系到计算结果的精度[1-4]。依照网格的拓扑结构，网格可分为结构网格、非结构网格和混合网格[5-7]：1）结构网格主要是指二维的四边形网格和三维的六面体网格，其优点是数据结构简单、计算简单快捷，其缺点是难以适应复杂外形；2）非结构网格主要是指二维的三角形网格和三维的四面体网格，其优点是易于处理复杂外形，并且网格易于自适应。算法主要有四/八叉树法、德洛内（Delaunay）三角化法和前沿推进法（Advancing-Front Method，AFM）[5]，其中Delaunay三角化算法以其完备的数学理论依据，成为目前应用和研究最为广泛的全自动三角形网格生成方法[6,8-9]；3）混合网格是结构化和非结构化网格之间的一个折中[10-15]，用来解决CFD计算中的边界层网格生成问题。它结合了结构化网格的正交性和方向性、非结构化网格的灵活性和易用性，已经成为一种处理复杂几何外形的新型、有效的网格生成技术。

降落伞稳态CFD计算的网格生成是一个带边界层的三维复杂域网格生成问题。边界层是指粘性流体沿固体表面流动或固体在流体中运动时，附于固体表面的一层流体。受到粘性作用的影响，在紧贴物体表面的地方，边界层的流速沿着物体法线方向从零开始逐渐增大，并在很短的距离内增大到跟边界层外流体相同的速度。为了满足边界层流速变化快的特性，在进行粘性流体计算时，需要在边界层生成细小紧密的层次结构化网格，这样不仅可以节省内存，还可以有效地捕捉很强的方向粘性应力。在边界层网格生成问题上，学者们做了大量的工作[16-22]，提出了多种用于边界层计算的混合网格生成方法。但是，对于很多极端复杂的实际外形而言，由于几何曲率变化剧烈，所生成的边界层网格常常出现网格相交的现象，为此需要付出很多人工劳动调整局部网格分布，极大地制约了CFD的应用[11]。而且，由于在边界附近所生成的层状单元（如三棱柱）纵横比较大，单元品质不高，整体上网格为混合单元结构，导致这些方法存在计算繁琐、应用范围窄、自动化差以及难以适应复杂外形的缺点。

针对混合网格生成方法存在的问题，本文提出了一种结合约束Delaunay网格三角化（Constrained Delaunay triangulation，CDT）和网格前沿推进技术的降落伞稳态CFD计算的流场网格生成方法，实现了降落伞网格和包括边界层区域在内的流场网格一体化、全自动、高品质的生成。Delaunay方法由于具有深刻的理论背景，能够有效保证算法的收敛性和网格品质，同时由于采用非结构化单元，易于适应复杂外形；而网格前沿推进法则用来计算边界层区域每个推进层的前沿节点集合，推进高度由一个几何级数来确定，边界层层节点集合在网格生成过程中作为约束Delaunay三角化的约束条件，避免了网格相交现象的发生，同时避免了复杂的网格求交计算和拓扑处理。

1 网格生成相关理论

1.1 分片线性复形

一个三维分片线性复形（Piecewise Linear Complex，PLC）[8]是指由一组统称为单元的顶点、边、多边形和多面体组成的集合，该集合满足如下条件：

（a）三维PLC实例（a）A 3D PLC（b）三维非PLC实例（b）3D non-PLCs

1）中每个单元的边界是中一组单元的并集；

2）两个不同的单元，∈的交集是中一组单元的并集。

一个三维PLC的底空间用||来表示，定义为的所有单元的并集，它是一个其拓扑由给定的拓扑空间。的边界复形是的所有单元集合的一个子集，其维数小于3。||就是本文要剖分的三维复杂域。||的边界是的边界复形的底空间。

1.2 约束Delaunay三角化

三维分片线性复形中的线段和多边形约束了进行三角化的方式。的三角化定义为这样一个复形：1）和有相同的顶点集；2）的每个单元是中一系列单纯形（顶点、线段、三角形、四面体）的并集；3）|| = ||。注意：||和||并不一定是凸的。的网格是∪的三角化，其中为插入中的不同于中已有顶点的新的点集，新增的点也称为施泰纳（Steiner）点。的三角化不允许增加新的点，但的网格生成允许增加新的点。

通常，的网格将中的每个多边形细分为许多三角形，将中的线段细分为许多边。中一条边如果是中一条线段的一部分，则称为一个子线段；中一个三角形如果是中一个多边形的一部分，则称为一个子多边形。

本文中描述的网格生成算法主要就是要构造一个输入PLC的Steiner CDT。

1.3 网格品质

图2 四面体τ的半径边比：最右边为sliver四体面

2 网格生成方法

2.1 问题描述

降落伞稳定下降阶段的数值模拟，可以将降落伞作刚体假设，降落伞的外形由实验数据给出，降落伞的流场不涉及流固耦合问题，仅视为纯粹的流场计算，本文称为降落伞稳态的CFD计算问题。

本文以平面圆形伞为例来进行降落伞稳态CFD计算网格生成研究。如图3所示，平面圆形伞的基本参数如下：伞顶孔直径0=90mm，伞外圆直径（即伞衣结构直径）1=1 275mm，伞衣幅数=12。稳态降落伞周围的流场设置为圆柱形，流场的位置与大小参数设置如下：圆柱流场的中心轴线与稳态降落伞的中心轴线重合；圆柱流场的底面圆直径2=51；圆柱流场的高为1+2，其中1=51为伞衣底边到圆柱顶面的高，2=51为伞衣底边到圆柱底面的高。

（a）伞衣参数（a）Parameters of canopy（b）流场参数（b）Parameter of flow field

降落伞稳态是指降落伞充气展开进入稳定后的降落伞形态。由于降落伞充气稳定后，降落伞伞衣形态变为空间曲面，因此采用三角网格来进行逼近描述。伞衣表面附近的边界层的网格节点采用前沿推进方法生成，推进高度由一个几何级数来确定，边界层节点集合在网格生成过程中作为约束Delaunay三角化的约束条件。流场的边界由两个顶面圆和侧面圆柱组成，可由两个等边的边形来逼近顶面和底面，由个矩形来逼近侧面，用来表示。

降落伞稳态三角网格、边界层节点集合和流场的逼近边界构成一个3维PLC，其底空间||即为本文的待剖分域。

降落伞稳态CFD计算网格生成问题是指：1）读入降落伞稳态三角网格，生成边界层节点集合和流场的逼近边界，并将、、构成一个三维分片线性复形，将其底空间||作为待剖分域；2）采用基于Delaunay的非结构网格生成和细化方法，将||剖分为一组高品质的四面体单元集合；3）将输出作为降落伞稳态CFD计算网格。

2.2 边界层Steiner点集生成

混合网格生成是边界层网格生成的主要方法之一，研究发现，目前的三棱柱/四面体混合网格生成方法中，对凹几何构造或由多几何区域组成的几何体，界面在向前推进时，对于产生的网格节点（单元）冲突或消失的情况没有给出有效的自动解决方法[10]。

针对降落伞稳态CFD网格生成问题，本文采用的非结构网格生成技术与混合网格生成方法的不同点在于：基于Delaunay的非结构网格生成技术，将边界层和外流场统一处理；不需要关注边界层的界面，仅仅关注边界层的顶点生成，边界层顶点生成后作为输入的分片线性复形的顶点的一部分，由网格生成算法统一处理，避免了处理边界复杂的拓扑和自相交等问题。边界层顶点生成技术整体上采用前沿推进法：

1）确定边界层网格推进方向。本文采用前沿推进法来确定边界层网格推进方向，由于伞衣内外都是流场，因此需要从伞衣表面开始向内和向外两个方向进行边界层网格推进。下面以向外推进为例来说明主要推进策略，描述如下：1）以伞衣表面的三角网格的顶点为初始位置进行推进：对两片伞衣连接处的顶点（如图4（b）的、点），同时往两个方向进行推进。推进方向由该顶点的一环邻域在同一片伞衣中的一环邻域三角形的法矢共同决定[24]，例如点的推进方向1由三角形1，2，3法矢的加权平均确定，推进方向2由三角形4，5，6法矢的加权平均确定；2）伞衣内部顶点和不在连接线处的伞衣边界顶点则向单法矢方向推进（如图4（a）的、点）。其推进方向由顶点的一环邻域三角形的法矢的加权平均确定，例如：的推进方向由三角形712的法矢确定，的推进方向由三角形9、10、13的法矢确定。

（a）单方向推进点：伞衣内部和不在伞衣连接处边界顶点（a）Unidirectional propulsion points: the vertices at the interior or disconnected boundary of a canopy（b）两方向推进点：两片伞衣连接处的顶点（b）Two-directional propulsion points: the vertices at the connection of two canopies

2.3 网格生成流程

本文以开源四面体网格生成（Tetrahedral mesh Generation，TetGen）[8]算法为基础来进行CFD网格生成，网格生成流程主要包括如下4个步骤：

1）3维PLC的三角化。此步为预处理，将中的所有边界面进行2维三角化，使得仅仅包含顶点、边、三角形单元；

2）点集的Delaunay四面体化。以的所有顶点集合为输入，进行点集Delaunay四面体化，生成顶点集合的四面体网格DT；

3）约束Delaunay网格生成。在DT作为初始的四面体网格的基础上，将中的所有边和三角形作为约束条件，进行约束Delaunay四面体网格生成。根据给定的规则，约束条件可以整体保留，或者进行细分保留；

4）保品质网格生成。此步骤主要考虑网格单元品质，根据用户指定的网格单元品质要求完成保品质的网格生成。又分为2个子步骤来完成：一是通过在网格中插入新点来提高单元品质；二是通过局部网格优化来提高单元品质。根据指定的规则不同，约束可以保持不被修改，或者可以进行细分保留来进一步提高网格单元品质。

3 网格生成实例与分析

（a）降落伞稳态曲面三角网格（a）Triangular mesh for parachute steady surface（b）边界层Steiner点集（b）Steiner point set of boundary layer

图6给出了本文方法生成的降落伞稳态CFD计算网格。从局部放大图可以看出：网格从伞衣表面开始向外、向内各为5层，并且所生成的四面体符合各向异性的特点；由于各个边界层并不需要维护界面的拓扑结构，因此生成方法全自动并且十分鲁棒。同时，伞衣附近的四面体单元尺寸较小，离降落伞越远，四面体尺寸越大，网格整体上呈现自适应特点。

（a）网格剖面视图1（a）Mesh section view 1（b）网格剖面视图2（b）Mesh section view 2（c）边界层网格局部放大（c）Local enlargement of boundary layer mesh（d）边界层网格局部放大（大倍数）（d）Local enlargement of boundary layer mesh（large multiplier）

4 结束语

本文提出了一种结合约束Delaunay网格三角化和网格前沿推进技术的降落伞稳态CFD计算的流场网格生成方法，实现了降落伞网格和包括边界层区域在内的流场网格一体化、全自动、高品质的生成。通过网格生成实例验证，能够全自动生成降落伞稳态CFD网格，生成的网格品质、网格规模等满足降落伞稳态CFD计算的需求，同时该方法具有通用性，适用于三维复杂域的非结构化网格生成。该方法具有如下优点：

1）整体上采用三维约束Delaunay网格三角化技术，具有深刻的理论背景，能够有效保证算法的收敛性和网格品质，同时由于采用非结构化单元，易于适应复杂外形；

2）采用网格前沿推进法实现了边界层区域每个推进层的前沿节点集合的计算，实现了降落伞网格和包括边界层区域在内的流场网格一体化、全自动的生成。前沿节点集合在网格生成过程中作为约束Delaunay三角化的约束条件，避免了网格相交现象的发生、以及复杂的网格求交计算和拓扑处理。

本文的下一步工作包括：1）进行CFD仿真计算，进一步验证网格品质对仿真计算的影响；2）开展和其他网格生成方法（如混合网格生成方法）的对比研究，并在此基础上对本文提出的方法进行改进提高。

[1] KALLINDERIS Y, WARD S. Prismatic Grid Generation for Three-dimensional Complex Geometries[J]. AIAA Journal, 1993, 31(10): 1850-1856.

[2] PIRZADEH S. Unstructured Viscous Grid Generation by the Advancing-layers Method[J]. AIAA Journal, 1994, 32(8): 1735-1737.

[3] 张来平, 赫新, 常兴华, 等. 复杂外形静动态混合网格生成技术研究新进展[J]. 气体物理, 2016, 1(1): 42-61. ZHANG Laiping, HE Xin, CHANG Xinghua, et al. Recent Progress of Static and Dynamic Hybrid Grid Generation Techniques over Complex Geometries[J]. Physics of Gases, 2016, 1(1): 42-61. (in Chinese)

[4] BAKER T J. Mesh Generation: Art or Science[J]. Progress in Aerospace Sciences, 2005, 41(1): 29-63.

[5] OWEN S J. A Survey of Unstructured Mesh Generation Technology[C]//In Proceedings of the 7th International Meshing Roundtable. Sandia National Laboratories, 1998: 239-267.

[6] CHENG S W, Dey T, SHEWCHUK J. Delaunay Mesh Generation[M]. CRC Press, 2012.

[7] 金颖. 基于Delaunay三角化的混合网格生成技术研究[D]. 北京: 北京航空航天大学, 2014. JIN Ying. The Research of Hybrid Grid Generation Technology Based on Delaunay Triangulation[D]. Beijing: Beihang University, 2014. (in Chinese)

[8] HANG S. TetGen, a Delaunay-based Quality Tetrahedral Mesh Generator[J]. ACM Transactions on Mathematical Software, 2015, 41(2): 1-36.

[9] SHEWCHUK J R, HANG S. Higher-quality Tetrahedral Mesh Generation for Domains with Small Angles by Constrained Delaunay Refinement[C]//Proceedings of the Thirtieth Annual Symposium on Computational Geometry. 2014: 290-299.

[10] SHAROV D, NAKAHASHI K. Hybrid Prismatic/Tetrahedral Grid Generation for Viscous Flow Applications[J]. AIAA Journal, 1998, 36(2): 157-162.

[11] 赵钟, 张来平, 赫新. 基于“各向异性”四面体网格聚合的复杂外形混合网格生成方法[J]. 空气动力学学报, 2013, 31(1): 34-39. ZHAO Zhong, ZHANG Laiping, HE Xin. Hybrid Grid Generation Technique for Complex Geometries Based on Agglomeration of Anisotropic Tetrahedrons[J]. Acta Aerodynamica Sinica, 2013, 31(1): 34-39. (in Chinese)

[12] 陈建军, 曹建, 徐彦, 等. 适应复杂外形的三维粘性混合网格生成算法[J]. 计算力学报, 2014, 31(3): 363-370. CHEN Jianjun, CAO Jian, XU Yan, et al. Hybrid Mesh of Generation Algorithm for Viscous Computations Complex Aerodynamics Configurations[J]. Chinese Journal of Computational Mechanics, 2014, 31(3): 363-370. (in Chinese)

[13] ATHANASIADIS A N, DECONINCK H. Object-oriented Three-dimensional Hybrid Grid Generation[J]. International Journal for Numerical Methods in Engineering, 2003, 58(2): 301-318.

[14] 王刚, 叶正寅, 陈迎春. 三维非结构粘性网格生成方法[J]. 计算物理, 2001, 18(5): 402-406. WANG Gang, YE Zhengyin, CHEN Yingchun. Generation of Three-dimensional Unstructured Viscous Grids[J]. Chinese Journal of Computational Physics, 2001, 18(5): 402-406. (in Chinese)

[15] DUSSIND, FOSSATI M, GUARDONE A, et al. Hybrid Grid Generation for Two-dimensional High-Reynolds Flaws[J]. Computers & Fluids, 2009, 38(10): 1863-1875.

[16] 甘洋科, 刘剑飞. 粘性边界层网格自动生成[J]. 力学学报, 2017, 49(5): 1029-1041. GAN Yangke, LIU Jianfei. Automatic Viscous Boundary Layer Mesh Generation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(5): 1029-1041. (in Chinese)

[17] AULARY R, DEY S, MESTREAU E, et al. An Entropy Satisfying Boundary Layer Surface Mesh Generation[J]. SEAM Journal on Scientific Computing, 2015, 37(4): 1957-1974.

[18] SALNI O, JANSEN K E, SHEPHARD M S, et al. Adaptive Boundary Layer Meshing for Viscous Flow Simulations[J]. Engineering with Computers, 2008, 24(3): 267-285.

[19] ATHANASIADIS A N, DECONINCK H. A Folding/Unfolding Algorithm for the Construction of Semi-structured Layers in Hybrid Grid Generation[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(48): 5051-5067.

[20] SHAROV D, LUO H, BAUM J D, et al. Unstructured Navier-stokes Grid Generation at Corners and Ridges[J]. International Journal for Numerical Methods in Fluids, 2003, 43(6): 717-728.

[21] DYEDOV V, EINSTEIN D R, JIAO X, et al. Variational Generation of Prismatic Boundary-layer Meshes for Biomedical Computing[J]. International Journal for Numerical Methods in Engineering, 2009, 79(8): 907-945.

[22] ALAUZET F, MARCUM D. A Closed Advancing-layer Method with Connectivity Optimization-based Mesh Movement for Viscous Mesh Generation[J]. Engineering with Computers, 2015, 31(3): 545-560.

[23] CHENG SW, DEY T K, EDELSBRUNNER H, et al. Sliver Exudation. [J]. Journal of the ACM, 2000, 47(5): 883-904.

[24] THURRNER G, WUTHRICH C A. Computing Vertex Normals from Polygoval Facets[J]. Journal of Graphics Tools, 1998, 3(1): 43-46.

Research on Mesh Generation Method for Parachute Steady-state CFD Computing

JIN Hongyu1WU Zhuangzhi1WANG Qi2JIA He2RONG Wei2,3

（1 School of Computer Science & Engineering, Beihang University, Beijing 100091, China）（2 Beijing Institue of Space Mechanics & Electricity, Beijing 100094, China）（3Beijing Engineering Research Center of Aerial Intelligent Remote Sensing Equipments, Beijing 100094, China）

The mesh generation of parachute steady-state CFD computation is a three-dimensional complicated domain mesh generation problem with boundary layers. At present, the hybrid mesh generation method for boundary layer computation often has the disadvantages of cumbersome calculation, narrow application scope, poor automation and difficult to adapt to complicated domain shape. This paper presents a flow field grid generation method for parachute steady-state CFD calculation based on constrained Delaunay mesh generation and advancing-front technology, which achieves full automatic and high quality mesh generation of parachute and flow field mesh including boundary layer region. In this method, mesh nodes in the boundary layer region are generated by the advancing-front method, and the three-dimensional constrained Delaunay mesh generation technology is adopted as a whole. The node set in the boundary layer region created by the advancing-front method is used as the constrained condition of the constrained Delaunay triangulation in the process of mesh generation, thus avoiding the complicated mesh intersection calculation and topological processing. Mesh generation examples show that the method can automatically generate parachute steady-state CFD meshes. The quality of the generated mesh and the number of the its elements meet the requirements of parachute steady-state CFD calculation. At the same time, the method has certain generality.

constrained Delaunay triangulation; mesh generation; boundary layer; parachute; spacecraft recovery

TP391, V11

A

1009-8518(2019)04-0030-08

10.3969/j.issn.1009-8518.2019.04.004

靳宏宇，男，1991年生，2018年获北京航空航天大学计算机应用技术专业工程硕士学位，软件开发工程师。研究方向为网格生成、计算机视觉、机器学习。E-mail：jinhongyu@buaa.edu.cn。

吴壮志，男，1969年生，2001年获北京航空航天大学计算机软件与理论专业博士学位，副教授。研究方向为网格生成、计算几何和计算机图形学。E-mail：zzwu@buaa.edu.cn。

2018-01-25

（编辑：庞冰）