辛大鈞， 薛琨

(北京理工大學 爆炸科學與技術國家重點實驗室， 北京 100081)

0 引言

戰(zhàn)斗部爆炸后，由于爆轟驅動的作用會形成大量高速飛行破片，初速度馬赫數(shù)在5以上，這些破片的速度在空氣中受空氣阻力的作用而逐漸衰減，落地之前將會有70%以上的時間處于亞聲速飛行狀態(tài)。因此，了解和掌握破片在超聲速至亞聲速區(qū)間內的速度衰減規(guī)律才能預測破片的彈道軌跡和不同距離處的殘余動能，對于評價破片式戰(zhàn)斗部毀傷效果以及評估安全距離均有重要意義。

破片在空氣中飛行時主要受重力及空氣阻力的作用，空氣阻力主要與破片的飛行速度、迎風面積、空氣密度以及空氣阻力系數(shù)有關。其中阻力系數(shù)作為影響空氣阻力的關鍵因素具有重要研究意義。之前的研究表明，阻力系數(shù)主要與破片形狀以及飛行速度有關，國內外學者對此開展過一系列研究。

一般認為破片速度與阻力系數(shù)呈現(xiàn)先增后減趨勢，峰值點出現(xiàn)在馬赫數(shù)1.4左右，國內外學者使用彈道槍試驗對不同形狀的典型破片(球形、立方形、圓柱形等)進行了一系列研究，大量的彈道槍研究結果表明，對于不同形狀的破片，阻力系數(shù)與馬赫數(shù)的關系僅僅是趨勢類似，不同馬赫數(shù)下阻力系數(shù)的數(shù)值最高點出現(xiàn)的位置均明顯不同，因此研究者一般認為阻力系數(shù)還受破片形狀的影響。

關于形狀對阻力系數(shù)的影響，Dunn等在對回收的炮彈破片進行阻力系數(shù)測試后，發(fā)現(xiàn)非規(guī)則破片的阻力系數(shù)顯著大于球形、立方形等規(guī)則破片，首次提及了破片形狀對阻力系數(shù)的影響。破片形狀一般用破片和球體的近似程度來表征。Dehn使用破片的最大迎風面積除以破片的平均迎風面積，即=/，或者破片體積的2/3次方除以破片平均迎風面積/作為描述破片形狀的參數(shù)，對于規(guī)則形狀破片，等于柯西迎風面積即1/4表面積，對于不規(guī)則破片，可通過正二十面體測量儀測出32個方向上的迎風面積取平均。McCleskey使用垂直風洞測試了大量翻滾非規(guī)則破片在馬赫數(shù)0.1下的阻力系數(shù)，首次比較系統(tǒng)地建立起了馬赫數(shù)0.1下阻力系數(shù)與Dehn定義的形狀參數(shù)之間的關系，并認為所有破片的阻力系數(shù)隨馬赫數(shù)的變化規(guī)律與球形破片一致，通過這種假設可以將不規(guī)則破片馬赫數(shù)0.1下的阻力系數(shù)試驗結果外推至跨音速及超音速區(qū)間。Miller進一步使用風洞及彈道槍對McClescky報告中的破片進行了阻力系數(shù)測試，得到跨聲速及超聲速段的阻力系數(shù)結果，與McClescky的外推結果相比明顯偏大。因此，簡單地對球形破片阻力系數(shù)規(guī)律進行外推，無法精確地獲得非規(guī)則破片跨聲速及超聲速速度下的阻力系數(shù)。Moxnes等重新分析了McClescky的試驗結果，認為破片阻力系數(shù)與的相關性比更好，應該用作為描述破片形狀的參數(shù)。

