劉建斌，郭德龍，任云燕，徐豫新,2,3,4，李旭東

（1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081；2.高能量密度材料教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081；3.北京理工大學(xué) 北京理工大學(xué)重慶創(chuàng)新中心，重慶 401120；4.北京理工大學(xué) 唐山研究院，河北，唐山 063000；5.中北大學(xué) 機(jī)電工程學(xué)院，山西，太原 030051）

破片式戰(zhàn)斗部主要是在高能炸藥爆炸作用下形 成大量高速破片，利用破片的高速撞擊、引燃和引爆作用毀傷目標(biāo)[1].在實(shí)戰(zhàn)中，破片通常在空氣中飛行一段距離后與目標(biāo)交會，由于空氣阻力的作用，破片飛行速度會逐漸減小.因此，在評估破片式戰(zhàn)斗部實(shí)戰(zhàn)威力時(shí)，通常會考慮破片的速度衰減因素[2?7]，GJB2425-95 中給出了破片阻力系數(shù)的測試要求[8].球形破片因具有氣動性能優(yōu)越、空氣阻力系數(shù)小等優(yōu)點(diǎn)，被廣泛運(yùn)用于預(yù)制破片戰(zhàn)斗部，且在空氣中不發(fā)生翻轉(zhuǎn)，易于在方法層面進(jìn)行空氣阻力系數(shù)特性研究，在此采用彈道槍對球形破片進(jìn)行空氣阻力系數(shù)特性研究，可為后續(xù)其他形狀破片的阻力系數(shù)特性研究提供技術(shù)支撐.

目前國內(nèi)外對球形破片的阻力系數(shù)特性早已開展了研究.例如，CHARTES 等[9]和王儒策等[10](1945)通過試驗(yàn)研究了?9/16 in(?14.3 mm)球形破片飛行馬赫數(shù)在0.29～3.96 范圍內(nèi)空氣阻力系數(shù)變化規(guī)律，分析了超音速下激波離體位置與馬赫數(shù)關(guān)系，獲得不同馬赫數(shù)區(qū)間下的空氣阻力系數(shù).HENDERSON[11](1976)研究了球體的阻力系數(shù)與馬赫數(shù)、雷諾數(shù)、球的溫度和大氣溫度之間的函數(shù)關(guān)系，并給出了函數(shù)模型.譚多望等[12?15](2002、2007、2008、2011)通過二級氫氣炮發(fā)射和靜爆試驗(yàn)方法測量了不同材料、不同尺寸的鋯珠、鋼珠和鎢珠的速度隨飛行距離的衰減情況，結(jié)果表明：對于理想球形鎢合金破片，在同一初始速度條件下，衰減系數(shù)為常數(shù)，空氣阻力系數(shù)與初始速度有關(guān)，兩者成線性關(guān)系.黃德雨等[16](2011)對低密度非金屬球形破片在長距離飛行下的速度衰減規(guī)律進(jìn)行了研究，結(jié)果表明破片在低速長距離飛行時(shí)的空氣阻力系數(shù)與馬赫數(shù)相關(guān).ALLEN[17](2018)在CHARTES 和THOMAS 數(shù)據(jù)基礎(chǔ)上基于量綱分析建立了球形破片飛行馬赫數(shù)在0.2～2.0 區(qū)間下存余馬赫數(shù)與飛行距離之間的函數(shù)關(guān)系式.此外，龐忠根[18](2009)、張華麗[19](2015)、李麗萍[20](2016)通過針對特定尺寸鎢球速度衰減規(guī)律進(jìn)行了試驗(yàn)研究并擬合了存余馬赫數(shù)與飛行距離函數(shù)關(guān)系式，適用馬赫數(shù)區(qū)間為3.0～6.0.綜上所述可見，球形破片空氣阻力系數(shù)研究仍十分活躍，目前對于球形破片空氣阻力系數(shù)研究通常通過試驗(yàn)獲得，一般認(rèn)為破片速度與阻力系數(shù)呈現(xiàn)先增后減趨勢，峰值點(diǎn)出現(xiàn)在馬赫數(shù)1.4 左右.在亞音速至高超音速速度區(qū)間內(nèi)，空氣阻力系數(shù)、空氣阻力隨馬赫數(shù)變化規(guī)律尚無系統(tǒng)研究，飛行距離、初始馬赫數(shù)和存余馬赫數(shù)三者間函數(shù)模型在工程上常用轉(zhuǎn)換現(xiàn)有模型實(shí)現(xiàn)，阻力系數(shù)近似于常數(shù)，難于應(yīng)用于初速2 000 m/s 以內(nèi)的速度衰減效應(yīng)計(jì)算以及爆炸安全領(lǐng)域的安全距離計(jì)算.

