吳壯志楊倩榮偉方世興

（1 北京航空航天大學(xué)計(jì)算機(jī)學(xué)院，北京 100191）

（2 北京空間機(jī)電研究所，北京 100094）

降落傘網(wǎng)格剖分技術(shù)研究

吳壯志1楊倩1榮偉2方世興2

（1 北京航空航天大學(xué)計(jì)算機(jī)學(xué)院，北京 100191）

（2 北京空間機(jī)電研究所，北京 100094）

降落傘幾何建模和網(wǎng)格剖分是降落傘流固耦合數(shù)值仿真的基礎(chǔ)。文章建立了包括傘衣主體、徑向帶、傘頂孔和底邊加強(qiáng)帶的降落傘精細(xì)幾何模型，并基于軟件開(kāi)發(fā)平臺(tái)VS2010使用C99語(yǔ)言設(shè)計(jì)并實(shí)現(xiàn)了針對(duì)該幾何模型的網(wǎng)格自動(dòng)剖分器。該剖分器采用約束Delaunay剖分算法，具有良好的理論基礎(chǔ)、生成的網(wǎng)格質(zhì)量（quality）高，滿(mǎn)足降落傘數(shù)值仿真對(duì)網(wǎng)格質(zhì)量的要求；并提供接口可將頂點(diǎn)和網(wǎng)格信息導(dǎo)入到商用軟件中做進(jìn)一步求解分析。實(shí)驗(yàn)證明，使用文中網(wǎng)格剖分器生成的網(wǎng)格，不但網(wǎng)格單元質(zhì)量?jī)?yōu)于使用商用軟件生成的網(wǎng)格，而且在網(wǎng)格單元質(zhì)量相似的情況下，網(wǎng)格單元數(shù)量相對(duì)較少。同時(shí)，該網(wǎng)格在數(shù)值仿真中表現(xiàn)良好。

網(wǎng)格剖分 約束德洛內(nèi)三角化 網(wǎng)格質(zhì)量 降落傘 航天返回

0 引言

隨著計(jì)算機(jī)仿真技術(shù)的不斷進(jìn)步，降落傘流固耦合技術(shù)也得到了快速發(fā)展。目前降落傘開(kāi)傘過(guò)程的仿真主要借助于通用的有限元軟件來(lái)完成的[1-5]，如HyperMesh，Ls-dyna，Abaqus，Ansys等，其仿真過(guò)程包括三個(gè)部分：前處理、數(shù)值求解和后處理。曾經(jīng)有人做過(guò)統(tǒng)計(jì)，在數(shù)值模擬仿真的三個(gè)階段中，前處理約占總時(shí)間的40%～60%，數(shù)值求解約占5%～20%，后處理約占30%[6]。前處理過(guò)程中的幾何建模和有限元網(wǎng)格剖分工作比較繁瑣，但是又十分重要，是進(jìn)行數(shù)值模擬仿真的基礎(chǔ)。

目前在對(duì)降落傘進(jìn)行幾何建模和網(wǎng)格剖分時(shí)，仍然有幾點(diǎn)不足之處：

第一，通常采用 CAD軟件對(duì)降落傘進(jìn)行幾何建模，雖然降落傘的形狀大小都比較類(lèi)似，但是只要形狀或尺寸有一點(diǎn)不同的地方，都需要重復(fù)建模過(guò)程，非常費(fèi)時(shí)費(fèi)力。

第二，通用的商用有限元軟件生成的網(wǎng)格單元是整齊劃一的，而仿真過(guò)程中傘頂孔位置屬于大變形區(qū)域，需要生成較稠密的網(wǎng)格。為了滿(mǎn)足這一需求，通常有兩種做法：一是將整體網(wǎng)格密度加大，達(dá)到大變形區(qū)域的網(wǎng)格密度，但是這樣會(huì)花費(fèi)成倍的計(jì)算時(shí)間，浪費(fèi)資源；二是人工分區(qū)域分別進(jìn)行不同網(wǎng)格密度的剖分，雖然可以達(dá)到對(duì)網(wǎng)格密度的要求，但是費(fèi)時(shí)費(fèi)力。