由于非球形破片在空氣中會發(fā)生不規(guī)則翻滾，當前研究隨機翻滾狀態(tài)下破片的平均阻力系數(shù)比較理想的方法是采用彈道槍試驗，但彈道槍試驗成本較高，且無法覆蓋所有的破片形狀，對破片尺寸也有要求，很難開展大規(guī)模試驗。因此眾多研究者借助數(shù)值模擬的方法來研究這一問題，Moxnes等用6自由度計算流體動力學(CFD)方法模擬了無約束破片在靜止流場中的自由飛行，并對阻力系數(shù)在時間上進行了微秒量級的平均。但由于計算資源有限，這種方法很難對破片飛行全過程進行模擬，且需要人為給破片施加初始的翻滾力矩，與實際破片經歷爆轟驅動后的初始運動狀態(tài)有一定區(qū)別。

為了研究大量非球形破片的形貌在超聲速至亞聲速范圍內對破片阻力系數(shù)的影響，本文選擇采用流動流場中固定破片的數(shù)值計算方法。該方法的關鍵是如何對破片不同迎風方向下的阻力系數(shù)進行平均，等效破片隨機翻滾狀態(tài)下的阻力系數(shù)。提出一種基于正二十面體的阻力系數(shù)平均方法，計算了大量不同形狀破片在不同馬赫數(shù)下的阻力系數(shù)，并采用與含義一致但學術上更通用的球形度代替，=π6/(π/)。發(fā)現(xiàn)阻力系數(shù)受形狀參數(shù)和馬赫數(shù)的影響規(guī)律很難使用單一函數(shù)準確描述。因此本文采用近年來發(fā)展迅速的深度神經網絡模型對阻力系數(shù)的影響因素進行擬合，獲得了相對準確實用的阻力系數(shù)預測模型。

1 數(shù)值模擬方法及驗證

1.1 球形破片阻力系數(shù)模擬方法

以直徑7 mm球形破片空氣阻力系數(shù)的數(shù)值模擬為例，介紹數(shù)值計算模型。

數(shù)值模擬采用破片固定，流場采用給定速度運動的方式。將進行模擬的球形破片劃分方形外流場區(qū)域，球形破片(直徑為)的球心與入流邊界距離為40，與出流邊界的距離為40，與各個側邊界的距離為20。網格劃分為四面體非結構化網格，流場中網格的最小尺寸定為1.00 mm，破片壁面附近的網格尺寸為0.10 mm。體網格總數(shù)約60 000個，數(shù)值模擬建模如圖1所示。

圖1 數(shù)值模擬Fig.1 Numerical modeling

求解器選擇Fluent 17.0，當馬赫數(shù)大于0.3時，設置為可壓縮流動。湍流模型選擇S-A湍流模型，主要用以求解一個有關渦黏性的運輸方程，計算量相對較小。對逆壓梯度的邊界層問題和壁面限制的流動問題有較好的計算結果，通常應用于空氣動力學問題中。設置流場域四周為壓力遠場邊界條件，用于模擬無窮遠處的自由流條件；流場材料設置為空氣，狀態(tài)為理想氣體，并用薩蘭德定律修正空氣黏度。球形破片表面設置為無滑移壁面，用于限定流體和破片區(qū)域。

差分格式選擇為壓力- 速度耦合方法；壓力插值采用標準格式；空間離散用Roe-FDS(Flux-Difference Splitting)格式，該離散方式具有很高的間斷分辨率和黏性分辨率；擴散項離散采用中心差分格式；對流項離散采用2階迎風格式。

根據上文的數(shù)值計算模型，對直徑7 mm的球形破片在馬赫數(shù)為0.1～6.0速度區(qū)間的阻力系數(shù)進行模擬。流場靜壓云圖如圖2所示。

圖2 直徑7 mm球形破片壓力云圖Fig.2 Pressure field around φ7 mm spherical fragment

由圖2可知：破片處于亞聲速流場時，在來流前端會形成一個高壓區(qū)，來流末端會形成一個低壓區(qū)，此時壓差阻力和摩擦阻力是空氣阻力的主要來源；隨著飛行馬赫數(shù)的提高，來流前端的壓力逐漸升高，高壓區(qū)會逐漸收窄并向破片兩側彎曲，形成一道弓形激波，此時波阻是阻力的主要來源。

