王佳穎,劉星雨

(1.武警工程大學(xué) 基礎(chǔ)部， 西安 710086；2.武警工程大學(xué) 裝備管理與保障學(xué)院， 西安 710086；3.上海交通大學(xué)海洋智能裝備與系統(tǒng)教育部實(shí)驗(yàn)室, 上海 200240)

對(duì)于不規(guī)則形狀破片迎風(fēng)面積(投影區(qū)域面積)的計(jì)算，通常采用二十面體均勻取向法的光學(xué)投影系統(tǒng)試驗(yàn)[1]。該參數(shù)可進(jìn)一步用于破片速度衰減分析、飛散特性研究以及比動(dòng)能評(píng)估[3-7]。為便于自動(dòng)化計(jì)算有限元破片迎風(fēng)面積并提高求解精度，根據(jù)蒙特卡洛剖分投影法建立的平均迎風(fēng)面積計(jì)算模型[8]提出一種最大邊長約束的凹域三角剖分算法。該算法以凸包Delaunay三角剖分算法為基礎(chǔ)[9-11]，以最大邊長約束為原則，對(duì)凹多邊形區(qū)域(簡稱凹域)邊界進(jìn)行搜索重構(gòu)，可進(jìn)一步提高平均迎風(fēng)面積計(jì)算模型的求解精度。目前基于最大邊長約束的凹域三角剖分自動(dòng)化求解破片迎風(fēng)面積方法，未見相關(guān)報(bào)道。雖然凸包Delaunay三角剖分算法也可直接用于破片迎風(fēng)面積的計(jì)算，但對(duì)于形狀不可預(yù)知的自然破片而言其誤差仍有待進(jìn)一步分析。

本文根據(jù)最大邊長約束的凹域三角剖分算法，通過凸包Delaunay三角剖分、剖分邊界凹凸性判別、最大邊長約束的剖分去除、剖分邊界邊搜索、邊界邊排序等步驟，得到任意形狀破片投影邊界并用于平均迎風(fēng)面積模型計(jì)算。同時(shí)，驗(yàn)證并對(duì)比了最大邊長約束的凹域三角剖分算法與凸包Delaunay三角剖分算法所得平均迎風(fēng)面積的精度與求解時(shí)間。

1 有限元破片

本文以37 mm爆震彈為例[12]，模型采用1/4對(duì)稱建模，網(wǎng)格大小1 mm×1 mm×2 mm，殼體破碎情況如圖1(a)所示。在40 μs時(shí)刻殼體已破碎完全，將該時(shí)刻各破片單元節(jié)點(diǎn)數(shù)據(jù)導(dǎo)出進(jìn)行投影剖分。為對(duì)比最大邊長約束的凹域三角剖分算法與凸包Delaunay三角剖分的迎風(fēng)面積計(jì)算精度與效率，本文重點(diǎn)選擇頂部方向、徑向、底部方向形狀差異較大破片進(jìn)行求解分析，其具體尺寸形狀如圖1(b)、圖1(c)、圖1(d)所示。

2 最大邊長約束的凹域三角剖分算法

2.1 算法概述

最大邊長約束的凹域三角剖分算法主要包括凸包Delaunay三角剖分、剖分邊界凹凸性判別、最大邊長約束的剖分去除、凹域邊界邊搜索、邊界邊排序5個(gè)步驟，其算法流程如圖2所示。

圖1 37 mm爆震彈破碎情況及不同方向上的有限元破片

圖2 最大邊長約束的凹域三角剖分算法流程框圖

凸包Delaunay三角剖分根據(jù)空外接圓規(guī)則可將離散投影節(jié)點(diǎn)集T構(gòu)建成Delaunay三角網(wǎng)，并獲得凸包邊界dm及三角形矩陣dt_clist；剖分邊界凹凸性判別根據(jù)凸包三角剖分是否存在大于R的邊為原則進(jìn)行遍歷，可提高破片投影為凸包時(shí)的計(jì)算效率；最大邊長約束的剖分去除以不大于凸包三角剖分平均邊長的λscale倍為原則，將超出約束邊長R的三角邊剔除；凹域邊界搜索根據(jù)公共邊、所有邊的差集運(yùn)算，可得到邊界邊集合；邊界邊排序根據(jù)進(jìn)棧、出棧的搜索遍歷，可得到經(jīng)排序后的多邊形頂點(diǎn)。而后采用叉乘法對(duì)任意多邊形面積進(jìn)行求解。