第三，真實(shí)的降落傘結(jié)構(gòu)包括傘衣主體、徑向帶，傘頂孔和底邊加強(qiáng)帶等，但是通常只針對(duì)傘衣主體進(jìn)行建模和網(wǎng)格剖分，簡(jiǎn)化徑向帶、傘頂孔和底邊加強(qiáng)帶為一維梁?jiǎn)卧猍3]（Beam）。

針對(duì)上述三點(diǎn)不足，本文建立了包括傘衣主體、徑向帶、傘頂孔和底邊加強(qiáng)帶的降落傘精細(xì)的幾何模型，并對(duì)降落傘的網(wǎng)格剖分技術(shù)進(jìn)行研究，基于軟件開(kāi)發(fā)平臺(tái)VS2010使用C99語(yǔ)言設(shè)計(jì)并實(shí)現(xiàn)了一個(gè)針對(duì)該幾何模型的網(wǎng)格自動(dòng)剖分器。該剖分器具有如下特點(diǎn)：采用約束Delaunay剖分算法，具有良好的理論基礎(chǔ)、生成的網(wǎng)格質(zhì)量（quality，下同）高，滿(mǎn)足降落傘數(shù)值仿真對(duì)網(wǎng)格質(zhì)量的要求；提供接口可以將生成的頂點(diǎn)和網(wǎng)格信息轉(zhuǎn)換成K文件，導(dǎo)入到HyperMesh等通用的商用軟件中做進(jìn)一步分析計(jì)算。

1 網(wǎng)格剖分相關(guān)理論

由于非結(jié)構(gòu)化網(wǎng)格可以很好適應(yīng)復(fù)雜自由邊界，并且可以很容易控制網(wǎng)格尺寸，所以本文將對(duì)降落傘生成非結(jié)構(gòu)化網(wǎng)格。典型通用的非結(jié)構(gòu)化網(wǎng)格生成方法有：四叉樹(shù)法、前沿推進(jìn)法和Delaunay三角剖分算法。文獻(xiàn)[7]曾對(duì)用于商業(yè)和科研項(xiàng)目的81種網(wǎng)格生成軟件產(chǎn)品進(jìn)行了非正式的調(diào)查，其中37種軟件使用的是 Delaunay算法。由此可以看出，Delaunay算法占據(jù)著重要的位置。所以本文擬采用約束Delaunay算法（Constrained Delaunay Triangulation）對(duì)降落傘進(jìn)行網(wǎng)格剖分。

1.1 Delaunay三角化

首次給出Delaunay三角網(wǎng)格概念的是俄國(guó)著名數(shù)學(xué)家，Boris Nikolaevich Delaunay。這里給出其精確定義：

設(shè)t是三角網(wǎng)格T的一個(gè)網(wǎng)格單元，若t的外接圓的內(nèi)部不包含任何頂點(diǎn)，則稱(chēng)t是Delaunay的，也稱(chēng)其滿(mǎn)足空?qǐng)A特性，如圖1所示。設(shè)P為二維空間R2中的有限點(diǎn)集，T為P的三角網(wǎng)格，若T中任何一個(gè)三角形t都是Delaunay的，則稱(chēng)T是P的一個(gè)Delaunay三角化，記為DT(P)。

基于以上定義可知，Delaunay三角化最重要的性質(zhì)就是空?qǐng)A特性。為了使三角網(wǎng)格單元滿(mǎn)足空?qǐng)A特性，Lawson提出了局部?jī)?yōu)化過(guò)程[8]（Local Optimization Procedure，LOP）。另外，Delaunay三角化還具有一些其他的優(yōu)良性質(zhì)，這里就不一一贅述[9-11]。

1.2 約束Delaunay算法

約束Delaunay算法，顧名思義是在Delaunay算法的基礎(chǔ)上加了一些約束條件，包括約束點(diǎn)集V和約束邊集 E，嚴(yán)格要求約束條件中的每條約束邊是三角剖分后某些三角形網(wǎng)格單元的邊，則三角網(wǎng)格應(yīng)滿(mǎn)足約束Delaunay屬性[12]。

