BP神經網絡和圓筒能量模型標定炸藥的JWL參數

崔 浩，郭 銳，顧曉輝，宋 浦，楊永亮，江 琳，俞旸暉

(1.南京理工大學 機械工程學院，江蘇 南京 210094; 2.西安近代化學研究所 燃燒與爆炸技術重點實驗室，陜西 西安 710065)

引 言

爆炸流體動力學程序現已廣泛應用于現代爆炸力學的多個研究領域，包括模擬爆炸、沖擊波、大應變和高應變率等各種復雜高能的物理過程，數值計算的精度和有效性很大程度上取決于所使用炸藥的狀態方程參數。目前存在多種炸藥爆轟產物狀態方程[1-2]，如VLW、Williamsburg、Lennard-Jones-Devonshire (LJD)、Becker-Kistiakowsky-Wilson (BKW)和Jones-Wilkins-Lee (JWL)等，其中精度較高、形式簡單的JWL狀態方程是目前最常用的經驗方程式，并已嵌入諸如AUTODYN、ABAQUS和LS-DYNA等眾多數值仿真軟件中。

LLNL實驗室的J. W. Kury等[3]設計的圓筒試驗用于表征炸藥爆轟時的能量密度，也是目前確定炸藥JWL狀態方程參數最常用的標定試驗之一。為了獲取高精度的JWL參數，通常使用流體動力學軟件模擬炸藥膨脹和驅動圓筒壁的過程，同時不斷手動調整JWL參數組合，直到仿真結果和試驗數據相吻合[4-6]。但由于JWL狀態方程中包含6個未知參數，不僅參數組合的復雜程度會導致手動調整過程變得繁瑣和費時，而且有時必須進行微小的調整和一定程度的直覺才能滿足參數組合精度要求，這增加了參數調整的不確定性。雖然CHEETCH、CHEQ等[7]專用標定參數軟件可用于快速確定JWL參數，但關于軟件方法未見公開。雖然Mortensen等[8]通過大量的數值仿真和圓筒試驗數據發現了JWL參數之間存在近似的線性關系，但通過這種線性關系確定的JWL參數精度有限，且線性關系無法適用于所有炸藥。此外，眾多學者還開展了利用智能算法求解JWL狀態方程參數的研究，其中王成等[9]提出了基于遺傳算法和γ律狀態方程的JWL參數計算法，發現擬合確定的JWL狀態方程p—V曲線比γ律狀態方程p—V曲線更接近標準的JWL狀態方程p—V曲線；楊晨琛等[10]開展了水下爆炸試驗，發展出了根據水汽界面確定爆轟產物JWL狀態方程的遺傳算法，獲得了多種炸藥的JWL狀態方程的近似全局最優解；Cui等[11]采用遺傳算法擬合得到了幾種理想炸藥的JWL狀態方程參數，將擬合參數代入數值仿真中顯示仿真結果和試驗數據相吻合。但由于遺傳算法中基因操作涉及到染色體適應度值的計算，可能會導致計算效率相對較低。

為此，本研究基于BP神經網絡和圓筒能量模型，提出了一種可以根據圓筒試驗數據快速準確標定炸藥JWL參數的三點標定法，即根據圓筒壁膨脹過程中特定3個點(爆轟產物相對體積為2.4、4.4和7.0)的數據進行標定的方法。三點標定法首先對BP神經網絡進行訓練，使其可以預測圓筒能量模型，隨后采用遺傳算法(Genetic Algorithm，GA)尋找訓練好的BP神經網絡的最優解。三點標定法將BP神經網絡強大的擬合預測功能與GA優秀的尋優能力相結合，省去了求解非線性守恒方程組和計算染色體適應度值的過程，提高了求解速度。

1 圓筒試驗和JWL狀態方程

1.1 圓筒試驗

圓筒試驗示意圖如圖1(a)所示，在圓筒試驗中，平面波發生器從一端起爆裝填在銅管內的炸藥，當爆轟波陣面通過觀察窗口時，銅管在爆轟產物的驅動下向外膨脹，銅管的徑向位移將遮蔽背光(由氬氣閃光彈提供照明)，高速掃描照相機或激光干涉儀則記錄下銅管的位移歷史。圖1(b)為高速掃描照相機記錄的爆轟產物驅動下膨脹的圓筒壁[12]。