通過試驗(yàn)與數(shù)值模擬相結(jié)合的方式獲得球形破片在飛行速度為6.0 馬赫數(shù)以內(nèi)的空氣阻力系數(shù)和空氣阻力隨馬赫數(shù)的變化規(guī)律，并在此基礎(chǔ)上建立了球形破片飛行距離計(jì)算模型.相關(guān)研究成果可為球形預(yù)制破片戰(zhàn)斗部威力場模型構(gòu)建，毀傷效能評估以及爆炸安全領(lǐng)域的安全距離計(jì)算提供支撐，也可為其他類型破片空氣阻力系數(shù)特性的研究提供方法借鑒.

1 破片速度測試試驗(yàn)

1.1 假設(shè)條件

首先，建立以下假設(shè)條件：

①試驗(yàn)測得的平均速度表示該區(qū)間的中點(diǎn)位置速度.

②通過試驗(yàn)測得的速度在計(jì)算空氣阻力系數(shù)過程中，空氣阻力系數(shù)Cd為常數(shù).

③空氣阻力系數(shù)Cd對應(yīng)的速度采用計(jì)算Cd值時(shí)的速度平均值表示.

1.2 破片速度測試方案及試驗(yàn)裝置

試驗(yàn)所用的破片速度加載裝置選用12.7 mm 口徑滑膛彈道槍，測時(shí)儀為南京理工大學(xué)LNG202G-2型，測時(shí)精度為1 μs，量程為1 μs～99 999 s.圖1 為試驗(yàn)裝置示意圖，包括彈道槍、破片、測速靶和測時(shí)儀，圖2 為試驗(yàn)現(xiàn)場圖，每發(fā)試驗(yàn)最多可測量4 組速度.根據(jù)測試儀器使用說明，斷通靶相對通斷靶具有更高的測試精度，為提高測試精度[18]，測速靶選用鋁箔斷通靶（簡稱鋁箔靶），試驗(yàn)過程中，破片依次穿過5 個(gè)鋁箔靶，靶間距分別為l1～l4；t0時(shí)刻，破片穿過Ⅰ靶，測時(shí)儀開始計(jì)時(shí)，接著依次穿過Ⅱ、Ⅲ、Ⅳ、Ⅴ靶，測時(shí)儀分別記錄穿過的時(shí)刻t1～t4，通過li/(ti?ti-1)計(jì)算得到破片穿過相鄰的2 個(gè)測速靶區(qū)間的平均速度v1～v4.靶間距用卷尺測量，卷尺最小刻度為1 mm.

圖1 破片速度測試裝置示意圖Fig.1 Schematic of test equipment

圖2 破片速度測試實(shí)驗(yàn)現(xiàn)場圖Fig.2 Layout of test site

破片速度衰減系數(shù)k由式（1）計(jì)算獲得：

式中：Si為Ⅰ靶至第i個(gè)測速區(qū)間中心位置距離；n為測速區(qū)間個(gè)數(shù)；vi為球形破片經(jīng)過第i個(gè)測速區(qū)間平均速度.空氣阻力系數(shù)Cd可通過式Cd= 2km/(ρa(bǔ)S)計(jì)算獲得，其中k為速衰減系數(shù)；ρa(bǔ)為當(dāng)?shù)乜諝饷芏?；m為破片質(zhì)量；S為破片迎風(fēng)面積.

1.3 試驗(yàn)結(jié)果

試驗(yàn)測量的球形破片包括：?6 mm 鋼球、?6 mm鎢球、?11 mm 鋼球和?11 mm 鎢球，球形破片和彈托實(shí)物如圖3 所示.空氣阻力系數(shù)計(jì)算過程中，不同尺寸、材質(zhì)的球形破片參數(shù)列于表1[21].試驗(yàn)共獲得了23 組?6 mm 和?11 mm球形破片空氣阻力系數(shù)列于表2，每組數(shù)據(jù)均由單枚破片測試獲得.

表1 試驗(yàn)用球形破片參數(shù)Tab.1 Diameter and mass of fragments for tests

表2 球形破片空氣阻力系數(shù)試驗(yàn)結(jié)果Tab.2 Test results of air resistance coefficient for fragments

圖3 球形破片和彈托實(shí)物圖Fig.3 Spherical fragments and sabots

將表2 中數(shù)據(jù)繪制成散點(diǎn)圖如圖4 所示.

圖4 球形破片空氣阻力系數(shù)試驗(yàn)結(jié)果Fig.4 Test result of air drag coefficient of spherical fragment