其中，T為無重復(fù)節(jié)點(diǎn)的投影點(diǎn)集；dm為凸包剖分邊界；dt_Clist為凸包三角剖分三角形節(jié)點(diǎn)編號(hào)；edges為凸包三角剖分升序排列的邊界；I為edges中不包含重復(fù)邊的行號(hào)集合；IN為edges的行號(hào)集合；IO為邊界邊行號(hào)；sort_vertices用于存放邊界頂點(diǎn)序列；polygon為按順(逆)時(shí)針排序的邊界邊節(jié)點(diǎn)坐標(biāo)。

R取值定義如下：

(1)

其中，R為邊長約束值；λscale為放大因子；li為剖分三角形第i條邊的長度；ne為三角邊總數(shù)。根據(jù)凸包Delaunay剖分空外接圓準(zhǔn)則及正六面體有限元節(jié)點(diǎn)距離特征，λscale取8～10時(shí)可較好地除去非凹域點(diǎn)集的三角邊，同時(shí)也可避免刪除凹域點(diǎn)集內(nèi)剖分長度較大的有效三角邊。本文λscale取8。

2.2 凸包Delaunay三角剖分

當(dāng)前隨機(jī)點(diǎn)集的Delaunay三角剖分方法主要有逐點(diǎn)插入法、分治算法、三角組網(wǎng)法等[11]。本文采用改進(jìn)Bowyer-Watson逐點(diǎn)插入算法，該算法相比于Lawson算法計(jì)算效率更高；相比于標(biāo)準(zhǔn)Bowyer-Watson算法，魯棒性更好。改進(jìn)的Bowyer-Watson算法通過引入放大因子，可避免標(biāo)準(zhǔn)超級(jí)三角形為大鈍角時(shí)外接圓心落在其外造成剖分邊界非凸包的情況。通過判定插入點(diǎn)在待剖分三角形右側(cè)、內(nèi)側(cè)或外側(cè)位置，分別對(duì)剖分三角形進(jìn)行入棧、再剖分、加入下次剖分等操作，最終可形成滿足空外接圓特性和最小角最大化特性的Delaunay三角剖分網(wǎng)格。其算法流程如圖3所示。

圖3 凸包Delaunay三角剖分算法流程框圖

判斷某點(diǎn)坐標(biāo)(xi,yi)與三角形外接圓位置，可先根據(jù)三角形三頂點(diǎn)坐標(biāo)(xi-1,yi-1)、(xi-2,yi-2)、(xi-3,yi-3)，求得外接圓圓心坐標(biāo)(xo,yo)及半徑r，而后通過該點(diǎn)與圓心的距離以及與半徑之間的關(guān)系判定該點(diǎn)位置在圓內(nèi)側(cè)、圓外側(cè)、圓右側(cè)。外接圓圓心坐標(biāo)可根據(jù)圓上三點(diǎn)至圓心距離兩兩相等推導(dǎo)得到，表達(dá)式如式(2)、式(3)所示。點(diǎn)(xi,yi)與外接圓位置判定表達(dá)式如式(4)所示。

(2)

(3)

(4)

2.3 凹域邊界搜索排序及編程實(shí)現(xiàn)

凹域邊界搜索的差集運(yùn)算表達(dá)式如式(5)：

EIO=EI-EIN-I

(5)

其中，EIO為凹多邊形區(qū)域邊界邊集合，有且僅有一條公共邊；EIN-I為凹多邊形區(qū)域的內(nèi)部邊集合，存在兩個(gè)公共邊；EI為不包含重復(fù)邊的凹多邊形邊集合。

凹域邊界邊搜索與排序的代碼如下：

edges數(shù)據(jù)結(jié)構(gòu)為一行存放一條角邊，連續(xù)三行為同一三角形；IA為edges中唯一值邊對(duì)應(yīng)的行號(hào)；通過差集運(yùn)算可獲得邊界邊點(diǎn)集并存入edges中；sort_vertices擬存放排序后的多邊形頂點(diǎn)序號(hào)，比當(dāng)前edges行數(shù)多1，確保首尾頂點(diǎn)相連。最終經(jīng)排序后的凹域邊界頂點(diǎn)坐標(biāo)存入polygon中。

3 結(jié)果與討論

3.1 凹域剖分算法的可行性驗(yàn)證

