苗飛超，周霖，張向榮，曹同堂

(1.北京理工大學 爆炸科學與技術(shù)國家重點實驗室, 北京 100081; 2.安徽東風機電科技股份有限公司, 安徽 合肥 230000)

0 引言

沖擊起爆特性是炸藥的重要性能之一，是作用可靠性和使用安全性評價的重要依據(jù)。Lee & Tarver點火增長反應(yīng)速率方程[1-2]已廣泛用于沖擊起爆的數(shù)值模擬,并成功重現(xiàn)了圓筒實驗[3]、拉氏分析實驗[4-6]、飛片沖擊實驗[7]、微通道裝藥爆轟實驗[8]的實驗結(jié)果。但點火增長反應(yīng)速率方程的參數(shù)較多(三項式的參數(shù)多達15個)且參數(shù)間存在關(guān)聯(lián)[9]，這導致點火增長反應(yīng)速率方程的標定存在較大的困難。為此，Murphy[9]對反應(yīng)速率方程的形式進行了修改，得到了參數(shù)更少，參數(shù)間無關(guān)聯(lián)的反應(yīng)速率方程。由于改進的反應(yīng)速率方程保留了Lee & Tarver模型中點火增長的概念，因此兩者在模擬沖擊起爆方面具有相同的功能。盡管改進的反應(yīng)速率方程具有以上優(yōu)點，但目前為止，大型商業(yè)軟件，包括LS-DYNA、AUTODYN和Abaqus等中并未包含該反應(yīng)速率方程。

因此，本文將改進的點火增長模型以自定義狀態(tài)方程的形式嵌入LS-DYNA軟件中，并采用此模型對DNAN基熔注炸藥的沖擊起爆過程進行了數(shù)值模擬，通過對比計算壓力曲線和實驗壓力曲線標定了改進的模型參數(shù)。

1 改進的點火增長反應(yīng)速率方程

1.1 模型簡介

三項式的Lee & Tarver反應(yīng)速率方程[2]的表達式為

(1)

式中：F為反應(yīng)分數(shù)；p為壓力；ρ和ρ0分別是當前時刻和初始時刻的密度；I、b、a、x、G1、c、d、y、G2、e、g、z、Figmax、FG1max和FG2min為15個可調(diào)參數(shù)；等號右側(cè)3項分別為點火項、慢速增長項和快速完成項。Murphy[9]保留Lee & Tarver反應(yīng)速率方程中點火增長的概念，對上述反應(yīng)速率方程進行了改進。改進的反應(yīng)速率方程表達式為

(2)

式中：η=ρ/ρ0-1；Freq、figmax、Ccrit、Eeta1、Grow1、Es1、em、Grow2、Es2、en為10個可調(diào)參數(shù)。與Lee & Tarver反應(yīng)速率方程相比，改進的點火- 增長反應(yīng)速率方程具有以下優(yōu)點：

1)慢速增長項和快速完成項的系數(shù)相關(guān)性弱。改進的點火增長反應(yīng)速率方程用sin(πFEs)代替了Lee & Tarver反應(yīng)速率方程慢速增長項和快速完成項的(1-F)mFn(m為c或e，n為d或g)。顯然，Es只影響sin(πFEs)曲線的形狀，峰值大小始終為1.0，而m和n的大小同時影響(1-F)mFn曲線的峰值和形狀。換言之，調(diào)整曲線的形狀時，Lee & Tarver模型需要同時改變2個參數(shù)(m和n)，而改進的模型只改變1個參數(shù)(Es)，簡化了參數(shù)的標定過程。

2) 反應(yīng)速率方程的參數(shù)較少。改進的反應(yīng)速率方程的參數(shù)有10個，比Lee & Tarver反應(yīng)速率方程少5個，進一步降低了參數(shù)標定的難度。

1.2 模型嵌入LS-DYNA軟件