圖1 圓筒試驗示意圖及高速掃描照相機圖Fig.1 Experimental set-up and camera image of the cylinder test

通過圓筒試驗可獲得一組圓筒外徑位移隨時間變化的離散點，為了計算圓筒外徑的位移速度和炸藥的作功能力，通常將這組離散點擬合成函數形式。位移函數形式多種多樣，包括但不限于以基于爆轟產物壓強指數下降假設的函數[13]、考慮沖擊波和爆轟產物作用的函數[14]、多項式函數[15]等，將這些位移函數微分即可獲得圓筒壁外壁徑向膨脹速度，可用于表征炸藥作功能力和計算炸藥JWL參數等。

1.2 JWL狀態方程

等熵條件下，JWL狀態方程表達式為：

p=Ae-R1V+Be-R2V+CV-(ω+1)

(1)

式中：p為爆轟產物壓力；V為爆轟產物相對比容；A、B、C、R1、R2、ω為6個常數。

在絕熱條件下，對式(1)積分可得爆轟產物內能為：

(2)

炸藥理想爆轟時，根據C-J條件和Hugoniot關系，JWL參數之間須滿足以下約束守恒方程組：

(3)

式中：ρ0為炸藥初始密度；D為炸藥爆速；VCJ為C-J點處爆轟產物的相對比容；pCJ為炸藥爆壓；E0為單位體積炸藥初始內能。由守恒方程組(3)可知，JWL狀態方程中的6個參數只有3個是獨立的，只需確定其中的3個參數便可確定剩余JWL參數。

2 圓筒能量模型

圓筒能量模型[16-18]描述了炸藥爆轟產物驅動圓筒壁時的能量分布和能量轉換關系，Sanchidrián[19]和Castedo[20]利用圓筒試驗和圓筒能量模型獲得了多種乳化炸藥的JWL方程參數，發現數值計算結果和試驗數據高度吻合，證明了可利用圓筒能量模型確定高精度的JWL參數。圖2為圓筒能量模型示意圖，其中θ為圓筒外壁和軸線的夾角；圓筒初始內半徑和壁厚分別為R0和x0，任意圓筒膨脹位移r時的圓筒內半徑和壁厚分別為R和R+x。

圖2 圓筒壁膨脹圖Fig.2 Diagram of the cylinder expansion

在炸藥爆轟產物膨脹過程中，其內能不斷轉換成圓筒壁和爆轟產物的動能，兩者的動能之和為：

(umcosθ)2

(4)

式中：Ed為單位體積炸藥爆轟產生的動能；ρm為圓筒壁密度；um為圓筒外壁移動速度。

圓筒能量模型假定初始單位長度半徑為R0的薄切片炸藥，其爆轟產物體積是以Rθ為母線的圓錐側面，因此當圓筒膨脹到內半徑為R時的爆轟產物相對比容為：

(5)

圓筒外徑徑向位移為：

r=(R+x)-(R0+x0)

(6)

假定圓筒密度在膨脹過程中保持不變，根據質量守恒可得：

(7)

結合式(6)、式(7)，爆轟產物的相對比容還可表示為：

(8)

除此之外，圓筒壁角θ可通過下式求出：

(9)

除此之外，還需考慮圓筒變形功，尤其是圓筒變形功在系統總能量的占比會隨著爆轟產物膨脹率的增大而增大[21]。圓筒膨脹內徑為R時的圓筒壁變形功可由下式求出：

(10)

式中：σf為銅管的流動應力，其值為89.63MPa。在圓筒試驗中，爆轟產物的能量轉換成圓筒壁和爆轟產物的動能以及圓筒壁的變形功，根據能量守恒定律可得：

E0-Es=Ed+Wdef

(11)

因此，定義圓筒試驗能量模型的能量差為:

f(V)=E0-Es-Ed-Wdef

(12)