結(jié)合表2 和圖4 可知，空氣阻力系數(shù)與速度大小相關(guān)；在速度近似相同條件下，4 種不同尺寸或材質(zhì)球形破片空氣阻力系數(shù)試驗(yàn)結(jié)果未顯著表現(xiàn)出與尺寸或材質(zhì)相關(guān).破片飛行馬赫數(shù)從0.70 增大至1.05 時(shí)，空氣阻力系數(shù)從0.59 增大至1.09，增大了約85%，可以看出球形破片飛行馬赫數(shù)在接近1 時(shí)空氣阻力系數(shù)急劇增大；馬赫數(shù)大于3 時(shí)，空氣阻力系數(shù)顯現(xiàn)出隨馬赫數(shù)增大而減小的趨勢.

2 數(shù)值模擬

2.1 模擬區(qū)域和網(wǎng)格劃分

選用適用范圍廣、物理模型豐富、對網(wǎng)格要求不太高并且計(jì)算精度尚可的Fluent 進(jìn)行計(jì)算，數(shù)值模擬區(qū)域形狀為圓柱形，圓柱的高度為60D(D為球形破片直徑)，直徑為30D，計(jì)算域尺寸設(shè)定如圖5 所示.入流邊界距離球心為20D，出流邊界距球心位置為40D，球心距圓柱形計(jì)算域側(cè)邊界各為15D.網(wǎng)格由ICEM CFD 生成，近壁面邊界層包含層流邊界層和湍流邊界層，通過計(jì)算不同馬赫數(shù)下壁面函數(shù)yplus 值選取合適的第一層網(wǎng)格的高度，第二層網(wǎng)格的高度為第一層網(wǎng)格高度的1.2 倍，依次類推；通過數(shù)值模擬試算，第一層網(wǎng)格的高度與馬赫數(shù)間的關(guān)系列于表3，網(wǎng)格類型為四面體非結(jié)構(gòu)化網(wǎng)格.

表3 第一層網(wǎng)格高度Tab.3 Inner-layer grid size under different Mach numbers

圖5 模擬區(qū)域尺寸與邊界設(shè)定示意圖Fig.5 Schematic of calculation area size and boundary setting

2.2 求解器設(shè)置

求解器為Fluent 15.0，設(shè)置流場域四周為壓力遠(yuǎn)場邊界條件(pressure far field)，用于模擬無窮遠(yuǎn)處的自由流條件；球形破片表面設(shè)置為壁面(wall)，用于限定Fluid 和Solid 區(qū)域，如圖5 所示；流體密度模型采用理想空氣(ideal-gas)，黏度模型采用經(jīng)典的Sutherland 模型；對于湍流模型的選取，通過數(shù)值模擬試算對比了分別采用一方程模型Spalart-Allmaras(S-A)模型和二方程模型k-epsilon 模型和k-omega 模型的計(jì)算結(jié)果，結(jié)果表明在其他參數(shù)相同的條件下，S-A模型、k-epsilon 模型和k-omega 模型的計(jì)算結(jié)果相差不大，考慮到數(shù)值計(jì)算的效率，因此湍流模型選用S-A 模型.

2.3 數(shù)值模擬結(jié)果

分別對?6 mm 和?11 mm 球形破片建立數(shù)值模擬模型，獲得2 種直徑球形破片飛行馬赫數(shù)小于6.0時(shí)的空氣阻力系數(shù)數(shù)值模擬結(jié)果如圖6 所示.選取部分馬赫數(shù)下?6 mm 和?11 mm 球形破片空氣阻力系數(shù)數(shù)值模擬結(jié)果列于表4.由圖6 和表4 可見：①馬赫數(shù)小于1.2 時(shí)，空氣阻力系數(shù)逐漸增大；馬赫數(shù)大于1.2 時(shí)，空氣阻力系數(shù)逐漸減小；② 馬赫數(shù)為0.7時(shí)，2 種直徑球形破片阻力系數(shù)差值最大為?3.70%，可認(rèn)為球形破片空氣阻力系數(shù)與直徑無關(guān).

表4 2 種直徑球形破片空氣阻力系數(shù)數(shù)值模擬結(jié)果Tab.4 Typical numerical simulation results of air resistance coefficient

圖6 2 種直徑球形破片空氣阻力系數(shù)數(shù)值模擬結(jié)果Fig.6 Comparison of air resistance coefficient for simulation results

表5 所示為空氣阻力系數(shù)數(shù)值模擬與試驗(yàn)結(jié)果.

表5 試驗(yàn)與數(shù)值模擬結(jié)果對比Tab.5 Comparison of air resistance coefficient between test results and numerical simulation results

