劉星雨，戰(zhàn)仁軍，歐陽的華，李永利

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

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

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

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

1 模型

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

1.1.1 計算模型

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

(1)

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

(2)

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

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

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

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

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

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

(10)

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

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

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

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

EIO=EI-EIN-I

(11)

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

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

1.2 殼體有限元破碎模型

1.2.1 有限元模型

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

圖3 爆震彈結(jié)構(gòu)原理及殼體幾何模型Fig.3 Structure principle and shell geometry model

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

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

(12)

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

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

空氣：采用Null材料本構(gòu)和線性多項式狀態(tài)方程Linear-Polynominal進行建立，即

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

(13)

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

表1 爆震劑材料模型參數(shù)

表 2 空氣材料模型參數(shù)

殼體：采用PLASTIC_KINEMATIC塑性隨動材料本構(gòu)，該模型應(yīng)變率采用了Cowper-Symonds模型。動態(tài)屈服應(yīng)力表達式為：

(14)

(15)

式中，Etan為切線模量，E為彈性模量。其材料模型參數(shù)見表3[2]，ρ為材料密度，V為泊松比，σ0為屈服極限，F(xiàn)S為最大有效塑性應(yīng)變。

表3 殼體材料模型參數(shù)

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

1.3 破片飛散模型

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

(16)

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

(17)

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

2 結(jié)果分析與討論

2.1 沖擊波超壓對比

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

2.2 質(zhì)量分布對比

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

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

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

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

兩類破片的質(zhì)量統(tǒng)計見表4。表4中，半預(yù)制破片質(zhì)量最大值、極差、均值、標準差均比自然破片小。

表4 破片質(zhì)量統(tǒng)計

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

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

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

自然破片仿真結(jié)果與試驗結(jié)果存在差異的原因是塑形隨動本構(gòu)模型是基于唯象理論并結(jié)合單元失效法出的破片失效與斷裂，難以真實模擬出材料微孔洞及應(yīng)力分布不均勻?qū)е碌臍んw裂紋擴展的不確定性，使得仿真自然破片與真實破片的形狀、質(zhì)量分布與試驗結(jié)果存在一定的差異。在后續(xù)研究中，可引入細觀尺度下的弱點隨機分布因子，對現(xiàn)有本構(gòu)進行二次開發(fā)，以使仿真自然破片與真實破片的形狀、質(zhì)量分布更加吻合。

2.3 初始速度對比

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

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

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

表5 不同方向的破片初始速度統(tǒng)計

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

2.4 平均迎風面積等初始參數(shù)

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

表6 MCSPS算法下的自然與半預(yù)制破片初始參數(shù)

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

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

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

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

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

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

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

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

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

表9 自然與半預(yù)制破片致死、致傷半徑

3 結(jié)論

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

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

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