蒙特卡洛剖分投影法的爆震彈破片平均比動能評估*

劉星雨，戰仁軍，歐陽的華，李永利

(1. 武警工程大學 裝備管理與保障學院， 陜西 西安 710086； 2. 上海交通大學 海洋智能裝備與系統教育部重點實驗室， 上海 200240；3. 西安建筑科技大學 資源工程學院， 陜西 西安 710055)

爆震彈是利用含能材料爆炸時產生的強烈聲響和眩目閃光，使有生目標暫時致聾致盲的非致命性彈種[1]。此類彈種采用手投、槍發或發射器方式進行發射，被廣泛用于部隊反恐、處突、海上維權執法等任務中。其封裝殼體多為塑料或鋁制材料，在爆炸過程中形成的較大破片易對人造成不可逆轉的傷害。為提升其安全性，國內外學者目前主要采用的是改善殼體材料或外置高強度開孔金屬殼的方法[2]。而對于采用殼體預刻槽方法控制破片的質量分布、飛散速度、比動能，進而減小破片場的致傷與致死半徑的研究則較少。特別是基于上述指標下破片“成型、飛散、著靶”全時域的理論模型與仿真方法，還少有研究。為此，建立破片全時域的仿真模型，自動化地求解破片質量分布、速度、平均比動能、致傷半徑、致死半徑等安全性指標，是高效地指導武器使用、優化武器性能、節約研發成本的客觀需求。

目前，國內外學者對殺傷性彈藥破片的研究相對較多，特別是殺傷性預制破片戰斗部威力場的數值仿真研究，提出了很多成熟的仿真模型與計算方法，如文獻[3-10]基于LS-DYNA、AUTODYN平臺對各類戰斗部破片場進行了建模、仿真評估、效能優化、快速計算等。但對非致命彈藥破片的安全特性研究，特別是全時域下的破片比動能研究還相對較少。上述研究中依托LS-DYNA、AUTODYN等通用有限元平臺，通常只能獲得爆炸初始時間段內的破片場飛散特性及沖擊波壓力，難以得到全時域全破片的平均比動能。而在建立全時域的飛散模型及平均比動能計算過程中，其重要參數平均迎風面積目前主要采用均勻取向理論的試驗方法獲得，而基于有限元結果的任意形狀破片平均迎風面積求解還未見相關報道。

故本文根據爆震彈破片的比動能定義、形成機理及飛散原理，建立了平均比動能計算模型，并結合有限元殼體破碎模型及破片飛散模型，對同一裝藥下的自然破片和半預制破片全時域平均比動能及安全半徑進行了評估。其中，平均比動能計算模型通過蒙特卡洛剖分投影(Monte Carlo Subdivision Projection Simulation，MCSPS)方法建立，可對任意形狀有限元破片的平均迎風面積及比動能進行求解；破碎模型采用了相容流固耦合方法進行數值模擬；飛散模型是經有限元質量損耗修正后的質點外彈道模型。通過模型求解，可得到距爆心任意位置上的破片垂直靶分布、平均比動能閾值，進而對自然破片與半預制破片的致傷、致死半徑進行了評估。

1 模型

1.1 基于MCSPS法的平均比動能計算模型

1.1.1 計算模型

根據破片的平均比動能定義，其表達式[11]為：

(1)

MCSPS方法通過對有限元破片的坐標點集進行N次平移、隨機旋轉、投影、三角剖分、凹域邊界搜索排序、叉乘法求面積并取平均，可得到N次投影面積Aij的平均值。當N值趨于無窮大時，所得均值結果可近似為該破片的平均迎風面積。即

(2)

各計算步驟如式(3)～(9)所示。

(3)

(4)

(5)

(6)

(7)

(8)

(9)

其中：f1，f2，f3分別為對坐標向量的平移、旋轉和投影運算；f4為平面坐標向量的Delaunay凸包三角剖分[12]，P為三角剖分凸包點集矩陣，Clist為凸包三角形矩陣；f5為帶最大邊長約束的三角剖分凹域邊界搜索與排序，Pb為排序后的邊界頂點向量；f6為叉乘法求多邊形面積，Aij為第i枚破片第j次MCSPS后的迎風面積；f7為進行N次MCSPS后所得迎風面積取平均，N為MCSPS總次數。

