降落伞网格剖分技术研究

吴壮志杨倩荣伟方世兴

（1 北京航空航天大学计算机学院，北京 100191）

（2 北京空间机电研究所，北京 100094）

降落伞网格剖分技术研究

吴壮志1杨倩1荣伟2方世兴2

（1 北京航空航天大学计算机学院，北京 100191）

（2 北京空间机电研究所，北京 100094）

降落伞几何建模和网格剖分是降落伞流固耦合数值仿真的基础。文章建立了包括伞衣主体、径向带、伞顶孔和底边加强带的降落伞精细几何模型，并基于软件开发平台VS2010使用C99语言设计并实现了针对该几何模型的网格自动剖分器。该剖分器采用约束Delaunay剖分算法，具有良好的理论基础、生成的网格质量（quality）高，满足降落伞数值仿真对网格质量的要求；并提供接口可将顶点和网格信息导入到商用软件中做进一步求解分析。实验证明，使用文中网格剖分器生成的网格，不但网格单元质量优于使用商用软件生成的网格，而且在网格单元质量相似的情况下，网格单元数量相对较少。同时，该网格在数值仿真中表现良好。

网格剖分 约束德洛内三角化 网格质量 降落伞 航天返回

0 引言

随着计算机仿真技术的不断进步，降落伞流固耦合技术也得到了快速发展。目前降落伞开伞过程的仿真主要借助于通用的有限元软件来完成的[1-5]，如HyperMesh，Ls-dyna，Abaqus，Ansys等，其仿真过程包括三个部分：前处理、数值求解和后处理。曾经有人做过统计，在数值模拟仿真的三个阶段中，前处理约占总时间的40%～60%，数值求解约占5%～20%，后处理约占30%[6]。前处理过程中的几何建模和有限元网格剖分工作比较繁琐，但是又十分重要，是进行数值模拟仿真的基础。

目前在对降落伞进行几何建模和网格剖分时，仍然有几点不足之处：

第一，通常采用 CAD软件对降落伞进行几何建模，虽然降落伞的形状大小都比较类似，但是只要形状或尺寸有一点不同的地方，都需要重复建模过程，非常费时费力。

第二，通用的商用有限元软件生成的网格单元是整齐划一的，而仿真过程中伞顶孔位置属于大变形区域，需要生成较稠密的网格。为了满足这一需求，通常有两种做法：一是将整体网格密度加大，达到大变形区域的网格密度，但是这样会花费成倍的计算时间，浪费资源；二是人工分区域分别进行不同网格密度的剖分，虽然可以达到对网格密度的要求，但是费时费力。

第三，真实的降落伞结构包括伞衣主体、径向带，伞顶孔和底边加强带等，但是通常只针对伞衣主体进行建模和网格剖分，简化径向带、伞顶孔和底边加强带为一维梁单元[3]（Beam）。

针对上述三点不足，本文建立了包括伞衣主体、径向带、伞顶孔和底边加强带的降落伞精细的几何模型，并对降落伞的网格剖分技术进行研究，基于软件开发平台VS2010使用C99语言设计并实现了一个针对该几何模型的网格自动剖分器。该剖分器具有如下特点：采用约束Delaunay剖分算法，具有良好的理论基础、生成的网格质量（quality，下同）高，满足降落伞数值仿真对网格质量的要求；提供接口可以将生成的顶点和网格信息转换成K文件，导入到HyperMesh等通用的商用软件中做进一步分析计算。

1 网格剖分相关理论

由于非结构化网格可以很好适应复杂自由边界，并且可以很容易控制网格尺寸，所以本文将对降落伞生成非结构化网格。典型通用的非结构化网格生成方法有：四叉树法、前沿推进法和Delaunay三角剖分算法。文献[7]曾对用于商业和科研项目的81种网格生成软件产品进行了非正式的调查，其中37种软件使用的是 Delaunay算法。由此可以看出，Delaunay算法占据着重要的位置。所以本文拟采用约束Delaunay算法（Constrained Delaunay Triangulation）对降落伞进行网格剖分。