理想情況下，爆轟產物膨脹過程中的系統能量差應始終等于或接近0。在實際應用中，當爆轟產物相對體積為2.4、4.4和7.0時的能量差接近0即可判定此組JWL參數適用于此炸藥[18,22]。因此定義圓筒模型的能量差函數為：

FE=exp[-0.5(f2(2.4)+f2(4.4)+f2(7.0))]

(13)

觀察上式可知，當FE值越接近1時表明系統能量差越小，此組JWL參數越適用于該炸藥。對于特定的炸藥，為了計算某一組JWL參數確定的圓筒能量差函數值，通常僅需確定JWL狀態方程中的3個變量即可(守恒方程組(3)使得JWL參數只有3個獨立變量)。

3 BP神經網絡預測圓筒能量模型

由前述可知，由于一組JWL參數的自變量為3個，為了計算一組JWL參數的FE值(GA中染色體適應度函數同樣選取FE)，需要求解一次非線性方程組，而由于GA中通常涉及大量的個體適應度值計算，包括種群進化過程中的基因操作過程需要根據個體適應度值來進行個體選擇、判斷交叉和變異概率等，每一次GA基因操作都需要求解大量的非線性方程組和計算各個染色體的適應度值，這無疑大大增加了計算量，導致計算效率較低。本研究期望通過訓練BP神經網絡擬合圓筒能量模型，可以直接預測任意一組JWL參數相對應的FE值，提高求解效率。

3.1 BP神經網絡結構

如圖3所示，為了改善網絡性能，提高辨識精度，本研究采用雙隱含層的網絡結構，包括單輸入層、雙隱含層和輸出層[23]。由上文可知僅需確定JWL狀態方程中的3個變量便可求出對應的FE值，由于R1、R2和ω的取值范圍較小，因此選取R1、R2和ω作為自變量，并設置其取值范圍分別為4～5、1～2和0.2～0.4[24]。本研究構建的神經網絡的輸入層節點數為3(自變量R1、R2、ω)、輸出層節點數為1，根據試算法確定第一、第二隱含層的節點數分別為20和15。

圖3 BP神經網絡結構圖Fig.3 Diagram of the BP neural network structure

雙隱含層BP神經網絡的訓練過程具體如下：

(1)網絡初始化，網絡參數(權值和閾值)隨機賦初值。

(2)信息正向傳播，計算各層的輸入和輸出。

設神經網絡有P個訓練樣本，對任意訓練樣本Xk[xk1,xk2,xk3]，期望輸出為T[tk]。計算各層神經元節點輸入、輸出分別為：

(14)

vj=f(uj)

(15)

式中：Wij為節點i與節點j之間的權值，j為與i節點相鄰的下一層的節點；bj為節點j的閾值；f為節點傳遞函數，本研究中的輸入層到第一隱含層、第一隱含層到第二隱含層神經元的傳遞函數均采用Log-Sigmoid函數，第二隱含層到輸出層的神經元傳遞函數則采用purnlin線性函數。

(3) 計算期望輸出與目標輸出之間的平方誤差為：

(16)

式中:ek為實際輸出與期望輸出的差值；n為迭代次數。

(4) 誤差反向傳播，計算節點i與節點j之間的權值修正量為：

(17)

由此可得到下一次迭代時的修正權值為：

Wij(n+1)=Wij(n)+ΔWij(n)

(18)

式中：η為學習步長。

(5)記憶訓練

重復步驟(2)～(4)，正向傳播和反向傳播交替進行，不斷逐層修正權值和閾值直至網絡達到精度要求后輸出網絡。

3.2 訓練結果

隨機生成2000個訓練樣本，訓練函數選用trainlm，網絡均方差精度設置成1×10-10。以TNT為例，訓練結束后取前100個樣本的期望輸出和實際輸出進行對比，如圖4所示。觀察圖4可發現，實際輸出與期望輸出基本重合，表明神經網絡訓練效果較好，可以把網絡預測輸出近似看成實際輸出。

圖4 期望輸出和實際輸出對比Fig.4 Comparison of expected outputs and actual outputs

4 BP-GA程序求解炸藥JWL參數