鑒于改進的點火- 增長反應(yīng)速率方程具有上述優(yōu)點，本文將改進的點火- 增長反應(yīng)速率方程以自定義狀態(tài)方程的形式嵌入LS-DYNA軟件中。下面對自定義狀態(tài)方程計算中，未反應(yīng)炸藥和爆轟產(chǎn)物的混合法則、狀態(tài)方程形式和各物理量的求解過程進行介紹。其中狀態(tài)方程混合采用Cochran等[10]的算法。

部分反應(yīng)炸藥是未反應(yīng)炸藥和爆轟氣體產(chǎn)物的混合物，假設(shè)兩種組分處于壓力p和溫度T平衡狀態(tài)，即

(3)

式中：下標m代表部分反應(yīng)炸藥；下標s代表未反應(yīng)炸藥；下標g代表爆轟氣體產(chǎn)物。部分反應(yīng)炸藥的相對體積Vm和內(nèi)能Em，滿足加和特性，

(4)

式中：相對體積V=v/v0，v和v0分別是當前時刻和初始時刻的比容；內(nèi)能E=e/e0，e和e0分別是當前時刻和初始時刻的比內(nèi)能。

未反應(yīng)炸藥和爆轟產(chǎn)物均采用含溫度的JWL狀態(tài)方程[2]，其形式為

(5)

式中：*為s或g，分別代表未反應(yīng)炸藥和爆轟產(chǎn)物；A、B、R1、R2、ω、CV均為可調(diào)參數(shù)。

溫度、壓力和反應(yīng)分數(shù)等物理量的求解過程如下：

1)溫度的計算。

為簡化計算，將能量方程轉(zhuǎn)化為溫度方程。溫度方程的表達式為

(6)

式中：Sij為偏應(yīng)力張量；εij為應(yīng)變張量；q為人工黏性項；CVm、J以及H分別為

(7)

采用預測- 校正格式對(6)式進行顯示差分計算，可得到ΔTm.

2)真實體積分數(shù)和壓力的計算。

定義β=(1-F)Vs/Vm為未反應(yīng)炸藥的真實體積分數(shù)[10]，則Vs和Vg可表示為

Vs=βVm/(1-F)，

(8)

Vg=(1-β)Vm/F.

(9)

由于未反應(yīng)炸藥和爆轟氣體產(chǎn)物處于壓力平衡，因此

pg-ps=δ(Tm,Vm,F,β)=0.

(10)

在給定Tm、Vm、F的條件下，采用牛頓迭代法對(10)式進行迭代求解，迭代變量為β. 在迭代過程中，通過(8)式和(9)式計算Vs和Vg，并將Vs和Vg代入單組分狀態(tài)方程中計算ps和pg.β通過迭代收斂時，得到pm=0.5(ps+pg)。

3)反應(yīng)分數(shù)F的計算。

反應(yīng)速率方程可表示為

(11)

式中：R為任意函數(shù)。本文采用Cochran等[10]的差分格式計算該方程，該格式兼具顯示格式和隱式格式的優(yōu)點，其形式為

(12)

式中：上標n代表第n個時間步；Δt為時間步長；

(13)

(13)式中λ的取值方法可保證反應(yīng)速率方程的差分求解過程是完全穩(wěn)定的。各物理量在一個時間步長內(nèi)的迭代求解過程如圖1所示。

圖1 各物理量迭代求解流程圖Fig.1 Flow chart of iteration for updating the variables

2 沖擊起爆數(shù)值模擬

圖2 一維拉格朗日分析測試系統(tǒng)示意圖[5]Fig.2 Schematic diagram of one-dimensional Lagrange analysis and measurement system[5]

Cao等[11]采用如圖2所示的一維拉格朗日分析測試系統(tǒng)，測量了DNAN基熔注炸藥(DNAN 29.5%/RDX 40%/Al 30%/N-甲基-4-硝基苯胺(MNA)0.5%)沖擊起爆過程中的壓力曲線。測試時，由透鏡產(chǎn)生的平面波經(jīng)空氣隙和鋁板衰減后入射到炸藥中。待測炸藥由3個厚度3 mm藥片和1個厚度20 mm的藥片組成，藥片直徑均為50 mm. 4個錳銅傳感器分別安裝在距離鋁衰減板0 mm(h1)、3 mm(h2)、6 mm(h3)和9 mm(h4)位置處。最終實驗得到圖2中4個傳感器記錄的壓力曲線。本文以該壓力曲線為基礎(chǔ)，標定改進的點火增長模型參數(shù)。通過對比計算壓力曲線和實驗壓力曲線，驗證改進的點火增長模型能否正確模擬DNAN基熔注炸藥的沖擊起爆過程。