1.1 Delaunay三角化

首次给出Delaunay三角网格概念的是俄国著名数学家，Boris Nikolaevich Delaunay。这里给出其精确定义：

设t是三角网格T的一个网格单元，若t的外接圆的内部不包含任何顶点，则称t是Delaunay的，也称其满足空圆特性，如图1所示。设P为二维空间R2中的有限点集，T为P的三角网格，若T中任何一个三角形t都是Delaunay的，则称T是P的一个Delaunay三角化，记为DT(P)。

基于以上定义可知，Delaunay三角化最重要的性质就是空圆特性。为了使三角网格单元满足空圆特性，Lawson提出了局部优化过程[8]（Local Optimization Procedure，LOP）。另外，Delaunay三角化还具有一些其他的优良性质，这里就不一一赘述[9-11]。

1.2 约束Delaunay算法

约束Delaunay算法，顾名思义是在Delaunay算法的基础上加了一些约束条件，包括约束点集V和约束边集 E，严格要求约束条件中的每条约束边是三角剖分后某些三角形网格单元的边，则三角网格应满足约束Delaunay属性[12]。

设点v′为任一点，t是一个三角网格单元，e′为约束边集E中的线段，若v′与t内部任何一点的连线都不与e′相交，则称点v′是对于t内部可见的点。如图2所示，w′为顶点，其中点v′对t是可见的，点w′对t是不可见的。

设t是T中的一个网格单元，若t的外接圆中不包含任何对于t可见的点，则称t是约束Delaunay的，也称约束Delaunay属性。

1.3 网格单元质量评价指标

在计算几何中，Miller等[13]指出，分析一个三角网格单元质量最好的指标是单元最小角、外接圆半径与最短边之比r1和纵横比r2。一个三角网格单元t的纵横比r2为t的外接圆半径R与内切圆半径r之比。其中最直观的就是三角网格单元的最小角，如果最小角越小，则该三角形越畸形。

2 降落伞的网格剖分

2.1 降落伞网格剖分特点

本文以某平面圆形伞为例，对该伞进行精细结构建模和网格剖分的研究。其中伞衣结构包括伞衣主体、径向带、伞顶孔和底边加强带，结构如图3所示。伞衣有12幅，径向带、伞顶孔和底边加强带的个数与伞衣个数相同，每幅伞衣都添加了一条边，即折痕，为之后的折叠做准备。具体的结构参数如表1所示，其中径向带宽度与加强带宽度值一样。

表1 降落伞结构参数表Tab.1 Structure parameters of a parachute

降落伞共包括N0幅伞衣，每幅伞衣精细的几何结构如图4所示（第i幅伞衣），主要包括：伞衣主体（四边形 v11v13v14v10），折痕（线段v1v8），伞顶孔加强带（四边形 v10v14v7v+）、底边加强带（四边形 v0v2v13v11）和径向带（四边形 v2v3v6v7）。将N0幅伞衣合成一体得到整个降落伞伞衣的几何模型。

2.2 网格剖分流程

本文拟采用约束 Delaunay算法对降落伞精细几何模型进行有限元网格剖分，算法描述如下：

输入：降落伞的约束点集V和约束边集E，网格尺寸控制指标S，网格质量控制指标B。

输出：降落伞网格剖分DT。

算法：1）根据输入的约束点集合V，生成初始Delaunay剖分DT0；2）检查约束边集E是否在初始网格剖分DT0中，恢复E中丢失的约束边集E′，得到网格剖分DT1；3）对恢复约束边集后的网格DT1调用约束Delaunay优化算法，得到最终的网格剖分DT；4）保存顶点和网格的信息。

2.3 初始Delaunay剖分