訓練好的BP神經網絡可以較好地擬合圓筒能量模型，并且能出快速準確地給出任意一組R1、R2和ω確定的FE值，但由于神經網絡不具備尋優功能，無法僅憑神經網絡尋求到最優的JWL參數組合，因此通常還需借助GA優秀的搜索尋優能力[25]。

4.1 GA步驟

本研究選取的染色體編碼方式、選擇機制、交叉算子和變異算子等均同文獻[11]，具體步驟如下：

(1) 選擇，從種群中選擇父代染色體進入配對庫；

(2) 交叉，配對庫中的父代兩兩交叉產生下一代；

(3) 變異，交叉生成的子代染色體以一定概率發生變異；

(4) 種群替換，執行步驟(1)～(3)后種群進化了一代，重復此步驟直到種群連續50代無進化。

4.2 BP-GA程序計算流程

三點標定法采用BP-GA程序標定炸藥JWL參數，具體流程圖如圖5所示。

圖5 BP-GA標定JWL參數流程圖Fig.5 Flow chart of JWL parameters calibrated by BP-GA

三點標定法首先創建BP神經網絡模塊，輸入圓筒的ρm、R0、x0、σf和炸藥的ρ0、D、E0、PCJ以及圓筒試驗中爆轟產物相對體積為2.4、4.4和7.0時的圓筒壁速度后隨機生成2000個樣本進行訓練。當訓練樣本均方差低于1×10-10時訓練完成，訓練結束后的BP神經網絡即可直接預測任意一條染色體(R1、R2、ω)的適應度值。BP神經網絡訓練結束后啟動GA模塊尋求訓練好的BP神經網絡的最優解，種群經過多代進化后進化成熟，挑選成熟后的種群中適應度值最大的個體作為最優解，此最優解即為炸藥JWL參數。

4.3 算例和計算效率

在微機CPU@3.60GHz，8.00GB內存的環境下，采用Matlab軟件進行編程。以TNT炸藥為算例，分別運行文獻[11]中的GA程序和本研究的BP-GA程序，其中GA程序均設置種群數量為60以及種群連續50代無進化為停止標準，而BP程序設置生成2000個訓練樣本。結果顯示文獻[11]中的GA程序單次運行耗時為8425.999s；而BP-GA程序單次運行所用總時間為396.019s，其中生成樣本耗時328.920s，訓練樣本耗時9.736s，GA尋優耗時57.363s，效率提升了大約21倍。

5 尋優結果及數值檢驗

5.1 尋優結果

將表1中4種炸藥的爆轟特性和文獻[26]中的圓筒試驗數據輸入BP-GA程序中進行計算，標定了4種常用炸藥的JWL方程參數，并和文獻[27]中采用圓筒試驗法得到的JWL參數進行對比，如表2所示。

表1 炸藥爆轟特性Table 1 Detonation properties of the explosives

表2 炸藥JWL狀態方程參數對比Table 2 Parameters of JWL EOS of explosives

文獻[27]中JWL參數和BP-GA程序標定的JWL參數所確定的p—V曲線之比如圖6所示，其中將文獻[27]確定的p—V曲線視為標準曲線。兩者的匹配度采用確定系數R2進行定量評估，該系數為擬合數據和標準數據的殘差之比。觀察圖6可發現，采用BP-GA程序標定得到的p—V曲線和文獻中標準p—V曲線相差較小，說明三點標定法的精度較高。

圖6 計算JWL狀態方程和標準JWL狀態方程的p—V曲線對比Fig.6 Comparison of p—V curves between calculated JWL EOS and standard JWL EOS

5.2 數值檢驗

采用非線性動力學分析軟件AUTODYN-2D對表2參數進行數值檢驗。如圖7所示，建立了標準圓筒試驗(圓筒內徑為Φ25.4mm)的二維軸對稱模型，其中圓筒壁采用拉格朗日算法以描述金屬運動特性，炸藥采用歐拉算法以模擬爆轟產物的膨脹過程。采用歐拉算法建立了尺寸為405mm×254mm的空氣域，炸藥和空氣的網格大小(周向×徑向)設置為0.5mm×0.488mm，圓筒壁網格大小設置為0.5mm×0.48mm。

