崔 浩, 郭 銳，宋 浦，顧曉輝，周 昊，邢柏陽

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

炸藥爆轟產(chǎn)物的狀態(tài)方程是描述爆轟波陣面過后爆轟產(chǎn)物的壓力、溫度和密度等物理量之間關系的復雜函數(shù)，目前主要使用經(jīng)驗或半經(jīng)驗的狀態(tài)方程式。其中，JWL(Jones-Wilkins-Lee)狀態(tài)方程[1]可以比較精確地描述炸藥爆轟產(chǎn)物膨脹驅動的過程，Kury等[2]提出的圓筒試驗可幫助確定該方程中的未知參數(shù)。炸藥爆炸產(chǎn)物JWL狀態(tài)方程參數(shù)的確立具有重要意義，是開展爆效毀傷、武器設計和安全防護等工作的基礎。

關于確定炸藥爆炸產(chǎn)物JWL狀態(tài)方程參數(shù)的方法學者們已進行了大量的研究，大多采用圓筒試驗和數(shù)值仿真相結合的方法，通常需要不斷手動調整參數(shù)組合直到仿真和試驗數(shù)據(jù)誤差小于1%[3-6]，但此種方法具有試驗成本高、參數(shù)獲取效率較低等缺陷。此外，還有將圓筒試驗數(shù)據(jù)代入理論模型來確定狀態(tài)方程參數(shù)的方法[7-9]，首先根據(jù)試驗數(shù)據(jù)計算出炸藥的一些爆轟參數(shù)，隨后代入守恒方程中依次計算JWL狀態(tài)方程中的低壓、中壓和高壓貢獻項的參數(shù)，此操作一定程度上提高了炸藥JWL參數(shù)的獲取效率。但由于圓筒試驗存在成本高、操作繁雜、數(shù)據(jù)讀取費時等缺陷，于是便有學者提出了理論擬合JWL參數(shù)的方法[10]。其中，趙崢等[11]提出了一種確定炸藥JWL參數(shù)的γ擬合法，并采用差分進化法得到了4種炸藥的JWL狀態(tài)方程參數(shù)，發(fā)現(xiàn)擬合JWL狀態(tài)方程和圓筒試驗法得到的JWL狀態(tài)方程的P-V曲線相差較??；王成等[12]基于γ律狀態(tài)方程和遺傳算法擬合了5種炸藥的JWL參數(shù)，并將其代入數(shù)值計算中進行檢驗，發(fā)現(xiàn)擬合JWL參數(shù)具有較高的精度；沈飛等[13]基于圓筒試驗中的能量轉換關系，建立了一種確定炸藥JWL參數(shù)的簡易計算流程，發(fā)現(xiàn)此種方法確定的JWL狀態(tài)方程的P-V曲線能夠滿足工程應用的需要。

本文結合遺傳算法優(yōu)秀的搜索能力和圓筒能量模型，提出了一種僅需已知炸藥的初始密度和爆速，便可通過自編程遺傳算法程序辨識炸藥JWL狀態(tài)方程參數(shù)的方法，而炸藥的初始密度和爆速通過簡單測量即可獲得。同時對炸藥驅動圓筒壁的過程進行數(shù)值計算，以檢驗辨識得到的JWL狀態(tài)方程參數(shù)的有效性。

1 JWL狀態(tài)方程及圓筒能量模型

1.1 JWL狀態(tài)方程

等熵形式下的JWL狀態(tài)方程為

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

(1)

式中：P為爆轟產(chǎn)物壓力；V為爆轟產(chǎn)物相對比容，其值為爆轟產(chǎn)物體積和炸藥初始體積之比；A、B、C、R1、R2、ω為6個待定參數(shù)。式(1)中右端依次是爆轟產(chǎn)物膨脹過程中高壓、中壓和低壓3個階段的貢獻項。

由熱力學關系對P積分，可得等熵線上的內(nèi)能為

(2)

炸藥理想爆轟時，根據(jù)C-J條件和Hugoniot關系式可得到如下3個相容方程

根據(jù)C-J條件-(?ps/?V)VCJ=ρ0D2可得

(3)

由爆轟產(chǎn)物狀態(tài)的Hugoniot關系式可得

(11)

因C-J等熵線通過C-J點，則

(4)

式中：ρ0為炸藥初始密度；D為炸藥爆速；VCJ為C-J點處爆轟產(chǎn)物的相對比容；PCJ為炸藥爆壓；E0為單位體積炸藥初始內(nèi)能。相容方程組中VCJ和PCJ的值可由下式確定