本文使用逐点插入的Lawson算法对约束点集V进行初始Delaunay剖分。基本原理为：首先建立一个大的三角形或多边形，把所有约束点包围起来，取V中一点，该点与包含它的三角形三个顶点相连，形成三个新的三角形，然后逐个对它们用LOP进行优化，以保证所生成的三角网格满足Delaunay的空圆属性。

经过上述步骤可得到包括约束点集V的初始Delaunay剖分DT0，如图5所示。但是DT0可能不完全包括约束边集E，对于丢失的约束边集E′，需要对其进行恢复。

2.4 约束边的恢复

通过在丢失的约束边中点处添加顶点，并经过LOP局部优化，恢复丢失的约束边E′。如果还存在丢失的子约束边，则继续迭代地在其中点处添加顶点，如图6所示，得到恢复约束边集后的网格DT1。约束边e可能被分割为k条子约束边 e1, e2, e3,… ,ek，约束边集E中将不包含e，而包含 e1, e2, e3,… ,ek。

2.5 约束Delaunay优化算法

经过约束边恢复的网格一般都很粗糙，而且比较细长，需要经过加密、匀称等网格优化处理才能用于数值分析。网格优化[14]是整个算法的核心，其中关键部分为：

1）如何控制生成的三角网格单元的质量：一般地，通过设置r1的上界B来控制网格单元的质量，本文设置B为1，即除输入条件中的尖角，生成的网格单元内角都会控制在[30°,120°]，很好的控制了三角网格单元的质量；

2）如何控制生成的三角网格单元的尺寸：通常设置三角网格单元边长L的上界Lmax来控制网格单元的尺寸，即所有网格单元的边长都须小于Lmax；

3）如何控制优化算法的终止：约束Delaunay优化算法是通过不断遍历三角网格单元，对网格单元质量和尺寸不符合要求的网格单元进行优化。这是一个迭代优化的过程。

2.5.1 顶点生成方法

不论是控制三角网格单元的质量还是控制三角网格单元的尺寸，都需要通过添加新的顶点来完成。新顶点的生成需要遵循两个原则：

1）以子约束边为直径的圆是包含子约束边的最小圆，若该圆内部包含其他顶点，则在该边中点处添加新的顶点，迭代添加顶点，直到各子约束边的直径圆内部不包含其他顶点；

2）如果r1＞B或边长大于Lmax，则说明该网格单元是畸形的，或该网格单元尺寸过大。那么在该网格单元的外接圆圆心处添加新的顶点，并去除该网格单元。然而，如果新顶点在某子约束边的直径圆内部，则不插入该顶点，并将相应的子约束边在中点处分割。该算法描述为：

输入：网格DT，约束边集合E；

输出：不包含畸形和尺寸过大的网格单元的网格剖分DT和新的约束边集合E；

遵循上述两个原则，一方面，约束插入的顶点符合约束Delaunay准则；另一方面，设置r1的上界B和边长的上界Lmax，很好地控制了网格单元的质量和尺寸。

2.5.2 输入条件中尖角的处理方法

对于尖角的处理方法，根据原则二插入其所在三角形的外接圆圆心来处理。但是若输入的约束条件中存在尖角[15-16]，根据原则二，会不断插入新的顶点，进入死循环，如图 7所示。对于这种情况，可按如下方法来处理。

若该顶点的插入半径rv小于其父顶点的插入半径 rpv时，停止继续插入新的顶点。

其中关于插入半径和父顶点的定义如下。根据插入半径分为四种情况：若顶点v为约束点，则其插入半径rv为v到对v可见的最近约束点的欧氏距离；若v为某约束边的中点，则rv为v到距其最近的网格顶点的距离；若v为被拒绝插入的某三角形外接圆圆心，则rv为该约束边长的一半；若v为某三角形外接圆圆心，则rv为外接圆半径。