馬赫數(shù)為1.00 時(shí)，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果相對誤差最大，為11.11%；馬赫數(shù)為0.70 時(shí)，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果相對誤差最小，為?8.47%，誤差范圍未超過20%[21]，但也說明在跨音速段，計(jì)算誤差整體較大.

選取?6 mm 球形破片空氣阻力和流場靜壓云圖數(shù)值模擬結(jié)果分別列于表6 和圖7.

表6 ?6 mm 球形破片空氣阻力數(shù)值模擬結(jié)果Tab.6 Numerical simulation results of air resistance for ?6 mm spherical fragment

從圖7 中可以看出，當(dāng)飛行馬赫數(shù)為0.10 時(shí)，低壓區(qū)（大氣壓為93 630 Pa）分布在球形破片的中部，隨著馬赫數(shù)增大，低壓區(qū)向尾部聚集.結(jié)合表6 可見，空氣阻力隨馬赫數(shù)的增大而增大，壓差阻力占比總空氣阻力達(dá)90%以上，摩擦阻力占總空氣阻力比例10%以下.將?6 mm 球形破片空氣阻力數(shù)值模擬結(jié)果繪制成曲線如圖8 所示.由圖8 可見，空氣阻力與馬赫數(shù)近似呈二次方關(guān)系[19]；在接近音速時(shí)，即馬赫數(shù)在0.8 左右時(shí)空氣阻力急劇增加.

圖8 ?6 mm 球形破片空氣阻力數(shù)值模擬結(jié)果Fig.8 Numerical simulation results of air resistance for ?6 mm spherical fragment

3 結(jié)果與分析

3.1 空氣阻力系數(shù)試驗(yàn)與數(shù)值模擬結(jié)果與分析

將試驗(yàn)與數(shù)值模擬獲得的?6 mm 和?11 mm 形破片的空氣阻力系數(shù)繪制成散點(diǎn)圖和點(diǎn)線圖，并與文獻(xiàn)[10]中給出的球形破片空氣阻力系數(shù)進(jìn)行對比，如圖9 所示.由圖9 可見，破片飛行馬赫數(shù)小于等于3.0 左右時(shí)，文獻(xiàn)[10]中給出的結(jié)果顯著偏小，大于3.0 左右時(shí)，數(shù)值模擬結(jié)果和文獻(xiàn)[10]中數(shù)據(jù)較為一致.

圖9 球形破片空氣阻力系數(shù)對比Fig.9 Comparison of air resistance coefficient for test, reference and simulation results

對?6 mm 球形破片空氣阻力系數(shù)試驗(yàn)和數(shù)值模擬結(jié)果進(jìn)行擬合，獲得空氣阻力系數(shù)與速度之間函數(shù)關(guān)系式，擬合結(jié)果如圖10 所示.

圖10 ?6 mm 球形破片空氣阻力系數(shù)擬合結(jié)果Fig.10 Fitting results of air resistance coefficient for ?6 mm spherical fragment

空氣阻力系數(shù)與馬赫數(shù)間函數(shù)關(guān)系式為

式中Ma為馬赫數(shù)，

式中：v為破片飛行時(shí)的速度；vs為破片飛行時(shí)當(dāng)?shù)匾羲?

3.2 破片飛行距離計(jì)算模型

破片在空氣中飛行時(shí)，空氣阻力大小與速度的平方成正比[22]，具體形式如式(4)所示，該式通常也被稱之為阻力二次方定律.

式中：F為空氣阻力；Cd為空氣阻力系數(shù)，為量綱一量；ρa(bǔ)為當(dāng)?shù)乜諝饷芏?；S為破片迎風(fēng)面積；v為破片飛行速度.在忽略重力影響的情況下，式(4)可變換為

式中m為破片的質(zhì)量.在式(5)的基礎(chǔ)上，結(jié)合式(2)建立破片飛行距離x與初速v0和存速vx間函數(shù)模型.在過程中，對相關(guān)變量進(jìn)行量綱一化，建立量綱一關(guān)系式[17].

首先，定義量綱一數(shù)z：

式中：D為球形破片的直徑；ρp為破片密度；ρa(bǔ)為空氣密度；x為破片飛行距離；kz為比例因子，與距離x單位相同.

對式(5)進(jìn)行變換：

結(jié)合式(3)：

球形破片質(zhì)量m跟密度ρp、直徑D之間存在如下關(guān)系式：

結(jié)合式(9)和(10)，式(8)可變換成：

當(dāng)Cd=A時(shí)，代入式(11)：