(5)

(6)

式中,γ為炸藥的絕熱指數(shù)(或多方指數(shù))，其值可由炸藥初始密度ρ0、爆速D、爆熱Q或炸藥初始體積能量E0[14-15]確定

(7)

1.2 圓筒能量模型

Gurney模型[16-17]描述了受限炸藥在無限長圓筒內(nèi)爆炸時的能量分布，但以前提出的Gurney模型假設圓筒壁在膨脹過程中像剛體一樣向外移動，圓筒壁厚度保持不變且移動速度相同，忽視了圓筒壁的塑性形變和厚度變化。為此，本文采用改進的Gurney模型[18-19]，此模型描述圓筒能量密度的精確性經(jīng)過了多種乳化炸藥的圓筒試驗驗證[20]，示意圖如圖1所示。

圖1 炸藥爆轟產(chǎn)物驅動圓筒壁膨脹圖

如圖1所示，在圓筒試驗中，從一端起爆銅管內(nèi)的炸藥，爆轟波以爆速D在炸藥中傳播，爆轟波過后炸藥轉換為爆轟產(chǎn)物并向外膨脹，圓筒壁在爆轟產(chǎn)物的驅動下向外移動，高速掃描照相機記錄下觀測點處的圓筒壁徑向位移。圖1中L為圓筒壁徑向位移，um為圓筒外壁徑向移動速度，θ為圓筒壁角(初始時θ為0°，隨著圓筒壁的移動逐漸增大)。R0為圓筒初始內(nèi)半徑，x0為圓筒初始壁厚，隨著圓筒壁向外移動，圓筒內(nèi)半徑變?yōu)镽，圓筒壁厚度變?yōu)閤。

在爆轟產(chǎn)物膨脹過程中，其內(nèi)能不斷轉換成圓筒壁和爆轟產(chǎn)物的動能，兩者的動能之和通過如下方程求出

(8)

式中,ρm為圓筒材料密度，在標準圓筒試驗中其值為8.93 g/cm3。

對于單位長度炸藥，初始體積是以R0為半徑的薄切片，而轉化為爆轟產(chǎn)物后的體積則是以Rθ為母線、R為半徑的圓錐的側面。此時爆轟產(chǎn)物的相對比容為

(9)

圓筒外徑徑向位移可表示為

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

(10)

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

(11)

結合式(10)和式(11)，式(9)可表示為

(12)

在Φ25.4 mm圓筒試驗[21]中R0=12.7 mm,x0=2.6 mm。此外，文獻[22]給出了計算CHNO型炸藥初始體積能量的經(jīng)驗公式

(13)

式中,Ed(7.0)為爆轟產(chǎn)物相對體積為7.0時的總動能。在圓筒試驗中，炸藥驅動圓筒壁移動是一個絕熱的過程，隨著圓筒壁的移動，爆轟產(chǎn)物的內(nèi)能僅轉換為圓筒壁和爆轟產(chǎn)物的動能，定義圓筒試驗能量模型的能量差為

(14)

對于一組JWL參數(shù)，若在爆轟產(chǎn)物驅動圓筒壁移動的各個時期，f(V)的值均接近0，則說明此組JWL參數(shù)適用于該炸藥。然而在Φ25.4 mm圓筒試驗中，一般只關注圓筒壁徑向位移為6 mm和19 mm時的圓筒能量，按式(12)計算出這兩處的爆轟產(chǎn)物相對體積分別為2.4和7.0，因此在實際計算中只需判斷f(2.4)和f(7.0)是否接近0即可。此外，Castedo等研究發(fā)現(xiàn)圓筒位移6 mm和19 mm時的圓筒壁角分別為10.0°和11.6°，沈飛等的研究和文獻[23]給出了計算CHNO型炸藥在這兩個位置處的圓筒外壁徑向移動速度經(jīng)驗公式

(15)

(16)

采用本節(jié)計算方法得到了4種炸藥的爆轟參數(shù)，和Urtiew等研究中的試驗數(shù)據(jù)進行對比，如表1所示。觀察表1可發(fā)現(xiàn)，除TNT的E0值外，計算結果和試驗數(shù)據(jù)相差較小，兩者的最大誤差不超過2%。

表1 炸藥爆轟參數(shù)對比

2 遺傳算法辨識JWL狀態(tài)方程參數(shù)過程