2.1 幾何模型

為節(jié)省計算內(nèi)存和計算時間，只建立了待測炸藥的幾何模型，并將實驗測得的圖2中衰減板與炸藥界面處(即圖2中h1)的壓力變化歷史作為炸藥的壓力邊界條件。圖3中4個觀測點位置對應(yīng)圖2中4個錳銅傳感器位置，即h1(0 mm)、h2(3 mm)、h3(6 mm)和h4(9 mm)，以便將計算值與實驗值對比。

圖3 數(shù)值模擬幾何模型Fig.3 Geometric model of numerical simulations

2.2 材料參數(shù)

未反應(yīng)炸藥和爆轟產(chǎn)物狀態(tài)方程參數(shù)來自文獻[11]。但文獻[11]中未反應(yīng)炸藥采用的是Gruneisen狀態(tài)方程。考慮到本文采用的是JWL狀態(tài)方程，因此需要將Gruneisen狀態(tài)方程擬合為JWL狀態(tài)方程。本文采用遺傳算法對其進行擬合，具體擬合方法見文獻[12]。擬合結(jié)果如圖4所示，擬合參數(shù)如表1所示。

圖4 JWL狀態(tài)方程參數(shù)擬合Fig.4 Fitting result of JWL equation of state

反應(yīng)速率方程參數(shù)通過試錯法優(yōu)化得到。通過不斷調(diào)整反應(yīng)速率方程參數(shù)，減小計算值相對實驗值的偏差。最終，優(yōu)化得到的反應(yīng)速率方程參數(shù)如表2所示。

2.3 計算結(jié)果與分析

采用表1和表2中的參數(shù)對圖3的模型計算，得到的計算值和實驗值如圖5所示。從圖5中可以看出，沖擊波到達時間和壓力駝峰對應(yīng)時間重合度較高，計算值與實驗值吻合較好。但是，h3和h4處計算的壓力峰值與實驗值相差較大。這可能是封裝傳感器的聚四氟乙烯與炸藥的阻抗不匹配導致的。因此，在對計算結(jié)果進行評價時，應(yīng)該更多的關(guān)注沖擊波到達時間、壓力變化趨勢和駝峰出現(xiàn)位置與實驗值是否吻合。從圖5來看，改進的點火增長模型能夠較好地描述炸藥的沖擊起爆過程。

表1 JWL狀態(tài)方程參數(shù)

注：DCJ為爆速，pCJ為爆壓，E0為內(nèi)能。

表2 改進的點火增長反應(yīng)速率方程參數(shù)

圖5 4個觀測點處的壓力時程曲線實驗值與計算值對比Fig.5 Comparison of experimental and simulated values of pressures at four fixed points

相對于軟件自帶模型，自定義模型可以輸出更多的計算信息，有助于判斷模型嵌入LS-DYNA軟件時采用的壓力更新算法是否可行。壓力更新算法中迭代變量的選擇尤為重要。本文2.2節(jié)中，采用牛頓法對(10)式進行迭代求解時取β為迭代變量，并沒有取Vs或Vg作為迭代變量，具體原因結(jié)合本次數(shù)值模擬計算結(jié)果分析如下。

圖6為本文數(shù)值模擬計算得到的4個觀測點處未反應(yīng)炸藥相對體積的變化曲線。圖6中h2位置處，初始時刻Vs為1，到計算結(jié)束時約為33，變化范圍較大。因此，如果以Vs作為迭代變量，則迭代次數(shù)較多，計算效率低。此外，在傳感器有效記錄時間內(nèi)，如果入射到炸藥中的沖擊波持續(xù)時間較短，則稀疏波的作用時間會比較長，Vs可能達到103量級。此時，在對Vs進行迭代時，容易出現(xiàn)不收斂情況，并且由于計算機存在舍入誤差，在計算時甚至會出現(xiàn)0/0型的錯誤。