基于MCSPS的平均迎風面積仿真算法流程如圖1所示。其中，Pe為破片單元數據庫，Pf為破片匯總數據庫，i2、i分為Pe庫行搜索標記的首端和末端，q為破片序號，mtkl為蒙特卡洛次數變量，N為蒙特卡洛總次數，設置為50 000次，area_mtkl為蒙特卡洛投影面積矩陣。

圖1 基于MCSPS的平均迎風面積仿真算法流程Fig.1 Flow chart of simulation algorithm of average windward area based on MCSPS

1.1.2 凹域三角剖分邊界搜索排序

對于投影點集為凸包的情況，可直接根據計算機圖形學原理，采用Delaunay三角剖分進行求解。而對于凹多邊形圍成圖形區域(簡稱凹域)，則無法直接進行求解。提出采用最大邊長約束的原則，將凸包剖分算法所得大于約束值R的三角邊去除，對凹域邊界重新進行搜索，并用叉乘法求解凹多邊形面積，可大大提高平均迎風面積的精度及計算效率。R取值為：

(10)

式中：R為邊長約束值；λscale為放大因子，取8；li為剖分三角形第i條邊的長度；ne為三角邊總數。

帶最大邊長約束的凹域三角剖分邊界搜索排序算法流程如圖2所示。圖2中：T為無重復節點的投影點集坐標矩陣；dm為凸包剖分邊界矩陣；dt_Clist為凸包三角剖分三角形節點編號矩陣；edges為凸包三角剖分升序排列的邊界矩陣；I為edges中不包含重復邊的行號向量；IN為edges的行號向量；IO為邊界邊行號向量；sort_vertices為用于存放邊界頂點序列的向量；polygon為按順(逆)時針排序的邊界邊節點坐標矩陣。

圖2 帶最大邊長約束的三角剖分邊界搜索算法流程圖Fig.2 Flow chart of boundary search-sort algorithm with maximum edge length constrain of triangulation

凹域三角剖分邊界搜索的差集運算表達式為：

EIO=EI-EIN-I

(11)

式中：EIO為凹多邊形區域邊界邊集合，有且僅有一條公共邊；EIN-I為凹多邊形區域的內部邊集合，存在兩個公共邊；EI為不包含重復邊的凹多邊形邊集合。

通過搜索算法得到邊界頂點坐標矩陣polygon，采用叉乘法求解多邊形面積可得到Aij，將結果代入式(2)和式(1)可得到破片的平均迎風面積及平均比動能。

1.2 殼體有限元破碎模型

1.2.1 有限元模型

某型爆震彈其結構原理如圖3(a)所示。全彈長度為110 mm(不含擊發機構)，直徑為38 mm；主裝藥為一種含鋁煙火藥，質量為60 g；殼體為ABS塑料材料。對擊發機構及殼體結合部作適當簡化，無刻槽與內刻槽殼體幾何模型如圖3(b)、圖3(c)所示。其中，1/4內刻槽殼體采用4條周向槽及2條軸向槽進行劃分，槽深2 mm，槽寬2 mm。

圖3 爆震彈結構原理及殼體幾何模型Fig.3 Structure principle and shell geometry model

通過相容流固耦合方法對殼體破碎過程進行1/4對稱建模，如圖4所示，爆震劑及空氣采用共節點Euler單元，連接座、殼體采用共節點的Lagrange單元。在對稱邊界面設置位移約束條件；在空氣外表面上設置無反射邊界。為了驗證爆震劑參數的可靠性，建立6.5 cm×20 cm×306 cm的空氣單元網格，用于計算3 m內的沖擊波壓力，以便與試驗進行對比，網格大小為0.1 cm×0.1 cm×0.2 cm。為了確保Euler單元應力應變輸運的精度，爆震劑與空氣單元網格采用相同大小進行劃分。為了提高殼體失效過程中的應力應變輸運精度，將靠近殼體附近單元網格尺寸采用0.1 cm×0.1 cm×0.1 cm進行劃分。

圖4 爆震彈流固耦合有限元模型Fig.4 Finite element model of fluid-solid coupling

彈體采用立式放置，炸高距地面1 m。爆震劑采用單點起爆，炸點位置距連接座下端面1.5 cm。