遺傳算法[24]于1975年由Holland教授等創(chuàng)立，是一種基于自然選擇和基因遺傳學原理的隨機并行搜索算法，且而不需要任何初始化信息即可尋求到全局最優(yōu)解[25]。由于遺傳算法不受搜索空間的約束，對目標函數(shù)沒有連續(xù)、可導或單峰的要求，早已廣泛應用于確定材料物態(tài)方程參數(shù)的研究[26]。

不同JWL參數(shù)組合代表了不同的圓筒模型能量差，本文的目的是尋求能使能量差最小的JWL參數(shù)組合。由于JWL方程含有3個獨立的自變量，僅參數(shù)組合的多樣性就使得尋優(yōu)過程變得困難，而遺傳算法被認為是解決具有較大、復雜搜索空間問題的最適合的搜索方法之一。本文希望借助遺傳算法優(yōu)秀的搜索尋優(yōu)能力尋求最優(yōu)的JWL參數(shù)組合，為了解決遺傳算法中的“早熟”現(xiàn)象及收斂速度慢的問題，對選擇、交叉和變異算子進行了優(yōu)化，同時采用了具有全局收斂性的精英保存策略。

2.1 染色體和適應度函數(shù)

每個染色體對應一個確定的解，本文采用求解精度更高、運算效率更快的浮點數(shù)編碼編輯染色體的基因信息。由3個相容方程可知，JWL狀態(tài)方程中的6個參數(shù)只有3個是獨立的，因此每條染色體攜帶3條基因。為了方便計算，本文選取取值范圍確定的R1、R2和ω作為自變量，并設置R1、R2和ω的取值范圍分別是4～5、1～2和0.2～0.4[27]。

染色體的適應性或一組JWL參數(shù)所確定的圓筒能量差，通過式(17)所示的適應度函數(shù)進行評估

F=|f(2.4)|+|f(7.0)|

(17)

由式(17)可知，當兩處的能量差絕對值之和越小時，適應度函數(shù)值越小，表明染色體的適應性越好，染色體所攜帶的遺傳信息越符合炸藥真實的JWL狀態(tài)方程。為了計算適應度函數(shù)值的大小，需在適應度子程序中先求出每條染色體(R1、R2、ω)對應的A、B、C，隨后代入式(14)中分別求出特定位置的能量差。

2.2 基因操作步驟

遺傳算法從一個由隨機生成的染色體組成的種群開始，通過應用一系列的選擇、交叉和變異等基因操作，使種群朝著適應性更好的染色體方向發(fā)展，本文使用的遺傳算法詳細步驟如下。

步驟1種群初始化。在自變量取值范圍值隨機生成80個染色體(個體)，80個染色體組成了初始種群，染色體數(shù)量的選定根據(jù)以往關于遺傳算法的工作經(jīng)驗。

步驟2選擇算子。通過選擇操作，80個父代染色體將被挑選入配對庫中進行交叉、變異等基因操作。本文采用錦標賽選擇機制，每次從種群中隨機取出sn個染色體(sn值通常取2)，然后選擇其中適應性最好的一個進入配對庫，重復此操作直到挑選出80個父代染色體。

錦標賽選擇機制可以使適應性較好的個體具有較大的“生存”機會，并且由于只使用適應性作為選擇的標準，不受適應度函數(shù)值的影響，因此避免了超級個體的影響，能在一定程度上避免過早收斂和種群停滯現(xiàn)象[28]。

步驟3交叉算子。交叉操作即為配對庫中的80個父代染色體兩兩配對并按照某種方式互換基因以產(chǎn)生下一代。本文的交叉機制采用算數(shù)交叉，其基本思想是兩個父代染色體線性組合產(chǎn)生出兩個新的個體，子代染色體的產(chǎn)生機制如式(18)所示

(18)

步驟4變異算子。為了減少種群多樣性的損失，降低陷入局部最優(yōu)解的風險，交叉生成的子代染色體將以一定的概率發(fā)生變異。本文采用的高斯變異機制的基本流程為：當隨機數(shù)小于變異概率時，隨機選取染色體的某一條基因進行變異，用符合均值為μ，方差為σ2的正態(tài)分布的隨機數(shù)來替換選取的基因值[29]。那么，變異后的基因值為

(19)

式中,μ=(xmin+xmax)/2,σ=(xmax-xmin)/6,xmin和xmax分別為變異基因值的上下限。