設(shè)點(diǎn)v′為任一點(diǎn)，t是一個(gè)三角網(wǎng)格單元，e′為約束邊集E中的線(xiàn)段，若v′與t內(nèi)部任何一點(diǎn)的連線(xiàn)都不與e′相交，則稱(chēng)點(diǎn)v′是對(duì)于t內(nèi)部可見(jiàn)的點(diǎn)。如圖2所示，w′為頂點(diǎn)，其中點(diǎn)v′對(duì)t是可見(jiàn)的，點(diǎn)w′對(duì)t是不可見(jiàn)的。

設(shè)t是T中的一個(gè)網(wǎng)格單元，若t的外接圓中不包含任何對(duì)于t可見(jiàn)的點(diǎn)，則稱(chēng)t是約束Delaunay的，也稱(chēng)約束Delaunay屬性。

1.3 網(wǎng)格單元質(zhì)量評(píng)價(jià)指標(biāo)

在計(jì)算幾何中，Miller等[13]指出，分析一個(gè)三角網(wǎng)格單元質(zhì)量最好的指標(biāo)是單元最小角、外接圓半徑與最短邊之比r1和縱橫比r2。一個(gè)三角網(wǎng)格單元t的縱橫比r2為t的外接圓半徑R與內(nèi)切圓半徑r之比。其中最直觀(guān)的就是三角網(wǎng)格單元的最小角，如果最小角越小，則該三角形越畸形。

2 降落傘的網(wǎng)格剖分

2.1 降落傘網(wǎng)格剖分特點(diǎn)

本文以某平面圓形傘為例，對(duì)該傘進(jìn)行精細(xì)結(jié)構(gòu)建模和網(wǎng)格剖分的研究。其中傘衣結(jié)構(gòu)包括傘衣主體、徑向帶、傘頂孔和底邊加強(qiáng)帶，結(jié)構(gòu)如圖3所示。傘衣有12幅，徑向帶、傘頂孔和底邊加強(qiáng)帶的個(gè)數(shù)與傘衣個(gè)數(shù)相同，每幅傘衣都添加了一條邊，即折痕，為之后的折疊做準(zhǔn)備。具體的結(jié)構(gòu)參數(shù)如表1所示，其中徑向帶寬度與加強(qiáng)帶寬度值一樣。

表1 降落傘結(jié)構(gòu)參數(shù)表Tab.1 Structure parameters of a parachute

降落傘共包括N0幅傘衣，每幅傘衣精細(xì)的幾何結(jié)構(gòu)如圖4所示（第i幅傘衣），主要包括：傘衣主體（四邊形 v11v13v14v10），折痕（線(xiàn)段v1v8），傘頂孔加強(qiáng)帶（四邊形 v10v14v7v+）、底邊加強(qiáng)帶（四邊形 v0v2v13v11）和徑向帶（四邊形 v2v3v6v7）。將N0幅傘衣合成一體得到整個(gè)降落傘傘衣的幾何模型。

2.2 網(wǎng)格剖分流程

本文擬采用約束 Delaunay算法對(duì)降落傘精細(xì)幾何模型進(jìn)行有限元網(wǎng)格剖分，算法描述如下：

輸入：降落傘的約束點(diǎn)集V和約束邊集E，網(wǎng)格尺寸控制指標(biāo)S，網(wǎng)格質(zhì)量控制指標(biāo)B。

輸出：降落傘網(wǎng)格剖分DT。

算法：1）根據(jù)輸入的約束點(diǎn)集合V，生成初始Delaunay剖分DT0；2）檢查約束邊集E是否在初始網(wǎng)格剖分DT0中，恢復(fù)E中丟失的約束邊集E′，得到網(wǎng)格剖分DT1；3）對(duì)恢復(fù)約束邊集后的網(wǎng)格DT1調(diào)用約束Delaunay優(yōu)化算法，得到最終的網(wǎng)格剖分DT；4）保存頂點(diǎn)和網(wǎng)格的信息。

2.3 初始Delaunay剖分

本文使用逐點(diǎn)插入的Lawson算法對(duì)約束點(diǎn)集V進(jìn)行初始Delaunay剖分?；驹頌椋菏紫冉⒁粋€(gè)大的三角形或多邊形，把所有約束點(diǎn)包圍起來(lái)，取V中一點(diǎn)，該點(diǎn)與包含它的三角形三個(gè)頂點(diǎn)相連，形成三個(gè)新的三角形，然后逐個(gè)對(duì)它們用LOP進(jìn)行優(yōu)化，以保證所生成的三角網(wǎng)格滿(mǎn)足Delaunay的空?qǐng)A屬性。