本文對文獻[4]球形破片彈道槍試驗結果進行了復現(xiàn)，破片的飛行速度范圍是馬赫數(shù)為1.5～6.0，將試驗結果以及文獻[4]中給出的彈道槍試驗結果與本文數(shù)值模擬結果進行對比，如圖3所示。

圖3 球形破片數(shù)值模擬結果與文獻[4]試驗以及本文復現(xiàn)試驗結果對比Fig.3 Comparison of numerical results of spherical fragments with experimental data in Ref. [4] and those repoduced in this paper

由圖3可見，模擬結果在破片飛行馬赫數(shù)小于0.5時，相比彈道槍試驗結果略微偏小，在破片處于跨聲速以及超聲速飛行狀態(tài)時，模擬結果與試驗結果匹配得較好?？傮w而言，破片飛行馬赫數(shù)小于0.5時阻力系數(shù)基本保持不變；飛行速度達到跨聲速時阻力系數(shù)迅速升高，在馬赫數(shù)為1.4附近達到最大值，此后飛行速度繼續(xù)增加，阻力系數(shù)略有下降。

球形破片數(shù)值模擬結果與試驗結果數(shù)值對比如表1所示。由表1可見：馬赫數(shù)為0.7時，數(shù)值模擬結果與試驗結果相對誤差最大，為-11.62%；馬赫數(shù)為3.0時，數(shù)值模擬結果與試驗結果相對誤差最小，為0.37%，誤差范圍未超過15%。

1.2 非球形破片阻力系數(shù)計算方法

非球形破片不同姿態(tài)下的迎風面積和阻力系數(shù)與破片姿態(tài)密切相關，如圖4所示。將非球形破片固定在笛卡爾坐標系中，迎風方向的單位矢量為(，)，其中表示與其在平面內投影向量的夾角，表示在平面內投影向量與軸的夾角。

破片處于隨機翻滾狀態(tài)時，不同迎風方向對應著不同的阻力系數(shù)(，)，而各個迎風方向在空間上出現(xiàn)的概率相同。因此理想狀態(tài)下，需要對所有迎風狀態(tài)下的阻力系數(shù)進行平均：

表1 球形破片模擬結果與試驗對比Tab.1 Comparison of numerical results of sphericalfragments with experimental data

圖4 非球形破片迎風方向示意圖Fig.4 Upwind direction of non-spherical fragment

(1)

圖5 立方體三類迎風狀態(tài)示意圖Fig.5 Three types of upwind state of the cube

(2)

上述平均方法可以在一定程度上估計方形破片的阻力系數(shù)，但對于任意形狀破片而言，該方法并不具有普適性。理論上，從正多面體中心指向各個面心或頂點的方向可以較好地實現(xiàn)對空間的平均劃分，而正多面體一共有5種，其中面最多，最接近球體的正多面體為正二十面體?；诖耍绹娪醚b備國際試驗操作規(guī)程(ITOP)4-2-813測量平均迎風面積的方法中采用正二十面體測量儀的方法得到不規(guī)則破片的平均迎風面積。具體做法如下：將破片放在虛擬正二十面體的中心，從正二十面體的中心指向其20個面心(，=1，2，…，20)與12個頂點(，=1，2，…，12)，產生32個方向。沿這32個方向對破片進行投影，可以得到對應的投影面積，即20個中心- 面心方向的投影面積和12個中心- 頂點方向的投影面積。這些投影面積求平均后得到的平均投影面積，可以等效破片的平均迎風面積。圖6所示為正二十面體平均方法示意圖。

圖6 正二十面體平均方法示意圖Fig.6 Schematic diagram of the icosahedron averaging method