當(dāng)Cd=A+BMa+CMa2時(shí)，代入式(11)：

式中：Ma0=v0/vs，Max=vx/vs，Max是指破片初始馬赫數(shù)為Ma0時(shí)，飛行距離x后的存余馬赫數(shù)；z0= 0；A、B、C為系數(shù).

通過對式(12)和式(13)進(jìn)行積分，結(jié)合式(6)可得空氣阻力系數(shù)為常數(shù)和與馬赫數(shù)成二次函數(shù)時(shí)的飛行距離x與初始馬赫數(shù)Ma0、存余馬赫數(shù)Max之間關(guān)系式：

選取試驗(yàn)數(shù)據(jù)和文獻(xiàn)[12]中給出的兩種材質(zhì)的球形破片長距離飛行時(shí)的試驗(yàn)數(shù)據(jù)，對比建立的破片飛行距離計(jì)算模型計(jì)算結(jié)果分別列于表7 和表8.從表7 和表8 中可以看出，通過建立的飛行距離計(jì)算模型計(jì)算的結(jié)果與試驗(yàn)結(jié)果及文獻(xiàn)[9]中給出的試驗(yàn)結(jié)果最大誤差在3%左右.

表7 本文試驗(yàn)結(jié)果與模型計(jì)算結(jié)果對比Tab.7 Comparisons of flight distance between test data and calculation results of this model

表8 文獻(xiàn)[12]試驗(yàn)結(jié)果結(jié)果與模型計(jì)算結(jié)果對比Tab.8 Comparisons of flight distance between references [12] test data and calculation results of this model

在戰(zhàn)斗部動爆威力計(jì)算時(shí)，球形破片空氣阻力系數(shù)通常取0.97[2?4]；此外，王儒策等[10]也給出了不同速度區(qū)間球形破片空氣阻力系數(shù).對比球形破片初始馬赫數(shù)為6.0，存余馬赫數(shù)分別為3.0、1.1 和0.5時(shí)，采用0.97 及文獻(xiàn)[10]中空氣阻力系數(shù)和上述建立的飛行距離計(jì)算模型計(jì)算的飛行距離結(jié)果列于表9，計(jì)算過程中所需的破片參數(shù)及大氣參數(shù)取表8 中鎢球?qū)?yīng)的參數(shù).從表9 中可知，在初始馬赫數(shù)取6.0條件下，當(dāng)存余馬赫數(shù)取3.0 時(shí)，空氣阻力系數(shù)取0.97計(jì)算的飛行距離相比上述模型計(jì)算的結(jié)果偏大0.11%，近似相等，采用文獻(xiàn)[10]中的空氣阻力系數(shù)計(jì)算的飛行距離相比上述模型計(jì)算的結(jié)果偏大4.56%；當(dāng)存余馬赫數(shù)取1.1 時(shí)，空氣阻力系數(shù)取0.97 計(jì)算的飛行距離相比上述模型計(jì)算的結(jié)果偏大6.85%，采用文獻(xiàn)[10]中的空氣阻力系數(shù)計(jì)算的飛行距離相比本文模型計(jì)算的結(jié)果偏大8.30%；當(dāng)存余馬赫數(shù)取0.5 時(shí)，空氣阻力系數(shù)取0.97 計(jì)算的飛行距離相比上述模型計(jì)算的結(jié)果偏小10.58%，采用文獻(xiàn)[10]中的空氣阻力系數(shù)計(jì)算的飛行距離相比上述模型計(jì)算的結(jié)果大9.00%.

表9 破片飛行距離計(jì)算結(jié)果對比Tab.9 Comparisons of flight distance results for different models

4 結(jié) 論

通過對?6 mm 和?11 mm 2 種直徑球形破片進(jìn)行空氣中運(yùn)動阻力系數(shù)特性試驗(yàn)和數(shù)值模擬研究，獲得了球形破片空氣中飛行馬赫數(shù)在6.0 內(nèi)空氣阻力及系數(shù)隨速度的變化規(guī)律，結(jié)論如下：

①相同飛行速度下，不同直徑球形破片空氣阻力系數(shù)相同；球形破片的空氣阻力系數(shù)隨馬赫數(shù)的增加先增大后減小，并非一恒定值，馬赫數(shù)為1.2 時(shí)最大.

②球形破片空氣阻力在馬赫數(shù)為0.8 左右時(shí)急劇增加；其中，壓差阻力占比90%以上，摩擦阻力占比10%以下.

③建立了球形破片空氣阻力系數(shù)隨馬赫數(shù)變化的函數(shù)和球形破片飛行距離計(jì)算模型，具有較高的計(jì)算精度.