經(jīng)過(guò)上述步驟可得到包括約束點(diǎn)集V的初始Delaunay剖分DT0，如圖5所示。但是DT0可能不完全包括約束邊集E，對(duì)于丟失的約束邊集E′，需要對(duì)其進(jìn)行恢復(fù)。

2.4 約束邊的恢復(fù)

通過(guò)在丟失的約束邊中點(diǎn)處添加頂點(diǎn)，并經(jīng)過(guò)LOP局部?jī)?yōu)化，恢復(fù)丟失的約束邊E′。如果還存在丟失的子約束邊，則繼續(xù)迭代地在其中點(diǎn)處添加頂點(diǎn)，如圖6所示，得到恢復(fù)約束邊集后的網(wǎng)格DT1。約束邊e可能被分割為k條子約束邊 e1, e2, e3,… ,ek，約束邊集E中將不包含e，而包含 e1, e2, e3,… ,ek。

2.5 約束Delaunay優(yōu)化算法

經(jīng)過(guò)約束邊恢復(fù)的網(wǎng)格一般都很粗糙，而且比較細(xì)長(zhǎng)，需要經(jīng)過(guò)加密、勻稱(chēng)等網(wǎng)格優(yōu)化處理才能用于數(shù)值分析。網(wǎng)格優(yōu)化[14]是整個(gè)算法的核心，其中關(guān)鍵部分為：

1）如何控制生成的三角網(wǎng)格單元的質(zhì)量：一般地，通過(guò)設(shè)置r1的上界B來(lái)控制網(wǎng)格單元的質(zhì)量，本文設(shè)置B為1，即除輸入條件中的尖角，生成的網(wǎng)格單元內(nèi)角都會(huì)控制在[30°,120°]，很好的控制了三角網(wǎng)格單元的質(zhì)量；

2）如何控制生成的三角網(wǎng)格單元的尺寸：通常設(shè)置三角網(wǎng)格單元邊長(zhǎng)L的上界Lmax來(lái)控制網(wǎng)格單元的尺寸，即所有網(wǎng)格單元的邊長(zhǎng)都須小于Lmax；

3）如何控制優(yōu)化算法的終止：約束Delaunay優(yōu)化算法是通過(guò)不斷遍歷三角網(wǎng)格單元，對(duì)網(wǎng)格單元質(zhì)量和尺寸不符合要求的網(wǎng)格單元進(jìn)行優(yōu)化。這是一個(gè)迭代優(yōu)化的過(guò)程。

2.5.1 頂點(diǎn)生成方法

不論是控制三角網(wǎng)格單元的質(zhì)量還是控制三角網(wǎng)格單元的尺寸，都需要通過(guò)添加新的頂點(diǎn)來(lái)完成。新頂點(diǎn)的生成需要遵循兩個(gè)原則：

1）以子約束邊為直徑的圓是包含子約束邊的最小圓，若該圓內(nèi)部包含其他頂點(diǎn)，則在該邊中點(diǎn)處添加新的頂點(diǎn)，迭代添加頂點(diǎn)，直到各子約束邊的直徑圓內(nèi)部不包含其他頂點(diǎn)；

2）如果r1＞B或邊長(zhǎng)大于Lmax，則說(shuō)明該網(wǎng)格單元是畸形的，或該網(wǎng)格單元尺寸過(guò)大。那么在該網(wǎng)格單元的外接圓圓心處添加新的頂點(diǎn)，并去除該網(wǎng)格單元。然而，如果新頂點(diǎn)在某子約束邊的直徑圓內(nèi)部，則不插入該頂點(diǎn)，并將相應(yīng)的子約束邊在中點(diǎn)處分割。該算法描述為：

輸入：網(wǎng)格DT，約束邊集合E；

輸出：不包含畸形和尺寸過(guò)大的網(wǎng)格單元的網(wǎng)格剖分DT和新的約束邊集合E；