圖6 4個觀測點處未反應(yīng)炸藥相對體積Fig.6 Relative volumes of unreacted explosives at four fixed points

圖7為本文數(shù)值模擬計算得到的4個觀測點處爆轟產(chǎn)物相對體積的變化曲線。由圖7可知，Vg的變化范圍較小，將其作為迭代變量似乎可行。但是，在波陣面處，Vg的變化極為劇烈，類似于跳躍間斷點。因此，若將Vg作為迭代變量，在波陣面處進行迭代時很難收斂。

圖7 4個觀測點處爆轟產(chǎn)物相對體積Fig.7 Relative volumes of detonation products at four fixed points

圖8為本文數(shù)值模擬計算得到的4個觀測點處迭代變量β的變化曲線。由圖8可知，在整個計算過程中β始終在0～1的范圍內(nèi)變化，β變化范圍較小。此外，整個計算過程中β光滑過渡，有助于快速收斂。本文計算中，每個網(wǎng)格的迭代次數(shù)不超過3次，證明本文所采用迭代算法是高效、可行的。因此，采用β作為迭代變量有利于減少迭代次數(shù)，提高計算效率。

圖8 4個觀測點處迭代變量βFig.8 Iterative variables β at four fixed points

圖9 4個觀測點處的壓力與相對體積路徑Fig.9 p-V paths at four fixed points

此外，自定義模型可以輸出更多的流場信息，有助于進一步觀察反應(yīng)區(qū)結(jié)構(gòu)，并對實驗設(shè)計提供指導。為了更詳細地描述沖擊起爆過程，可以對圖9中4個觀測點在p-V圖上的路徑與雨貢紐曲線和等熵線的關(guān)系進行分析。以h1處觀測點為例，沖擊波入射到炸藥中，壓力從0 GPa升高到7.5 GPa，落在雨貢紐曲線；隨后，炸藥開始反應(yīng)，反應(yīng)分數(shù)逐漸增加，(V，p)點逐漸向等熵線靠近；當相對體積為1.0(圖10中“+”)時，對應(yīng)的反應(yīng)分數(shù)為0.64，隨后繼續(xù)膨脹。其他3條曲線的過程類似。值得注意的是，h4處觀測點反應(yīng)分數(shù)達到1時，相對體積為0.95. 而Chapman-Jouguet點對應(yīng)的相對體積為0.72，所以在h4處仍未建立穩(wěn)定爆轟。因此，如果想在實驗中記錄到穩(wěn)定爆轟的狀態(tài)，可以適當增加最后一個傳感器位置到邊界的距離。

圖10 4個觀測點處反應(yīng)分數(shù)Fig.10 Reaction fractions at four fixed points

3 結(jié)論

本文通過自定義狀態(tài)方程的形式將改進的點火增長模型嵌入LS-DYNA軟件，對DNAN基含鋁熔注炸藥沖擊起爆進行了數(shù)值模擬，并與實驗值進行了對比。得出以下結(jié)論：

1)計算得到的4個傳感器位置處的沖擊波到達時間與實驗值相差不超過0.2 μs，改進的點火增長模型能夠很好地描述炸藥的沖擊起爆過程。

2)以β作為迭代變量的壓力更新算法計算效率高。本文計算中，每個網(wǎng)格的迭代次數(shù)不超過3次。

3)將Lagrange位置處的壓力與相對體積路徑線與未反應(yīng)炸藥的雨貢紐曲線、爆轟產(chǎn)物的等熵線畫在一張圖上有助于理解炸藥在沖擊起爆過程中的狀態(tài)變化、判斷該Lagrange位置處的炸藥是否達到穩(wěn)定爆轟和深化對沖擊起爆過程的認識。