1.2.2 材料本構及失效準則

爆震劑：采用HIGH_EXPLOSIVE_BURN高能炸藥本構，其瓊斯-威爾金斯-李(Jones-Wilkins-Lee, JWL)狀態方程為：

(12)

式中：p為壓力；Vg為爆轟產物體積和炸藥初始體積之比；E0為爆轟產物單位體積的初始內能；A，B，R1，R2，ω為特征參數。

密度ρ根據裝藥質量及體積計算得到，爆速D根據文獻[1]得到，JWL方程各參數借鑒含鋁炸藥文獻[13]，并作以適當調整，可使1 m、2 m、3 m處的沖擊波超壓與試驗基本吻合。各參數取值見表1。

空氣：采用Null材料本構和線性多項式狀態方程Linear-Polynominal進行建立，即

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E0

(13)

式中，C0～C6為與氣體性質有關的狀態方程參數，C0、C1、C2、C3為壓強量綱，C4、C5、C6無量綱。地面空氣μ=0，排除大氣相對壓力，由此得出的空氣材料模型參數見表2。

表1 爆震劑材料模型參數

表 2 空氣材料模型參數

殼體：采用PLASTIC_KINEMATIC塑性隨動材料本構，該模型應變率采用了Cowper-Symonds模型。動態屈服應力表達式為：

(14)

(15)

式中，Etan為切線模量，E為彈性模量。其材料模型參數見表3[2]，ρ為材料密度，V為泊松比，σ0為屈服極限，FS為最大有效塑性應變。

表3 殼體材料模型參數

LS-DYNA中add_erosion關鍵字提供了最大主應變失效、最大主應力失效、等效塑性應變失效、最大靜水拉壓失效等失效準則。經過多次試算對比，并結合材料本構中最大有效塑性應變FS取值，發現采用最大主應變失效準則可較好地解決殼體裂紋附近單元出現大面積失效的情況。

1.3 破片飛散模型

當殼體完全破碎后，沖擊波對破片的速度影響逐漸減弱，直至某一時刻速度不再顯著增加。此后，各破片在空氣阻力FR、自身重力G及旋轉力矩M作用下做自由飛行。根據質點外彈道學原理[14]并考慮有限單元失效法中質量、平均迎風面積的損耗因素，引入破碎前殼體質量mr進行修正。可推得經修正后的第i枚破片在二維平面內速度vi、水平夾角θi、水平位移xi、垂直位移yi和比動能edi隨時間t變化的常微分方程組。

(16)

設置迭代終止條件為yin≤0或tin≥Tmax(Tmax為時間閾值，取10 s)。迭代初始值yi0從破碎模型中直接讀取，vi0、θi0、xi0需從破碎模型當前幀及前一幀破片單元的速度、坐標分量數據中進行計算，其表達式為：

(17)

根據空間坐標矩陣變換關系，將各破片二維軌跡矩陣[xi,yi]按各破片初始方位角繞坐標軸y旋轉，可獲得各破片的三維飛散軌跡矩陣[x′i,y′i,z′i]以及垂直靶分布情況。

2 結果分析與討論

2.1 沖擊波超壓對比

為了驗證爆震劑參數的可靠性，進行自由場沖擊波超壓測試試驗，測得了距爆心1 m、2 m、3 m處的沖擊波超壓峰值并與仿真超壓進行了對比，不同距離上的仿真與試驗結果對比如圖5所示。圖5中，無刻槽與內刻槽殼體在1 m、2 m、3 m處的沖擊波超壓峰值仿真結果基本相同；兩類殼體在1 m、2 m處的超壓峰值與試驗結果基本吻合，但在3 m處試驗結果比仿真結果更低。3 m處仿真結果與試驗結果產生誤差的原因主要有兩個方面：一是JWL狀態方程及高能炸藥本構模型主要是基于TNT炸藥而建立，對于含鋁煙火藥的爆炸數值模擬會存在一定誤差；二是LS-DYNA等商用有限元軟件無反射邊界面附近仍會存在沖擊波超壓反射的情況，使得邊界附近單元壓力值偏大。

2.2 質量分布對比