遵循上述兩個(gè)原則，一方面，約束插入的頂點(diǎn)符合約束Delaunay準(zhǔn)則；另一方面，設(shè)置r1的上界B和邊長(zhǎng)的上界Lmax，很好地控制了網(wǎng)格單元的質(zhì)量和尺寸。

2.5.2 輸入條件中尖角的處理方法

對(duì)于尖角的處理方法，根據(jù)原則二插入其所在三角形的外接圓圓心來(lái)處理。但是若輸入的約束條件中存在尖角[15-16]，根據(jù)原則二，會(huì)不斷插入新的頂點(diǎn)，進(jìn)入死循環(huán)，如圖 7所示。對(duì)于這種情況，可按如下方法來(lái)處理。

若該頂點(diǎn)的插入半徑rv小于其父頂點(diǎn)的插入半徑 rpv時(shí)，停止繼續(xù)插入新的頂點(diǎn)。

其中關(guān)于插入半徑和父頂點(diǎn)的定義如下。根據(jù)插入半徑分為四種情況：若頂點(diǎn)v為約束點(diǎn)，則其插入半徑rv為v到對(duì)v可見(jiàn)的最近約束點(diǎn)的歐氏距離；若v為某約束邊的中點(diǎn)，則rv為v到距其最近的網(wǎng)格頂點(diǎn)的距離；若v為被拒絕插入的某三角形外接圓圓心，則rv為該約束邊長(zhǎng)的一半；若v為某三角形外接圓圓心，則rv為外接圓半徑。

輸入條件中的約束點(diǎn)沒(méi)有父頂點(diǎn)，其余每個(gè)頂點(diǎn) v，包括被拒絕插入的頂點(diǎn)，都有父頂點(diǎn)pv。若 v為某約束邊的中點(diǎn)，則pv為該約束邊任一頂點(diǎn)；若v為插入或被拒絕的某三角形的外接圓圓心，則pv為該三角形最短邊的兩頂點(diǎn)中最后插入的頂點(diǎn)。

那么，根據(jù)上述條件，約束邊分割算法為：

輸入：三角網(wǎng)格DT，約束邊集合E。

輸出：子約束邊直徑圓中不包含頂點(diǎn)的三角網(wǎng)格DT和新的約束邊集合E。

算法：1）for e∈Edo；2）while e的直徑圓內(nèi)存在其他頂點(diǎn)將該邊e進(jìn)行分割，形成新的三角網(wǎng)格；更新約束邊集合E；4）調(diào)用LOP優(yōu)化算法；5）endwhile；6）endfor。

2.5.3 終止條件

約束Delaunay優(yōu)化是一個(gè)迭代的過(guò)程，優(yōu)化算法的終止條件為所有網(wǎng)格單元t滿(mǎn)足：；②邊長(zhǎng)。

綜上所述，約束Delaunay優(yōu)化算法可表述為：

輸入：三角網(wǎng)格DT，約束邊集合E。

輸出：優(yōu)化后的約束Delaunay三角網(wǎng)格DT。

算法：1）while不滿(mǎn)足條件①②do；2）調(diào)用約束邊分割算法；3）調(diào)用畸形網(wǎng)格單元調(diào)整算法；4）endwhile；5）保存所有頂點(diǎn)和網(wǎng)格信息。

2.6 網(wǎng)格輸出處理

Keyword文件簡(jiǎn)稱(chēng)K文件，是商用有限元軟件可接收的文件格式之一。K文件就是在一組關(guān)鍵詞組織下的數(shù)據(jù)塊組成的。關(guān)鍵詞必須包含符號(hào)“*”作為標(biāo)識(shí)，“*KEYWORD”和“*END”分別是K文件的開(kāi)始符和結(jié)束符。

常用的關(guān)鍵詞包括：*NODE，*ELEMENT，*PART等。*NODE表示頂點(diǎn)信息，NID表示頂點(diǎn)編號(hào)，X，Y，Z表示頂點(diǎn)坐標(biāo)。*ELEMENT表示網(wǎng)格信息，EID是其編號(hào)，PID表示*PART編號(hào)，N1，N2，N3分別表示網(wǎng)格單元頂點(diǎn)的NID。圖8展示了K文件的組織方式以及各種實(shí)體之間的關(guān)系。