步驟5精英保存策略。為了防止選擇、交叉、變異等遺傳操作的隨機性將當前群體中最好的個體破壞掉，本文采用具有全局收斂性的精英保存策略。當前群體中適應性最好的個體不參與交叉運算和變異運算，而是用它來替換掉本代群體經(jīng)過交叉、變異等基因操作后產(chǎn)生的適應性最差的個體。該策略的實施是保證遺傳算法收斂的重要條件之一。

步驟6種群替換。執(zhí)行步驟2～步驟5后種群進化了一代，形成了新的種群，重復此步驟直到種群連續(xù)30代在適應性上沒有任何提高，連續(xù)30代無進化則認為種群進化成熟。

2.3 算法流程及進化過程分析

遺傳算法具體的操作流程如圖2，輸入炸藥初始密度和爆速后，可根據(jù)式(15)和式(16)計算出u6 mm和u19 mm，并代入式(13)中計算出E0，然后根據(jù)式(7)計算出γ的值，隨后可通過式(5)、式(6)計算出VCJ和PCJ的值。最后通過基因操作種群進化成熟后挑選出最優(yōu)個體，讀取最優(yōu)個體攜帶的遺傳信息獲得R1、R2、ω的值并代入相容方程中計算A、B、C的值，得到的參數(shù)組合即為炸藥的JWL參數(shù)組合。

圖2 基因遺傳算法辨識JWL狀態(tài)方程參數(shù)流程圖

圖3為TNT炸藥爆轟產(chǎn)物內(nèi)能隨相對體積的變化曲線，觀察圖中可發(fā)現(xiàn)JWL爆轟產(chǎn)物內(nèi)能方程中3項的占比和相互之間關系隨相對體積而變化，為了確保JWL參數(shù)在物理上的合理性，對JWL參數(shù)添加如下限定條件：

圖3 TNT炸藥爆轟產(chǎn)物比內(nèi)能分布

(1)當爆轟產(chǎn)物相對體積大于某一特定值(通常為V>7[20，22])時，前兩項之和必須小于第三項，即

(20)

(2)當爆轟產(chǎn)物的相對體積大于2時，第一項要小于第二項

(21)

(3)同樣地，高壓時，JWL方程第一項必須大于第二項(通常選取C-J狀態(tài))

(22)

圖4為4種炸藥最優(yōu)個體的適應度函數(shù)值隨進化次數(shù)的變化曲線。觀察圖4可知，種群進化成熟后，4種炸藥最優(yōu)個體的適應度值均接近0，說明進化后圓筒模型能量差很小。由于采用了精英保存策略，種群中最優(yōu)個體在進化過程中無“退化”現(xiàn)象，遺傳算法收斂。

圖4 最優(yōu)個體的適應度值隨進化次數(shù)變化曲線

3 遺傳算法辨識結果及數(shù)值檢驗

3.1 遺傳算法辨識結果

在僅已知炸藥密度和爆速的前提下，利用遺傳算法程序辨識得到了4種常用炸藥的JWL狀態(tài)方程參數(shù)，表2列出了辨識JWL參數(shù)和標準JWL參數(shù)的對比。其中標準JWL參數(shù)取自AUTODYN材料庫,均是通過圓筒試驗標定，其確定的JWL狀態(tài)方程可認為是炸藥的標準JWL狀態(tài)方程。4種炸藥的JWL狀態(tài)方程P-V曲線對比如圖5所示，觀察圖5可發(fā)現(xiàn)，對于這4種炸藥，兩條P-V曲線偏差均較小，為了定量分析辨識曲線和標準曲線的相似度，引入決定系數(shù)R2

圖5 辨識JWL狀態(tài)方程和標準JWL狀態(tài)方程的P-V曲線對比

(23)

式中：SSE為誤差平方和；SST為試驗數(shù)據(jù)和均值之差平方和；R2的取值范圍為[0,1]，其值越接近1，表明辨識曲線越接近標準曲線。計算發(fā)現(xiàn)TNT、HMX、Comp B和PBX-9404的辨識JWL狀態(tài)方程P-V曲線的R2值分別為0.997 7、0.993 8、0.999 3和0.998 6，其值均大于0.99，說明辨識JWL狀態(tài)方程和標準JWL狀態(tài)方程吻合較好。觀察表2可發(fā)現(xiàn)，采用本文方法辨識得到的6個JWL參數(shù)均不同于標準JWL參數(shù)，但是兩者定義的P-V曲線卻非常接近，這是因為JWL狀態(tài)方程的多項式和非線性允許不同參數(shù)組合可以同時高精度地定義同一個方程。

