基于人工神經網絡的非球形破片阻力系數預測模型

辛大鈞， 薛琨

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

0 引言

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

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

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

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

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

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

1 數值模擬方法及驗證

1.1 球形破片阻力系數模擬方法

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

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

圖1 數值模擬Fig.1 Numerical modeling

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

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

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

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

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

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

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

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

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

1.2 非球形破片阻力系數計算方法

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

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

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

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

(1)

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

(2)

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

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

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

(3)

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

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

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

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

2 破片形狀對阻力系數影響規律

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

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

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

表3 數值模擬破片外形統計

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

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

從圖10、圖11中可以看出，阻力系數與球形度明顯相關，但二者很難通過簡單的函數關系進行擬合，且不同飛行速度下球形度和阻力系數的關系也明顯不同。基于此，本文采取人工神經網絡模型建立阻力系數綜合預測模型。

3 阻力系數預測模型的建立

3.1 人工神經網絡模型

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

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

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

圖13 神經元Fig.13 Neuron

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

(4)

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

=()

(5)

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

()=ln(1+e)

(6)

圖14 Softplus激活函數Fig.14 Softplus activation function

(7)

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

3.2 阻力系數預測模型

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

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

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

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

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

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

3.3 模型泛化能力評價

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

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

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

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

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

4 結束語

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

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

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

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