本文根據(jù)*NODE和*ELEMENT內(nèi)容以及編號(hào)、坐標(biāo)等格式組織網(wǎng)格信息，最后生成的文件后綴設(shè)為“.k”，即為K文件，可導(dǎo)入HyperMesh等商用有限元軟件做進(jìn)一步仿真分析。

3 網(wǎng)格優(yōu)劣評(píng)價(jià)

本文以某平面圓形傘為例，建立其精細(xì)幾何模型，并生成三角網(wǎng)格。作為對(duì)比試驗(yàn)，使用HyperMesh軟件對(duì)該幾何模型生成三角網(wǎng)格。

3.1 網(wǎng)格剖分結(jié)果

圖9是采用本文方法對(duì)整個(gè)降落傘生成的網(wǎng)格，其中黑色部分為傘衣，藍(lán)色部分為加強(qiáng)帶，圖9（b）為圖9（a）中紫色區(qū)域的放大圖。

3.2 網(wǎng)格優(yōu)劣分析

整體網(wǎng)格的優(yōu)劣主要是由網(wǎng)格單元質(zhì)量和網(wǎng)格單元數(shù)量來(lái)評(píng)估[17]，其中網(wǎng)格單元質(zhì)量采用網(wǎng)格單元的最小角來(lái)評(píng)估。一般地，網(wǎng)格單元最小角處于[30°,60°]這一范圍，定義網(wǎng)格單元質(zhì)量較好。

在比較兩種網(wǎng)格優(yōu)劣時(shí)，有兩種評(píng)判標(biāo)準(zhǔn)：

1）當(dāng)網(wǎng)格單元數(shù)量一定時(shí)，通過(guò)網(wǎng)格單元最小角的分布和網(wǎng)格中最小角度來(lái)評(píng)估整體網(wǎng)格優(yōu)劣；網(wǎng)格單元最小角分布在[30°,60°]，并且網(wǎng)格最小角越大，表示整體網(wǎng)格質(zhì)量越優(yōu)；

2）當(dāng)網(wǎng)格單元質(zhì)量類(lèi)似時(shí)，即網(wǎng)格最小角大于30°，且網(wǎng)格單元最小角都分布在[30°,60°]，此時(shí)網(wǎng)格單元數(shù)量越少，說(shuō)明整體網(wǎng)格質(zhì)量越優(yōu)。

基于上述兩種評(píng)判標(biāo)準(zhǔn)，本文設(shè)計(jì)了兩個(gè)實(shí)驗(yàn)來(lái)驗(yàn)證本文方法生成網(wǎng)格的優(yōu)越性：

1）網(wǎng)格單元數(shù)量一定時(shí)，比較網(wǎng)格質(zhì)量。

圖10是網(wǎng)格單元數(shù)量約為5 920時(shí)，即本文和HyperMesh網(wǎng)格尺寸分別設(shè)為0.04m和0.02m時(shí)所生成的網(wǎng)格，其中圖10（a）和圖10（b）為本文方法生成的網(wǎng)格，圖10（c）和圖10（d）為HyperMesh生成的網(wǎng)格，圖10（b）和圖10（d）分別是圖10（a）和圖10（c）中紫色區(qū)域的放大圖，紅色區(qū)域表示最小角小于30°的畸形網(wǎng)格單元。

由表2和圖11可以看出，本文方法生成的網(wǎng)格最小角為30°，網(wǎng)格單元的最小角大多分布在42°左右，此時(shí)整體網(wǎng)格質(zhì)量良好；而HyperMesh生成的網(wǎng)格單元最小角大多分布在55°左右，但是網(wǎng)格最小角為 17°，雖然最小角小于 30°的網(wǎng)格單元數(shù)量不多，但是如圖 10（d），畸形的網(wǎng)格單元集中存在于傘頂孔區(qū)域，由于該區(qū)域在仿真過(guò)程中屬于大變形區(qū)域，則很容易引起負(fù)體積，使得計(jì)算結(jié)果發(fā)散，仿真失敗。

2）網(wǎng)格單元質(zhì)量類(lèi)似時(shí)，比較網(wǎng)格數(shù)量。