输入条件中的约束点没有父顶点，其余每个顶点 v，包括被拒绝插入的顶点，都有父顶点pv。若 v为某约束边的中点，则pv为该约束边任一顶点；若v为插入或被拒绝的某三角形的外接圆圆心，则pv为该三角形最短边的两顶点中最后插入的顶点。

那么，根据上述条件，约束边分割算法为：

输入：三角网格DT，约束边集合E。

输出：子约束边直径圆中不包含顶点的三角网格DT和新的约束边集合E。

算法：1）for e∈Edo；2）while e的直径圆内存在其他顶点将该边e进行分割，形成新的三角网格；更新约束边集合E；4）调用LOP优化算法；5）endwhile；6）endfor。

2.5.3 终止条件

约束Delaunay优化是一个迭代的过程，优化算法的终止条件为所有网格单元t满足：；②边长。

综上所述，约束Delaunay优化算法可表述为：

输入：三角网格DT，约束边集合E。

输出：优化后的约束Delaunay三角网格DT。

算法：1）while不满足条件①②do；2）调用约束边分割算法；3）调用畸形网格单元调整算法；4）endwhile；5）保存所有顶点和网格信息。

2.6 网格输出处理

Keyword文件简称K文件，是商用有限元软件可接收的文件格式之一。K文件就是在一组关键词组织下的数据块组成的。关键词必须包含符号“*”作为标识，“*KEYWORD”和“*END”分别是K文件的开始符和结束符。

常用的关键词包括：*NODE，*ELEMENT，*PART等。*NODE表示顶点信息，NID表示顶点编号，X，Y，Z表示顶点坐标。*ELEMENT表示网格信息，EID是其编号，PID表示*PART编号，N1，N2，N3分别表示网格单元顶点的NID。图8展示了K文件的组织方式以及各种实体之间的关系。

本文根据*NODE和*ELEMENT内容以及编号、坐标等格式组织网格信息，最后生成的文件后缀设为“.k”，即为K文件，可导入HyperMesh等商用有限元软件做进一步仿真分析。

3 网格优劣评价

本文以某平面圆形伞为例，建立其精细几何模型，并生成三角网格。作为对比试验，使用HyperMesh软件对该几何模型生成三角网格。

3.1 网格剖分结果

图9是采用本文方法对整个降落伞生成的网格，其中黑色部分为伞衣，蓝色部分为加强带，图9（b）为图9（a）中紫色区域的放大图。

3.2 网格优劣分析

整体网格的优劣主要是由网格单元质量和网格单元数量来评估[17]，其中网格单元质量采用网格单元的最小角来评估。一般地，网格单元最小角处于[30°,60°]这一范围，定义网格单元质量较好。

在比较两种网格优劣时，有两种评判标准：

1）当网格单元数量一定时，通过网格单元最小角的分布和网格中最小角度来评估整体网格优劣；网格单元最小角分布在[30°,60°]，并且网格最小角越大，表示整体网格质量越优；

2）当网格单元质量类似时，即网格最小角大于30°，且网格单元最小角都分布在[30°,60°]，此时网格单元数量越少，说明整体网格质量越优。

基于上述两种评判标准，本文设计了两个实验来验证本文方法生成网格的优越性：

1）网格单元数量一定时，比较网格质量。

图10是网格单元数量约为5 920时，即本文和HyperMesh网格尺寸分别设为0.04m和0.02m时所生成的网格，其中图10（a）和图10（b）为本文方法生成的网格，图10（c）和图10（d）为HyperMesh生成的网格，图10（b）和图10（d）分别是图10（a）和图10（c）中紫色区域的放大图，红色区域表示最小角小于30°的畸形网格单元。

由表2和图11可以看出，本文方法生成的网格最小角为30°，网格单元的最小角大多分布在42°左右，此时整体网格质量良好；而HyperMesh生成的网格单元最小角大多分布在55°左右，但是网格最小角为 17°，虽然最小角小于 30°的网格单元数量不多，但是如图 10（d），畸形的网格单元集中存在于伞顶孔区域，由于该区域在仿真过程中属于大变形区域，则很容易引起负体积，使得计算结果发散，仿真失败。