圖7 圓筒試驗計算模型Fig.7 Calculation model of cylinder test

在仿真中設置拉格朗日算法和歐拉算法耦合，空氣域表面設置成無反射邊界。此外，與標準試驗中一樣，在距離起爆端一定距離的外壁表面設立一個高斯點，以記錄圓筒壁的徑向位移和膨脹速度。

炸藥材料采用JWL狀態方程，銅管材料選取密度為8.93g/cm3的OHFC無氧銅并采用Steinberg-Guinan強度模型，其剪切模型G和屈服應力Y與有效塑形應變ε、壓力p和溫度T的函數關系式為[28]：

(19)

Y=Y0(1+βε)nG/G0

(20)

須符合受限條件：

Y0(1+βε)n≤Ymax

(21)

式中：η為壓縮率；β和n分別表示硬化常數和硬化指數。下標0表示參考狀態(T=300K，p=0，ε=0)下G和Y的值，而帶下標p和T的素數參數是該參數相對于參考狀態下的壓力和溫度的導數。上述參數值均取自AUTODYN材料庫(見表3)[29]。此外，采用Mie-Grüneisen狀態方程[30]描述銅管在高壓的狀態，取沖擊Hugoniot參數c=3940m/s，Grüneisen系數為2.02。

表3 OHFC無氧銅材料特性(Steinberg-Guinan模型)Table 3 Properties of the OHFC copper (Steinberg-Guinan)

圖8和圖9分別為數值計算得到的圓筒壁徑向位移和速度變化曲線，和試驗結果對比可發現，對于這4種炸藥，仿真結果較好地重現了試驗數據，證明了BP-GA程序標定的JWL參數的有效性和高精度，能夠滿足工程應用需求。

圖8 徑向位移時程曲線對比圖Fig.8 Time-history curves of radial displacement

圖9 圓筒外壁移動速度對比圖Fig.9 Comparison of velocities of cylinder outer wall

5.3 不同尺寸圓筒試驗的適用性

上述數值仿真和試驗結果的一致性證明了BP-GA程序處理標準圓筒試驗數據的高精度。根據藥柱直徑的不同，圓筒試驗通常分為Φ25.4mm標準圓筒試驗和Φ50mm圓筒試驗，為了探討本方法是否同樣適用于Φ50mm圓筒試驗，采用三點標定法對文獻[31]中某RDX基PBX炸藥的Φ50mm圓筒試驗數據進行處理，計算得到其JWL參數為：A=650.08GPa、B=8.994GPa、C=0.937GPa、R1=4.441、R2=1.086、ω=0.400。將上述參數代入仿真程序中進行Φ50mm圓筒試驗工況計算，除圓筒和炸藥模型尺寸不同外，網格大小和仿真條件設置均同上一節。圖10為仿真結果和試驗數據對比，同樣地，仿真得到的圓筒外壁位移和速度曲線均和試驗相吻合，證明了三點標定法可適用于兩種不同直徑的圓筒試驗。

圖10 50mm圓筒試驗數值模擬與實驗結果對比Fig.10 Comparison of numerical simulation and experimental results of 50mm cylinder test

6 結 論

(1)基于圓筒試驗中3個特定點的數據和圓筒能量模型，隨機生成2000個樣本對雙隱含層結構的BP神經網絡進行訓練，訓練好的BP神經網絡可較好地預測任意一組R1、R2和ω代表的能量差值。

(2) 三點標定法采用GA尋求訓練好的BP神經網絡的最優值，此種方法結合了BP神經網絡強大的擬合預測功能與GA優秀的尋優能力，省去了求解非線性守恒方程組和計算染色體適應度值的過程，算例表明相比GA算法，BP-GA程序的計算效率大約提高了21倍。

(3)將4種炸藥的Φ25.4mm標準圓筒以及一種炸藥的Φ50mm圓筒試驗數據代入BP-GA程序中進行計算，將獲得的JWL狀態方程參數代入數值仿真中，發現仿真和試驗吻合較好，表明三點標定法可適用于兩種不同直徑的圓筒試驗。