由表3和圖12可以看出，本文和HyperMesh網(wǎng)格尺寸分別設(shè)為0.04m和0.002m時(shí)所生成的網(wǎng)格，兩種網(wǎng)格單元質(zhì)量類(lèi)似，此時(shí)本文方法生成的網(wǎng)格數(shù)量為5 922個(gè)，而HyperMesh為104 016個(gè)，比本文方法生成的網(wǎng)格單元數(shù)量多很多。實(shí)驗(yàn)結(jié)果表明，當(dāng)網(wǎng)格質(zhì)量類(lèi)似時(shí)，本文方法比HyperMesh生成的網(wǎng)格單元數(shù)量少。

4 仿真結(jié)果

由于HyperMesh生成的傘頂孔區(qū)域的網(wǎng)格單元質(zhì)量較差，所以在仿真求解過(guò)程中出現(xiàn)了負(fù)體積現(xiàn)象，使得結(jié)果發(fā)散。如圖13所示是使用本文網(wǎng)格剖分器生成的網(wǎng)格（網(wǎng)格尺寸設(shè)為0.04m），數(shù)值仿真降落傘的充氣過(guò)程。

開(kāi)始充氣時(shí)，充氣形狀變化比較慢。從0.06s到0.012s傘衣頂部變化顯著，氣流將整個(gè)傘頂部完全頂開(kāi)，形成“烏賊”狀。在0.019s時(shí)，充氣達(dá)到最大外形形狀，即充滿(mǎn)狀態(tài)。從0.019s到0.034s傘衣的充氣形狀發(fā)生幾次較為明顯的收縮。從0.034s到0.044s充氣形狀會(huì)不斷發(fā)生微小的變化，即“呼吸現(xiàn)象”。最終在0.044s時(shí)充氣過(guò)程基本完成。

5 結(jié)束語(yǔ)

本文通過(guò)建立降落傘精細(xì)幾何模型，設(shè)計(jì)并實(shí)現(xiàn)了針對(duì)該模型的有限元網(wǎng)格自動(dòng)剖分器，并對(duì)該網(wǎng)格剖分器生成的網(wǎng)格進(jìn)行降落傘充氣過(guò)程的數(shù)值仿真。實(shí)驗(yàn)驗(yàn)證：

1）當(dāng)網(wǎng)格單元數(shù)量一定時(shí)，本文方法生成的網(wǎng)格不存在畸形單元，優(yōu)于HyperMesh生成的網(wǎng)格；

2）當(dāng)網(wǎng)格單元質(zhì)量相似時(shí)，本文方法生成的網(wǎng)格單元數(shù)量明顯少于HyperMesh生成的網(wǎng)格；

3）本文網(wǎng)格剖分器在降落傘充氣過(guò)程模擬仿真中是有效的，并且表現(xiàn)良好。

References)

[1] 王利榮. 降落傘理論與應(yīng)用[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] 賈賀, 榮偉, 陳國(guó)良. 基于LS-DYNA軟件的降落傘充氣過(guò)程仿真研究[J]. 航天器環(huán)境工程, 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] 郭鵬. 大型降落傘開(kāi)傘過(guò)程研究[D]. 長(zhǎng)沙: 國(guó)防科學(xué)技術(shù)大學(xué), 2012.

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

[5] 余莉, 史獻(xiàn)林, 明曉. 降落傘充氣過(guò)程的數(shù)值模擬[J]. 航空學(xué)報(bào), 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] 孟憲海, 李吉?jiǎng)? 楊欽, 等. 復(fù)雜限定Delaunay三角化算法[J]. 中國(guó)科學(xué): 信息科學(xué), 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] 杜平安. 有限元網(wǎng)格劃分的基本原則[J]. 機(jī)械設(shè)計(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三角化的實(shí)現(xiàn)[J]. 工程圖學(xué)學(xué)報(bào), 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年獲北京航空航天大學(xué)計(jì)算機(jī)軟件與理論博士學(xué)位，副教授。研究方向?yàn)橛?jì)算機(jī)圖形學(xué)、計(jì)算幾何和人體測(cè)量學(xué)等。E-mail：zzwu@buaa.edu.cn。

（編輯：陳艷霞）

2016-10-16

國(guó)家重大科技專(zhuān)項(xiàng)工程