2）网格单元质量类似时，比较网格数量。

由表3和图12可以看出，本文和HyperMesh网格尺寸分别设为0.04m和0.002m时所生成的网格，两种网格单元质量类似，此时本文方法生成的网格数量为5 922个，而HyperMesh为104 016个，比本文方法生成的网格单元数量多很多。实验结果表明，当网格质量类似时，本文方法比HyperMesh生成的网格单元数量少。

4 仿真结果

由于HyperMesh生成的伞顶孔区域的网格单元质量较差，所以在仿真求解过程中出现了负体积现象，使得结果发散。如图13所示是使用本文网格剖分器生成的网格（网格尺寸设为0.04m），数值仿真降落伞的充气过程。

开始充气时，充气形状变化比较慢。从0.06s到0.012s伞衣顶部变化显著，气流将整个伞顶部完全顶开，形成“乌贼”状。在0.019s时，充气达到最大外形形状，即充满状态。从0.019s到0.034s伞衣的充气形状发生几次较为明显的收缩。从0.034s到0.044s充气形状会不断发生微小的变化，即“呼吸现象”。最终在0.044s时充气过程基本完成。

5 结束语

本文通过建立降落伞精细几何模型，设计并实现了针对该模型的有限元网格自动剖分器，并对该网格剖分器生成的网格进行降落伞充气过程的数值仿真。实验验证：

1）当网格单元数量一定时，本文方法生成的网格不存在畸形单元，优于HyperMesh生成的网格；

2）当网格单元质量相似时，本文方法生成的网格单元数量明显少于HyperMesh生成的网格；

3）本文网格剖分器在降落伞充气过程模拟仿真中是有效的，并且表现良好。

References)

[1] 王利荣. 降落伞理论与应用[M]. 北京: 宇航出版社, 1997: 76-100. WANG Lirong. Theory and Application of Parachute[M]. Beijing: China Astronautic Publishing House, 1997: 76-100. (in Chinese)

[2] CHENG Hen. Study of Velocity Effects on Parachute Inflation Performance Based on Fluid-structure Interaction Method[J]. Applied Mathematics & Mechanics, 2014, 35(9): 1177-1188.

[3] 贾贺, 荣伟, 陈国良. 基于LS-DYNA软件的降落伞充气过程仿真研究[J]. 航天器环境工程, 2010, 27(3): 266, 367-373.

JIA He, RONG Wei, CHEN Guoliang. Simulation of Parachute Inflation Based on LS–DYNA[J]. Spacecraft Environment Engineering, 2010, 27(3): 266, 367-373. (in Chinese)

[4] 郭鹏. 大型降落伞开伞过程研究[D]. 长沙: 国防科学技术大学, 2012.

GUO Peng. Study on the Parachute Opening Process of Large Parachute[D]. Changsha: National University of Defense Technology, 2012. (in Chinese)

[5] 余莉, 史献林, 明晓. 降落伞充气过程的数值模拟[J]. 航空学报, 2007, 28(1): 52-57.

YU Li, SHI Xianlin, MING Xiao. Numerical Simulation of Parachute Inflation[J]. Chinese Journal of Aeronautics, 2007, 28(1): 52-57. (in Chinese)

[6] SHEPHARD M S. Automatic Generation and Control of Finite Element Meshes[C]. International Conference on Vehicle Structural Mechanics, Warrendale, 1988: 223-236.

[7] OWEN S J. A Survey of Unstructured Mesh Generation Technology[C]. 7th International Meshing Roundtable, Dearborn, 1998.

[8] LAWSON B C L. Software for C Surface Interpolation[C]. Mathematical Software III, Academic Press, New York, 1977: 161-194.

