點(diǎn)火增長(zhǎng)反應(yīng)速率方程在LS-DYNA軟件中嵌入及應(yīng)用

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

(1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100081; 2.安徽東風(fēng)機(jī)電科技股份有限公司, 安徽 合肥 230000)

0 引言

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

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

1 改進(jìn)的點(diǎn)火增長(zhǎng)反應(yīng)速率方程

1.1 模型簡(jiǎn)介

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

(1)

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

(2)

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

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

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

1.2 模型嵌入LS-DYNA軟件

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

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

(3)

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

(4)

式中：相對(duì)體積V=v/v0，v和v0分別是當(dāng)前時(shí)刻和初始時(shí)刻的比容；內(nèi)能E=e/e0，e和e0分別是當(dāng)前時(shí)刻和初始時(shí)刻的比內(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)分?jǐn)?shù)等物理量的求解過(guò)程如下：

1)溫度的計(jì)算。

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

(6)

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

(7)

采用預(yù)測(cè)- 校正格式對(duì)(6)式進(jìn)行顯示差分計(jì)算，可得到ΔTm.

2)真實(shí)體積分?jǐn)?shù)和壓力的計(jì)算。

定義β=(1-F)Vs/Vm為未反應(yīng)炸藥的真實(shí)體積分?jǐn)?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的條件下，采用牛頓迭代法對(duì)(10)式進(jìn)行迭代求解，迭代變量為β. 在迭代過(guò)程中，通過(guò)(8)式和(9)式計(jì)算Vs和Vg，并將Vs和Vg代入單組分狀態(tài)方程中計(jì)算ps和pg.β通過(guò)迭代收斂時(shí)，得到pm=0.5(ps+pg)。

3)反應(yīng)分?jǐn)?shù)F的計(jì)算。

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

(11)

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

(12)

式中：上標(biāo)n代表第n個(gè)時(shí)間步；Δt為時(shí)間步長(zhǎng)；

(13)

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

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

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

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

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

2.1 幾何模型

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

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

2.2 材料參數(shù)

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

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

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

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

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

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

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

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

圖5 4個(gè)觀測(cè)點(diǎn)處的壓力時(shí)程曲線實(shí)驗(yàn)值與計(jì)算值對(duì)比Fig.5 Comparison of experimental and simulated values of pressures at four fixed points

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

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

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

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

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

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

圖8 4個(gè)觀測(cè)點(diǎn)處迭代變量βFig.8 Iterative variables β at four fixed points

圖9 4個(gè)觀測(cè)點(diǎn)處的壓力與相對(duì)體積路徑Fig.9 p-V paths at four fixed points

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

圖10 4個(gè)觀測(cè)點(diǎn)處反應(yīng)分?jǐn)?shù)Fig.10 Reaction fractions at four fixed points

3 結(jié)論

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

1)計(jì)算得到的4個(gè)傳感器位置處的沖擊波到達(dá)時(shí)間與實(shí)驗(yàn)值相差不超過(guò)0.2 μs，改進(jìn)的點(diǎn)火增長(zhǎng)模型能夠很好地描述炸藥的沖擊起爆過(guò)程。

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

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