為驗(yàn)證最大邊長約束的凹域三角剖分算法可行性，對(duì)圖1中的頂部破片、徑向破片及底部破片進(jìn)行求解分析。經(jīng)過一次投影、凸包三角剖分、帶最大邊長約束的凹域剖分結(jié)果分別如圖4(a)、(b)、(c)所示。從圖5可知三枚破片形狀尺寸差異較大，但通過最大邊長約束的凹域三角剖分算法均可得到破片投影圖形的邊界，并與投影點(diǎn)集吻合較好；凸包Delaunay三角剖分所得邊界比投影點(diǎn)集真實(shí)邊界面積更大。在單次投影結(jié)果中，圖4(a)中凸包剖分所得面積與凹域剖分面積相差不大，而在圖4(b)、(c)中與凹域剖分所得結(jié)果則相差較大。由此說明，相比于凸包Delaunay三角剖分算法，最大邊長約束的凹域三角剖分算法能更準(zhǔn)確地對(duì)任意形狀破片投影邊界進(jìn)行搜索，可進(jìn)一步提高自然破片平均迎風(fēng)面積的求解精度。

3.2 平均迎風(fēng)面積精度對(duì)比

對(duì)比3枚破片5 000次投影后所得平均迎風(fēng)面積結(jié)果，兩種算法下的迎風(fēng)面積均值、標(biāo)準(zhǔn)差見表1。

圖4 不同方向破片的一次投影、凸包三角剖分、帶最大邊長約束的凹域剖分結(jié)果圖

表1 凸包Delaunay算法與最大邊長約束的凹域剖分算法均值和標(biāo)準(zhǔn)差

從表1可知：凸包Delaunay算法頂部破片、底部破片標(biāo)準(zhǔn)差均比凹域剖分算法要小，但徑向破片標(biāo)準(zhǔn)差更大；凸包Delaunay算法不同方向破片所得迎風(fēng)面積均值均比凹域剖分算法要大，其中頂部破片、底部破片均值誤差相差不大，僅為5.4%和3.6%的差異，但徑向破片均值誤差可達(dá)66.5%。分析造成凸包Delaunay算法各破片均值誤差差異較大的原因，主要由各破片形狀特征量單元節(jié)點(diǎn)與質(zhì)心距離的標(biāo)準(zhǔn)差決定。對(duì)于節(jié)點(diǎn)與質(zhì)心距離標(biāo)準(zhǔn)差較小的頂部、底部破片，其投影圖形更加緊湊，使得非凹域點(diǎn)集的無效三角邊更少，誤差相對(duì)變??；而對(duì)于距離標(biāo)準(zhǔn)差更大的徑向破片，其投影圖形更加狹長，使得非凹域點(diǎn)集的無效三角邊增多，進(jìn)而誤差顯著增大。由此說明，凸包Delaunay算法的平均迎風(fēng)面積精度受到破片形狀特征量影響較大，對(duì)于節(jié)點(diǎn)與質(zhì)心距離標(biāo)準(zhǔn)差更小的破片而言凸包算法與凹域剖分算法精度相差較小，而對(duì)于節(jié)點(diǎn)與質(zhì)心距離標(biāo)準(zhǔn)差較大的破片精度則相差較大。對(duì)自然破片的平均迎風(fēng)面積求解，采用最大邊長約束的凹域剖分算法更能保證計(jì)算精度。

3.3 平均迎風(fēng)面積求解時(shí)間對(duì)比

為避免CPU時(shí)鐘的誤差影響，選擇100次循環(huán)所得平均時(shí)間作為求解時(shí)間，兩種算法計(jì)算得到的時(shí)間見表2。運(yùn)行環(huán)境為Intel (R) Core (TM) i7-8750h CPU @ 2.20 GHz，16 GB，Win10 x64。從表2可知：對(duì)于同一算法而言，不同方向的各破片求解時(shí)間相差不大，說明兩種算法的破片單元數(shù)對(duì)求解時(shí)間并不敏感；而對(duì)于同一破片而言，最大邊長約束的凹域剖分算法比凸包Delaunay剖分算法求解時(shí)間更長，約為30倍。若已知破片形狀為凸包或類凸包破片時(shí)，可優(yōu)先選用凸包Delaunay剖分算法進(jìn)行計(jì)算；而對(duì)于計(jì)算精度優(yōu)先級(jí)高于時(shí)效的情況，則選用最大邊長約束的凹域剖分算法更為適合。

表2 凸包Delaunay算法與最大邊長約束的凹域剖分算法求解得到時(shí)間 s

4 結(jié)論

基于最大邊長約束的凹域三角剖分算法可對(duì)任意形狀破片的迎風(fēng)面積進(jìn)行求解。相比于凸包Delaunay算法，凹域三角剖分算法計(jì)算精度更高，但求解時(shí)間更長。在已知破片為凸包或類凸包形狀時(shí)，可優(yōu)先采用凸包Delaunay算法進(jìn)行求解。但對(duì)于形狀未知的有限元破片而言，采用凹域三角剖分算法更為適合。