破片在彈道飛行時的迎風面積和阻力系數(shù)均與破片姿態(tài)存在一一對應的強依賴關系，因此破片自由翻滾狀態(tài)下平均迎風面積的等效方法也適用于平均阻力系數(shù)。借鑒美國軍用裝備國際試驗操作規(guī)程(ITOP)4-2-813測量平均迎風面積的方法，同樣采用正二十面體法對非球形破片的平均阻力系數(shù)進行等效計算。具體做法如下，首先采用11節(jié)中給出的數(shù)值方法計算得到流場方向沿20個中心- 面心連線和12個中心- 頂點連線的破片阻力系數(shù)，d和d，然后將這些阻力系數(shù)平均，以等效破片隨機翻滾狀態(tài)下的平均阻力系數(shù)，

(3)

圖7 立方體破片模擬結果與試驗結果對比Fig.7 Comparison of numerical results of cubic fragments with experimental data

表2 不同平均方法平均阻力系數(shù)對比Tab.2 Comparison of different averaging methods

圖8比較了正二十面體法計算得到的圓柱形破片(長徑比為1)在不同馬赫數(shù)下的平均阻力系數(shù)和文獻[6-7]的試驗結果。由圖8可見，等效平均阻力系數(shù)曲線可以較好地穿過彈道槍試驗的散點，再次驗證了正二十面體平均方法對非球形破片在隨機翻滾狀態(tài)下平均阻力系數(shù)的等效有效性。

圖8 圓柱體破片模擬結果與試驗結果對比Fig.8 Comparison of numerical results of cylindrical fragments with experimental data

2 破片形狀對阻力系數(shù)影響規(guī)律

為考察破片形狀對隨機翻滾狀態(tài)下破片平均阻力系數(shù)的影響，采用正二十面體法計算得到不同形狀破片的平均阻力系數(shù)。為了擴大結論的適用范圍，破片形狀包括了球形、長方體、圓柱體、三棱柱等規(guī)則形狀，也包括一部分靜爆試驗回收的不規(guī)則破片，圖9展示了其外觀、建模和網格化的過程。為進一步豐富樣本的類型，還通過球諧函數(shù)生成了一系列形狀不規(guī)則的破片外形并引入樣本中，通過球形度描述破片的形狀，所有的破片樣本統(tǒng)計結果如表3所示。圖10所示為馬赫數(shù)01下不同球形度破片的平均阻力系數(shù)。顯然，在破片馬赫數(shù)為01時，平均阻力系數(shù)與球形度存在顯著的負相關。圖10中同時給出了McClescky垂直風洞試驗的結果。值得注意的是，McClescky垂直風洞試驗采用的是非規(guī)則破片，非規(guī)則破片在馬赫數(shù)01時平均阻力系數(shù)與球形度的相關性與各類形狀破片相同，表明球形度是影響破片平均阻力系數(shù)重要的形狀參量，可以作為一種統(tǒng)一的參數(shù)對破片形狀進行描述。

圖11給出了不同球形度破片在跨聲速和超聲速飛行時的平均阻力系數(shù)。由圖11可見：在跨聲速和超聲速飛行時，破片平均阻力系數(shù)與球形度的負相關性仍然存在，但相關性明顯減弱；特別是球形度低于06的破片，平均阻力系數(shù)與球形度的相關性變得非常不顯著。

圖9 爆炸回收非規(guī)則破片建模過程Fig.9 Modeling process of irregular fragment recovered from the detonation

表3 數(shù)值模擬破片外形統(tǒng)計

圖10 不同形狀破片模擬結果與垂直風洞試驗結果比較Fig.10 Comparison of numerical results of different shapes of fragments with vertical wind tunnel experimental data

圖11 不同馬赫數(shù)下阻力系數(shù)與球形度的關系Fig.11 Drag coefficient versus sphericity under different Mach numbers

從圖10、圖11中可以看出，阻力系數(shù)與球形度明顯相關，但二者很難通過簡單的函數(shù)關系進行擬合，且不同飛行速度下球形度和阻力系數(shù)的關系也明顯不同?；诖耍疚牟扇∪斯ど窠浘W絡模型建立阻力系數(shù)綜合預測模型。

3 阻力系數(shù)預測模型的建立