通過破碎模型得到的未刻槽、內刻槽殼體破碎情況如圖6(a)、圖6(b)所示。兩類殼體在40 μs 時刻已完全破碎。根據飛散方向不同，將破片分為頂部破片、底部破片和徑向破片，并進行編號，如圖6(c)、圖6(d)。從圖6可知，內刻槽半預制破片比無刻槽自然破片數量多，徑向破片尺寸更小，質量分布更加集中。

(a) 未刻槽殼體破碎情況(a) Breaking situation of non-grooved shell fragments

通過同一破片搜索排序算法，可得到以破片為單位的數量、質量、單元數、體積等。其質量頻數分布情況如圖7所示。圖7中，兩類破片質量分布呈現指數衰減的規律，大多數破片主要集中于0～1 g之間；除頂部擊發座破片超過3 g外，半預制破片1～3 g的破片數比自然破片更少；破片均值和標準差更小，質量分布更加集中。

圖7 仿真破片質量分布折線圖Fig.7 Line graph of mass-frequency distribution of simulation fragments

兩類破片的質量統計見表4。表4中，半預制破片質量最大值、極差、均值、標準差均比自然破片小。

表4 破片質量統計

綜上可知，兩類殼體結構的不同導致了應力集中點位置的不同，使得各破片單元失效的擴展方向產生差異，進而影響了破片的形狀及質量。僅從質量因素上看，半預制破片可在一定程度上降低破片比動能，但具體情況還需和破片速度及平均迎風面積一起進行分析。

另外，通過破片回收試驗，得到了同一殼體參數下的自然破片與半預制破片的質量分布情況，如圖8所示。從質量分布上看，兩類殼體的試驗破片質量分布情況與仿真破片結果較為相似，絕大多數破片質量在0～1 g區間內，且兩類破片質量呈指數衰減規律。但從破片形狀上看，內刻槽仿真破片與試驗破片形狀較為吻合，但無刻槽仿真破片No.7、No.9、No.11、No.13與試驗破片形狀仍存在一定差異。

圖8 試驗破片質量分布折線圖Fig.8 Line graph of mass-frequency distribution of test fragments

自然破片仿真結果與試驗結果存在差異的原因是塑形隨動本構模型是基于唯象理論并結合單元失效法出的破片失效與斷裂，難以真實模擬出材料微孔洞及應力分布不均勻導致的殼體裂紋擴展的不確定性，使得仿真自然破片與真實破片的形狀、質量分布與試驗結果存在一定的差異。在后續研究中，可引入細觀尺度下的弱點隨機分布因子，對現有本構進行二次開發，以使仿真自然破片與真實破片的形狀、質量分布更加吻合。

2.3 初始速度對比

根據破碎模型結果，通過LS-PREPOST后處理平臺對頂部、底部、徑向的各破片建立跟蹤單元，其初始速度的時間歷程如圖9所示。圖9中：兩類破片在50 μs之后速度逐漸趨于穩定，不再顯著增加；底部方向破片單元受旋轉力矩影響較大，速度曲線呈現周期振蕩，需采用中值濾波算法對其進行平滑處理后再進行數據讀取。

(a) 自然破片初始速度(a) Initial velocity of natural fragments

不同方向上的破片初始速度統計見表5。從表5可知，兩類殼體底部破片速度均大于徑向破片速度，而徑向破片速度又大于頂部破片速度。造成不同方向上速度差異的原因主要有兩個方面：一是頂部破片的卸載效應使其方向上的應力波顯著降低；二是底部端面距起爆點較遠，應力波輸運至此處時會在端面附近產生疊加，使應力波顯著增大。

表5 不同方向的破片初始速度統計

其中，頂部方向半預制破片與自然破片速度差異不大，而徑向、底部方向半預制破片的速度均值、標準差均比自然破片更低。僅從初始速度上看，半預制破片速度閾值更低且更加集中，可在一定程度上降低破片比動能。

2.4 平均迎風面積等初始參數

將破碎模型有限元數據導出，通過MCSPS平均迎風面積仿真程序計算，可得兩類破片的數量、質心坐標、質心速度、單元數、平均迎風面積等初始參數，其中各破片迎風面積數據見表6。表6中，自然破片4、6、8、10、14、15單元數不足3個，節點投影后會出現無法構建三角剖分的超級三角形情況。對于不足3個節點的破片飛屑采用了近似計算的方法。將表6中各破片初始參數代入飛散模型中進行求解，可得到各破片平均比動能隨距離的衰減情況。