表2 炸藥JWL狀態(tài)方程參數(shù)對比

趙錚等的研究采用差分進化法擬合得到了4種炸藥的JWL參數(shù)，其中有兩種炸藥與本文辨識的炸藥相同，即TNT和HMX。將采用本文遺傳算法辨識得到的JWL方程和趙錚等研究中的差分進化法擬合的JWL方程的P-V曲線進行對比，結果如圖6所示。觀察圖6可發(fā)現(xiàn)，對于這兩種炸藥，采用遺傳算法辨識得到的JWL方程P-V曲線更接近標準JWL方程P-V曲線。

(a)TNT

3.2 數(shù)值檢驗

為了進一步檢驗本文辨識方法的有效性和精度，利用顯式動力學軟件AUTODYN-2D對Φ25.4 mm標準圓筒試驗過程進行數(shù)值計算，考慮到圓筒模型的對稱性，建立二維X軸對稱模型。從一端起爆裝填在銅管內(nèi)的炸藥，起爆方式為面起爆，同時在距離起爆端200 mm圓筒外徑處設置一個高斯點記錄圓筒徑向位移時程曲線。圓筒試驗仿真涉及到的材料包括炸藥、空氣和銅管，其中銅管材料為OHFC無氧銅，采用Steinberg-Guinan強度模型和Mie-Grüneisen狀態(tài)方程，材料參數(shù)均取自AUTODYN材料庫，炸藥JWL狀態(tài)方程參數(shù)則分別采用辨識JWL參數(shù)和標準JWL參數(shù)。炸藥和空氣采用歐拉網(wǎng)格，空氣域大小為405.0 mm×76.2 mm，銅管采用拉格朗日網(wǎng)格，并在仿真中設置歐拉和拉格朗日網(wǎng)格互相耦合。為了保證計算精度的同時提高計算速度，銅管網(wǎng)格大小劃分成0.5 mm×0.5 mm，空氣和炸藥網(wǎng)格大小劃分為0.500 mm×0.488 mm。為了防止沖擊波在邊界反射，空氣域的邊界條件設置成無反射流出Flow-out。

圖7和圖8分別為圓筒外壁徑向位移時程曲線和圓筒外壁速度變化對比圖。從圖7可看出，采用辨識JWL參數(shù)和標準JWL參數(shù)計算得到的圓筒外壁徑向位移曲線整體偏差較小，尤其是圓筒壁膨脹前期，實線和虛線基本重合。TNT、HMX、Comp B和PBX-9404的圓筒徑向位移曲線的R2值分別為0.999 9、0.997 4、0.997 9和0.998 0，其值均接近1。圖8顯示對于所有炸藥的圓筒外壁速度變化曲線，實線在波形、起跳點和速度值等方面較好地重現(xiàn)了虛線，且R2值均大于0.99。兩組參數(shù)數(shù)值計算結果的高相似性證明了本文提出的利用遺傳算法和圓筒能量模型辨識炸藥JWL狀態(tài)方程參數(shù)的有效性和高精度。

圖7 圓筒外壁徑向位移時程曲線對比圖

圖8 圓筒外壁移動速度對比圖

4 結 論

本文提出了一種適用于計算密度大于1.0 g/cm3凝聚態(tài)CHNO型炸藥JWL狀態(tài)方程參數(shù)的方法，該方法基于基因遺傳算法和圓筒能量模型，相比于傳統(tǒng)的圓筒試驗法，只需已知炸藥的初始密度和爆速。主要結論如下：

(1)已知炸藥的初始密度和爆速，可計算出E0、γ、VCJ、pCJ、u6 mm和u19 mm的值，將這些參數(shù)代入自編程遺傳算法程序中可辨識得到爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)，此種方法經(jīng)濟、安全、方便、準確。

(2)基于本文所提出的方法，所得的4種炸藥的JWL狀態(tài)方程的P-V曲線與標準JWL狀態(tài)方程的P-V曲線相吻合，決定系數(shù)R2值均大于0.99。

(3)仿真結果顯示采用辨識JWL參數(shù)計算得到的4種炸藥的圓筒外壁徑向位移曲線決定系數(shù)R2值均大于0.99，證明了采用遺傳算法辨識炸藥JWL狀態(tài)方程參數(shù)的有效性和高精度，能夠滿足工程應用。