3.1 人工神經網絡模型

由引言可知，影響破片阻力系數(shù)的參數(shù)主要包括破片形狀和飛行馬赫數(shù)，因此本文采用人工神經網絡擬合的方式構建(，)。

人工神經網絡是指在輸入參數(shù)和輸出參數(shù)之間有若干層隱藏的神經網絡，人工神經網絡具有很強的復雜非線性系統(tǒng)建模能力，能夠刻畫難以用解析式表達的高維映射。一個典型人工神經網絡結構如圖12所示，包含一個輸入層、一個輸出層以及兩層隱藏層，每一個節(jié)點都是一個神經元(見圖13)，每個神經元接受來自前一層的神經元不同權重的值，經神經元內激活函數(shù)處理后向下一層傳播。圖12中，、、、分別為神經網絡輸入參數(shù)，、、分別為不同層之間神經元連接的權重向量，、、分別為不同層之間神經元連接的偏移量，(=1，2，3)、分別為第1層、第2層神經元的輸入值，、分別為第1層、第2層神經元的輸出值，為神經網絡模型的輸出值。圖13中，、、、分別表示該神經元接受上一層不同神經元輸出的權重值，表示偏移量，為該神經元的輸入值，為該神經元輸出值，()表示激活函數(shù)。

圖12 人工神經網絡結構Fig.12 Structure of artificial neural network

圖13 神經元Fig.13 Neuron

神經元的輸入為前一層神經元的輸出的線性組合：

(4)

輸入神經元的值經過激活函數(shù)()處理后繼續(xù)向后傳播：

=()

(5)

本文使用Softplus函數(shù)作為神經元的激活函數(shù)，是因為當輸入值>0時Softplus函數(shù)梯度較大，迭代過程可以較快達到收斂，Softplus函數(shù)表達式如(6)式，函數(shù)圖像如圖14所示。

()=ln(1+e)

(6)

圖14 Softplus激活函數(shù)Fig.14 Softplus activation function

(7)

損失函數(shù)是神經網絡參數(shù)、的函數(shù)，因此神經網絡需要迭代找到使損失函數(shù)最小的、，采用的迭代方法是梯度下降算法，損失函數(shù)對神經網絡參數(shù)的梯度通過鏈式求導法則獲得，即神經網絡的向后傳播算法。本文采用Adam優(yōu)化算法對神經網絡進行優(yōu)化，在向后傳播基礎上，以一定學習率不斷更新模型的權重，使其梯度不斷下降，損失函數(shù)的值不斷減小。

3.2 阻力系數(shù)預測模型

本文搭建的人工神經網絡模型，輸入參數(shù)為影響破片阻力系數(shù)的因素，包括球形度以及馬赫數(shù)，輸出參數(shù)為阻力系數(shù)。對數(shù)值模擬計算得到的712個數(shù)據樣本進行擬合，同時為了防止過擬合，隨機選擇其中的600個樣本點進行模型訓練，其余樣本點用于模型驗證。訓練樣本的空間分布如圖15所示。搭建包含3個隱藏層，每個隱藏層包含20個神經元的神經網絡模型。圖16所示為模型訓練過程中訓練集以及驗證集的損失函數(shù)的變化曲線，可以看出迭代到100 000步之后，訓練集和驗證集的損失函數(shù)實現(xiàn)了收斂，損失函數(shù)值均降低到0003以下。表明該神經網絡模型沒有出現(xiàn)過擬合或欠擬合狀態(tài)，擬合效果比較理想。訓練完成后得到的阻力系數(shù)預測模型如圖17所示，訓練得到的不同球形度下阻力系數(shù)與馬赫數(shù)的關系如圖18所示。

圖15 訓練樣本的空間分布Fig.15 Spatial distribution of training samples

圖16 損失函數(shù)變化曲線Fig.16 Changing curves of loss function

圖17 神經網絡阻力系數(shù)模型訓練結果Fig.17 Neural network-based drag coefficient prediction model