[9] D′AZEVEDO E F, SIMPSON R B. On Optimal Interpolation Triangle Incidences[J]. SIAM Journal Scientific Computing, 1989, 10(6), 1063-1075.

[10] MUSIN O R. Properties of the Delaunay Triangulation[C]. 13th Annual Symposium on Computational Geometry, Nice, 1997: 424-426.

[11] RAJAN V T. Optimality of the Delaunay Triangulation in Rd[C]. 7th Annual Symposium on Computational Geometry, North Conway, 1991: 357-363.

[12] KLEIN R. Construction of the Constrained Delaunay Triangulation of a Polygonal Domain[M]. CAD Systems Development, Springer Berlin Heidelberg, 1997: 313-326.

[13] MILLER G L, TALMOR D, TENG S H, et al. A Delaunay Based Numerical Method for Three Dimensions: Generation, Formulation, and Partition[J]. Antimicrobial Agents & Chemotherapy, 1997, 43(10): 2412-2416.

[14] SHEWCHUK J R. Delaunay Refinement Algorithms for Triangular Mesh Generation[J]. Computational Geometry Theory & Applications, 2010, 22(1): 21-74.

[15] 孟宪海, 李吉刚, 杨钦, 等. 复杂限定Delaunay三角化算法[J]. 中国科学: 信息科学, 2010, 40(3): 381-392.

MENG Xianhai, LI Jigang, YANG Qin, et al. Complex Delaunay Triangulation Algorithm[J]. Science China Information Sciences, 2010, 40(3): 381-392. (in Chinese)

[16] SHEWCHUK J R. Mesh Generation for Domains with Small Angles[C]. 16th ACM Symposium on Computational Geometry, Hong Kong, 2000: 1-10.

[17] 杜平安. 有限元网格划分的基本原则[J]. 机械设计与制造, 2000, (1): 34-36. DU Ping′an. Fundamental Principles of Finite Element Meshing[J]. Mechanical Design and Manufacturing, 2000, (1): 34-36. (in Chinese)

[18] SHEWCHUK J R. Delaunay Refinement Algorithms for Triangular Mesh Generation[J]. Computational Geometry, 2002, 47(7): 741-778.

[19] 徐永安, 谭建荣, 杨钦, 等. 二维任意域约束Delaunay三角化的实现[J]. 工程图学学报, 1999, 20(1): 51-55.

XU Yong′an, TAN Jianrong, YANG Qin, et al. Realization of Delaunay Triangulation in 2-D Arbitrary Domain[J]. Journal of Engineering Graphics, 1999, 20(1): 51-55. (in Chinese)

Research on Mesh Generation Technique of Parachute

WU Zhuangzhi1YANG Qian1RONG Wei2FANG Shixing2

（1 School of Computer Science and Technology, Beihang University, Beijing 100191, China）（2 Beijing Institute of Mechanics & Electricity, Beijing 100094, China）

Geometry modeling and meshing of parachute are the basis of the fluid-solid coupling numerical simulation of parachute. In this paper, the refined geometry model of parachute, including the body of canopy, radial belts, umbrella top and bottom reinforcing belts, is established, and based on software VS2010, an automatic mesh generator for this geometric model is designed and implemented, which uses C99. The generator uses constrained Delaunay algorithm, which has a good theoretical foundation and can generate high quality meshes, to meet the requirements of parachute numerical simulation. And the interface is provided to import the information of vertex and mesh into the commercial software for further simulation and analysis. Experiments show that the quality of meshes generated by this generator is superior to the quality of meshes generated by common commercial software, and the number of meshes is relatively less under similar quality of meshes. In addition the meshes behave well in numerical simulation.

mesh generation; constrained delaunay triangulation; mesh quality; parachute; spacecraft recovery

TP391, V11

: A

: 1009-8518(2017)01-0006-10

10.3969/j.issn.1009-8518.2017.01.002

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

（编辑：陈艳霞）

2016-10-16

国家重大科技专项工程