表6 MCSPS算法下的自然與半預制破片初始參數

為驗證MCSPS算法求迎風面積的可靠性，選擇長方體、圓柱體預制破片進行求解并與理論公式[15]結果進行對比。設長方體尺寸為1 cm×2 cm×4 cm，網格尺寸為2 mm×2 mm×2 mm；圓柱體為半徑r=1 cm，高度h=4 cm，網格尺寸為2 mm×2 mm×1 mm。當循環次數N=16 384(214)時，MCSPS算法所得長方體破片、圓柱體破片迎風面積結果與理論公式對比見表7。從表7中可知， MCSPS算法結果與理論公式結果基本吻合，其中長方體破片誤差僅為0.12%，而圓柱體破片誤差為5.29%。

表7 MCSPS算法與理論公式平均迎風面積誤差對比

MCSPS算法與理論公式結果存在誤差的原因主要有兩點：一是偽隨機數誤差，即由軟件產生的均勻分布樣本是基于線性同余算法建立的，與真實均勻分布樣本不可避免地存在一定誤差，但此類誤差在16 384次循環下對結果影響較小，如長方體破片的誤差結果；二是剖分精度誤差，即有限元網格是一種對真實物體的剖分重構，網格尺寸的大小會直接影響物體幾何特征的精度，特別是對于圓柱體這類具有曲線、曲面的幾何體而言，網格尺寸過大會導致圓柱體被分割成正多棱柱體，進而所得平均迎風面積將比真實的理論值要小，如圓柱體的誤差結果，達到了5.29%。

2.5 平均比動能及安全半徑對比

根據表6數據及飛散模型，得到了與爆心不同距離的破片垂直靶分布和平均比動能閾值，其中兩類破片1.5 m、5 m處垂直靶分布如圖10所示。從分布上看，兩類破片飛散角相差不大，分布均呈現不完全對稱，其中半預制破片比自然破片橫向分布更稀疏、縱向分布更密集。從比動能數值上看，在1.5 m處自然破片的比動能閾值比半預制破片閾值更大，其徑向破片6、7、8、10、13、14均超過半預制破片的比動能閾值；而在5 m處其比動能閾值反而遠低于半預制破片的比動能閾值。造成此現象的原因主要有兩點：一是自然破片7、13的尺寸較大且形狀更加不規則，單位質量的平均迎風面積相比于半預制破片更大，飛散過程中速度衰減更快，進而使得自然破片比動能衰減得更快；二是由于自然破片6、8、10、14初始射角為負，在不足5 m時就已著地。

(a) 1.5 m處的自然破片(a) Natural fragments at 1.5 m

兩類破片1 m至5 m內的比動能最大值、最小值、均值、標準差等結果見表8。在小于等于2.5 m范圍內，自然破片比動能最大值、均值均高于半預制破片；超過2.5 m后，半預制破片比動能最大值反而高于自然破片。

表8 不同距離全破片的平均比動能評估結果

根據彈丸、破片、小箭殺傷人員的最小比動能標準160 J/cm2以及擦傷皮膚的最小比動能標準9.8 J/cm2[11]，對表7數據進行插值，得到兩類破片致死半徑和致傷半徑，見表9。從表9中可知，相比于自然破片，半預制破片的致死半徑有顯著減小，但致傷半徑并未有顯著改善。根據非致命防暴彈的安全性要求，殼體破片比動能不宜高于擦傷皮膚標準，即自然破片、半預制破片安全半徑分別為4.7 m、4.6 m。

表9 自然與半預制破片致死、致傷半徑

3 結論

1) 基于MCSPS方法建立的平均比動能計算模型，可求解任意形狀的有限元破片平均迎風面積和平均比動能。

2)全時域下的破片比動能評估相對于初始時刻的比動能評估更全面。

3)相比于自然破片，采用半預制破片方案可顯著減小爆震彈破片的致死半徑，但并不能對其安全半徑有顯著改善。要進一步降低其致傷半徑，還應考慮起爆點位置、殼體厚度、材料等因素。