圖18 訓練得到的不同球形度下阻力系數(shù)與馬赫數(shù)的關系Fig.18 Drag coefficient versus Mach number for different sphericities trained by ANN

圖17中阻力系數(shù)預測模型可以較好地對(，)進行擬合。任意形狀的破片在任意飛行速度下的阻力系數(shù)預測值可以從圖17所示的結果中得到。圖18為訓練得到的不同球形度下阻力系數(shù)與馬赫數(shù)的關系。從圖18中可以看出，不同球形度破片的阻力系數(shù)隨馬赫數(shù)變化的曲線形狀類似但數(shù)值明顯不同：當破片飛行速度較小時，破片形狀會顯著影響阻力系數(shù)；當飛行馬赫數(shù)大于2時，形狀對阻力系數(shù)的影響變得不明顯。

3.3 模型泛化能力評價

神經網絡的泛化能力是指神經網絡在訓練完成以后神經網絡對測試樣本或工作樣本做出正確反應的能力。在完成模型訓練后，本文進一步計算了3種破片(球形度分別為046、076和091)在超聲速至亞聲速區(qū)間內的阻力系數(shù)作為該預測模型的測試樣本。不同的破片阻力系數(shù)模擬值與模型的預測值的對比如圖19(a)、圖19(b)、圖19(c)所示，圖19(d)為所有測試點的預測值與模擬值的相對誤差。

從圖19中可以看出，神經網絡阻力系數(shù)預測模型可以很好地預測3種破片的阻力系數(shù)曲線，測試結果相對誤差均小于3，表明該模型具有較強的泛化能力。

神經網絡訓練得到的阻力系數(shù)預測模型可用于破片的外彈道軌跡計算中，實現(xiàn)根據破片的實際形狀賦予其相應的阻力系數(shù)曲線。為進一步驗證阻力系數(shù)預測模型的適用性，特別是針對非規(guī)則破片的適用性，將預測模型應用于文獻[11]中非規(guī)則破片的軌跡計算，并與彈道槍實測阻力系數(shù)計算出的軌跡進行對比(見圖20)。該破片為爆炸測試場回收的不規(guī)則破片，經測試，球形度為058，質量為314 g，初始速度1 524 m/s，飛散角度為20°。從圖20中可以看出，使用神經網絡預測模型得到的軌跡計算結果與彈道槍實測結果的吻合程度相當好，最終神經網絡模型計算得到的落點位置為429 m，彈道槍實測阻力系數(shù)結果得到的落點位置為432 m，相對誤差07，計算精度完全可以滿足破片彈道計算的要求。

圖19 不同球形度破片模型預測阻力系數(shù)與數(shù)值模擬結果對比Fig.19 Comparison between the predicted drag coefficient and the numerically simulated results of fragments with different sphericities

圖20 使用不同阻力模型數(shù)值計算的破片運行軌跡對比Fig.20 Trajectories of fragments calculated with different drag models

4 結束語

本文采用正二十面體平均方法對非球形破片高速飛行狀態(tài)的空氣阻力系數(shù)進行數(shù)值模擬，得到形狀以及馬赫數(shù)對破片阻力系數(shù)的影響規(guī)律，通過人工神經網絡建立了預測破片阻力系數(shù)的模型。得出主要結論如下：

1) 各類形狀破片的空氣阻力系數(shù)隨馬赫數(shù)的增加均呈現(xiàn)先增大后逐漸減小的規(guī)律，但數(shù)值明顯不同，形狀對阻力系數(shù)存在明顯影響。

2) 對于非球形破片，正二十面體平均方法是一種估計破片隨機翻滾狀態(tài)下的阻力系數(shù)的有效方法，可應用于任意形狀破片的阻力系數(shù)計算中。

3) 本文建立的人工神經網絡阻力系數(shù)預測模型可以較好地描述馬赫數(shù)以及破片形狀對阻力系數(shù)的影響，且具有良好的泛化能力，可以利用其進行其他非球形破片阻力系數(shù)的預測。