鈍感炸藥點火增長模型的歐拉數(shù)值模擬＊

郝鵬程，馮其京，洪 滔，王言金

（1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100094；2.中國工程物理研究院研究生部，北京 100088）

塑性炸藥的爆轟機理很復(fù)雜，對爆轟波結(jié)構(gòu)的研究，早在一個多世紀(jì)以前就提出了平面、定常的CJ爆轟模型，即認(rèn)為爆轟波由前導(dǎo)沖擊波及緊跟的稀疏波構(gòu)成，能量在沖擊波后瞬間釋放；ZND模型則在前導(dǎo)沖擊波后引入化學(xué)反應(yīng)區(qū)［1］。實際上，炸藥內(nèi)部微觀結(jié)構(gòu)復(fù)雜，由大量晶粒以及塑性粘結(jié)劑、缺陷、氣孔等構(gòu)成，在化學(xué)反應(yīng)區(qū)內(nèi)存在復(fù)雜的作用機制?，F(xiàn)今大多認(rèn)為，當(dāng)炸藥受到強烈沖擊，由于摩擦、空洞塌陷、剪切帶的形成和局部塑性變形等，使炸藥局部急劇升溫，即所謂的熱點形成，消耗部分炸藥，進而引發(fā)周圍炸藥點火起爆連成一片，在化學(xué)反應(yīng)區(qū)內(nèi)炸藥逐漸轉(zhuǎn)化為產(chǎn)物。在這一物理圖景下，提出了許多唯象反應(yīng)率模型，如JTF模型［2］、HVRB模型［3］、點火增長模型［4－5］和統(tǒng)計熱點模型［6］等。在這些模型中，引入產(chǎn)物質(zhì)量分額的概念，反應(yīng)區(qū)內(nèi)的炸藥為反應(yīng)物與產(chǎn)物的混合物，在一定假設(shè)之下得到混合物的狀態(tài)。點火增長模型被廣泛應(yīng)用于模擬各種爆轟以及炸藥與惰性材料的相互作用等現(xiàn)象。

唯象反應(yīng)率模型從較簡單的物理特征出發(fā)，常被推廣應(yīng)用于復(fù)雜的爆轟狀況下，如炸藥直徑效應(yīng)的模擬，也能得到令人滿意的結(jié)果。凝聚炸藥存在減敏現(xiàn)象，即受到弱沖擊的鈍感炸藥與未受沖擊的炸藥相比，起爆炸藥需要更強的沖擊壓力，并且，這一效應(yīng)是與時間相關(guān)的過程［7］。當(dāng)爆轟波遇到球形發(fā)散或拐角時，出現(xiàn)全部或部分炸藥未被起爆的區(qū)域［8］，形成死區(qū)。單純的反應(yīng)率模型不能很好地模擬這一現(xiàn)象，特別是當(dāng)爆轟波達到拐角后，由于稀疏效應(yīng)壓力減弱，無法起爆拐角后的炸藥，但當(dāng)爆轟波繞過拐角，壓力增強會引爆先前未爆的炸藥，先前形成的死區(qū)逐漸消失；引入減敏模型［9］后，模擬結(jié)果與實驗現(xiàn)象接近。

自主研發(fā)的多介質(zhì)歐拉彈塑性流體力學(xué)程序（Multicomponent Eulerian elastic－plastic hydrodynamics code，MEPH2Y）能有效地模擬流體、彈塑性介質(zhì)高速碰撞等問題［10］，并行網(wǎng)格自適應(yīng)加密（A－daptive mesh refinement，AMR）技術(shù)有效地擴展了MEPH2Y的適用范圍［11］。為了研究鈍感炸藥的各種爆轟現(xiàn)象，本文中擬在MEPH2Y中加入反應(yīng)進程變量和減敏進程變量的控制方程，通過對鈍感炸藥的平面爆轟波傳播、直徑效應(yīng)和死區(qū)形成等爆轟現(xiàn)象的數(shù)值模擬，驗證該程序用于研究鈍感炸藥的適用性和可靠性。

1 計算方法

1.1 控制方程

為了模擬鈍感炸藥的各種爆轟現(xiàn)象，在彈塑性流體力學(xué)控制方程中，引入炸藥反應(yīng)率和減敏模型，即加入反應(yīng)進程變量F和減敏進程變量ψ的控制方程

式中：t為時間；SF和Sψ分別為與F和ψ對應(yīng)的源項，F(xiàn)和ψ在0～1之間；F＝1，表示炸藥均為產(chǎn)物；ψ＝1，表示炸藥敏感性最低。式（1）右端的反應(yīng)率源項SF采用三項點火增長反應(yīng)率模型［5］

式中：H（x）為階躍函數(shù)，I、G1、G2、a、b、c、d、e、g、x、y、z、FIG，max、FG1，max和FG2，min為炸藥參數(shù)，ρ0和ρ分別為初始密度和密度，p為壓力。RIG、和分別為源項SF的點火項和2個增長項。在炸藥沖擊點火和爆轟波傳播等不同的物理過程中，RIG、和等3項的意義有所差別［5］。

式（2）定義的減敏進程變量ψ通過如下方式進入反應(yīng)進程變量的控制方程

式（4）中的a為式（3）中點火項RIG的壓縮率的閾值，a0、a1和Fc為參數(shù)。式（2）右端的減敏函數(shù)的源項可以采用不同的形式，本文中選取如下表達式［9］

式中：A和ε為參數(shù)。則定常壓力p下，若減敏進程變量ψ從0積分到ψ1，需時t，積分式（7），可得

1.2 數(shù)值計算方法

MEPH2Y采用交替方向方法解無源項的方程組，每個方向為拉氏步計算緊跟輸運步計算，當(dāng)2個方向的無源項方程組解完后計算源項。AMR技術(shù)［11］在邏輯上采用嵌套的歐拉型網(wǎng)格塊，網(wǎng)格塊細(xì)分采用在每個方向上一分為二的方式，不同加密級的網(wǎng)格塊可選用相等或不等的時間步長，僅計算受到擾動的計算塊。在二維爆轟問題計算中，跟蹤加密反應(yīng)區(qū)的計算網(wǎng)格。

1.3 狀態(tài)方程

對炸藥反應(yīng)物、產(chǎn)物均采用JWL狀態(tài)方程

式中：i為s或g，分別代表反應(yīng)物和產(chǎn)物；Ai、Bi、、、ωi和為參數(shù)；V＝ρ0／ρ，E＝ρ0e，e為比內(nèi)能，T為溫度?；旌衔餇顟B(tài)方程采用體積可加、內(nèi)能可加、溫度平衡和壓力平衡等假設(shè)得到。若定義

由內(nèi)能可加、等溫假設(shè)和狀態(tài)方程，可得

由等壓假設(shè)、式（11）和狀態(tài)方程，并考慮到體積可加假設(shè)，得到非線性方程

解上述非線性方程，可采用Newton迭代法或二分法。模擬死區(qū)時，在稀疏區(qū)內(nèi)對初值的選取要求較嚴(yán)格，Newton迭代往往不收斂到正確解，此時采用二分法較合適。

2 鈍感炸藥的爆轟數(shù)值模擬

2.1 計算參數(shù)

算例涉及鈍感炸藥TATB、LX－17和PBX9502。炸藥反應(yīng)物及產(chǎn)物JWL狀態(tài)方程和點火增長模型的參數(shù)分別采用文獻［4，9，5］中相應(yīng)炸藥的參數(shù)。炸藥LX－17減敏模型的相關(guān)參數(shù)分別為：A＝8.364 2s－1·kPa－1，ε＝0.001，a0＝0.22，a1＝0.5，F(xiàn)c＝0.01；其中，A通過將p＝1GPa，ψ1＝0.98，t＝1.29μs代入式（8）得到［9］。

2.2 鈍感炸藥的沖擊點火

以11.38GPa的壓力持續(xù)沖擊TATB，采用1／100mm網(wǎng)格計算。表1為穩(wěn)定爆轟波爆速、von Neumann尖點和CJ狀態(tài)，計算值與實驗值［4］一致。N和CJ分別表示尖點和CJ狀態(tài)。

表1 TATB炸藥爆轟波的von Neumann尖點和CJ狀態(tài)參數(shù)Table 1Parameters for detonation von Neumann spike and CJ states of TATB explosive

2.3 引入減敏模型的鈍感炸藥沖擊點火

以12.55和14.97GPa的壓力持續(xù)沖擊LX－17炸藥，計算中采用1／200mm網(wǎng)格。

圖1為不同時刻的壓力曲線，從0時刻開始，從左向右依次間隔Δt＝0.1μs，實線為標(biāo)準(zhǔn)模型結(jié)果，而虛線為減敏模型結(jié)果。采用標(biāo)準(zhǔn)點火增長模型模擬，均能形成自持爆轟；而采用減敏模型模擬，較弱的沖擊則不能起爆炸藥，這與文獻［9］結(jié)果一致。

圖1 LX－17炸藥在不同持續(xù)沖擊壓力下的壓力曲線Fig.1Pressure curves of LX－17explosive under different sustaining impact pressures

圖2為某時刻不同沖擊壓力下反應(yīng)進程變量和閾值的空間分布，引入減敏函數(shù)后，在較弱的沖擊壓力下（見圖2（a）），當(dāng)ρ／ρ0＜1＋a時，反應(yīng)率中點火項被抑制；F＜時，增長項的發(fā)展也被抑制，使炸藥無法起爆。而較強的沖擊壓力下增長項未被抑制（見圖2（b）），如圖1（b）所示，壓力曲線無太大影響，爆轟能夠正常形成。

圖2 用減敏模型得到的LX－17炸藥的壓縮率和反應(yīng)率進程變量的空間分布Fig.2Compression rates and reactive progress variables of LX－17explosive using the desensitizing model

2.4 鈍感炸藥直徑效應(yīng)

計算中采用軸對稱幾何，采用AMR技術(shù)跟蹤反應(yīng)區(qū)，反應(yīng)區(qū)內(nèi)網(wǎng)格寬度Δx為1／10、1／40mm，初始時刻在柱形 PBX9502炸藥的前端設(shè)置31.5GPa的高壓。圖3中給出了藥柱內(nèi)爆轟波傳播速度與藥柱曲率的關(guān)系，實心點線為實驗結(jié)果［12－13］，空心點線為計算結(jié)果。從圖中可以看到計算結(jié)果有收斂到實驗值的趨勢。另外，文獻［12］中報道在爆速低于7.3km／s時爆轟波將逐漸熄滅。本文中也觀察到：較小藥柱在較高初始壓力（如31.5GPa）下，能夠持續(xù)較長時間的爆轟波傳播；但在較低壓力（如14.0GPa）下則會熄滅。

2.5 鈍感炸藥死區(qū)數(shù)值模擬

研究帶凸井的鈍感炸藥LX－17，模擬爆轟波在經(jīng)過拐角后失敗，即存在死區(qū)的現(xiàn)象。采用文獻［9］中模型的幾何尺寸，如圖4所示，使用柱坐標(biāo)系。拐角兩側(cè)壁面和對稱軸為反射邊界條件，其他為出口邊界條件。初始時刻，在凸井中心R＝7.78mm的半球內(nèi)，取為爆轟產(chǎn)物（ρ＝ρ0，p＝41.46GPa，F(xiàn)＝1，ψ＝1），之外為反應(yīng)物（ρ＝ρ0，p＝0，F(xiàn)＝ψ＝0）。采用AMR技術(shù)計算，跟蹤反應(yīng)區(qū)，反應(yīng)區(qū)內(nèi)采用1／40mm網(wǎng)格。

圖3 PBX9502炸藥的爆速－曲率關(guān)系Fig.3Detonation speed－curvature curves of PBX9502explosive

圖4 炸藥幾何形狀Fig.4Geometric shape of explosives

圖5 LX－17炸藥帶拐角的爆轟Fig.5Detonation with the corner of LX－17explosive

圖5為數(shù)值模擬結(jié)果，其中，上排為標(biāo)準(zhǔn)反應(yīng)率模型計算結(jié)果，而下排為帶減敏反應(yīng)率模型的計算結(jié)果；每張圖的上半部分為壓力等值云圖，下半部分為反應(yīng)進程變量等值線。t＜1.6μs（見圖5（a）），球面爆轟波逐漸形成并呈球形向外傳播，2模型計算結(jié)果無明顯差別。1.6μs＜t＜2.6μs（見圖5（b）），沖擊波在拐角處轉(zhuǎn)彎稀疏并開始減弱，壓力降低，在拐角附近開始彎曲，在受稀疏的前沿，爆轟波無法形成，而在遠離拐角的炸藥仍然維持柱面爆轟，2模型結(jié)果無明顯差別。而后如圖5（c）～（f）所示，沖擊波前沿嚴(yán)重彎曲，在標(biāo)準(zhǔn)模型中，先前受弱沖擊的一部分炸藥在二次沖擊壓力下，點火并形成爆轟，在t＝3.8μs彎曲的沖擊波與凸井碰撞，此后如圖5（g）中標(biāo)準(zhǔn)模型計算結(jié)果所示，大部分炸藥被彎曲的沖擊波沖擊起爆，僅在拐角處有極少量未起爆炸藥，即死區(qū)僅局部地、暫時地存在；然而，在減敏模型計算中，受弱沖擊的炸藥敏感性降低，彎曲的較弱的沖擊無法使炸藥起爆形成爆轟，最終在拐角附近形成死區(qū)，直到計算結(jié)束，死區(qū)仍然存在。值得注意的是，在標(biāo)準(zhǔn)模型計算中，部分受二次沖擊的炸藥雖然起爆，但該區(qū)域的溫度始終僅500K左右，對于爆轟來講不大合理，這也說明減敏模型有一定的合理性，也許反應(yīng)率函數(shù)應(yīng)該考慮更多的因素。另外，平面幾何與軸對稱幾何的結(jié)果相似，標(biāo)準(zhǔn)反應(yīng)率模型亦僅形成局部、暫時的死區(qū)而后所剩無幾；而在考慮了減敏效應(yīng)的模型中死區(qū)能夠永久存在，如圖5（h）所示，只是死區(qū)形狀略有不同。

3 結(jié) 論

在彈塑性流體力學(xué)程序中引入點火增長的反應(yīng)率模型和炸藥減敏模型，借助網(wǎng)格自適應(yīng)加密技術(shù)，對鈍感炸藥的沖擊點火、炸藥直徑效應(yīng)和帶拐角炸藥的死區(qū)形成等爆轟現(xiàn)象進行了數(shù)值模擬。數(shù)值結(jié)果表明，標(biāo)準(zhǔn)反應(yīng)率模型能夠正確模擬平面爆轟波的特征狀態(tài)和鈍感炸藥的直徑效應(yīng)，在引入減敏模型后能夠較好地模擬鈍感炸藥的死區(qū)形成過程。鈍感炸藥的數(shù)值研究表明，自主研發(fā)的彈塑性流體力學(xué)程序應(yīng)用于炸藥研究是可行的。

［1］Fickett W，Davis W C.Detonation：Theory and experiment［M］.Berkeley：University of California Press，1979：13－54.

［2］Johnson J N，Tang P K，F(xiàn)orest C A.Shock wave initiation of heterogeneous reactive solids［J］.Journal of Applied Physics，1985，57（9）：4323－4334.

［3］Starkenberg J.Modeling detonation propagation and failure using explosive initiation models in a conventional hydrocode［C］∥Wardell J F，Maienschein J L.Proceedings of 12th International Detonation Symposium.San Diego，CA：Office of Naval Research，2002：1381－1390.

［4］Lee E L，Tarver C M.Phenomenological model of shock initiation in heterogeneous explosives［J］.Physics of Fluids，1980，23（12）：2362－2372.

［5］Tarver C M，McGuire E M.Reactive flow modeling of the interaction of TATB detonation waves with inert materials［C］∥Wardell J F，Maienschein J L.Proceedings of 12th International Detonation Symposium.San Diego，CA：Office of Naval Research，2002：971－980.

［6］Nichols III A L，Tarver C M.A statistical hot spot reactive flow model for shock initiation and detonation of solid high explosives［C］∥Wardell J F，Maienschein J L.Proceedings of 12th International Detonation Symposium.San Diego，CA：Office of Naval Research，2002：779－788.

［7］Campbell A W，Travis J R.The shock desensitization of PBX－9404and composition B3［C］∥Davis W C.Proceedings of the 8th International Detonation Symposium.Albuquerque，NM：Naval Surface Warfare Center，1985：1057－1068.

［8］Cox M，Campbell A W.Corner－turning in TATB［C］∥Short J M.Proceedings of the 7th International Detonation Symposium.White Oak，MD：Naval Weapons Center，1981：624－633.

［9］DeOliveira G，Kapila A K，Schwendeman D W，et al.Detonation diffraction dead zones and the ignition－and－growth model［C］∥Peiris S M.Proceedings of the 13th International Detonation Symposium.Norfolk，VA：Office of Na－val Research，2006：13－22.

［10］馮其京，何長江，張敏，等.二維高速碰撞問題歐拉數(shù)值模擬的混合網(wǎng)格計算［J］.計算物理，2003，20（2）：147－152.

FENG Qi－jing，HE Chang－jiang，ZHANG Min，et al.Calculation of mixed cells in simulating high－speed collision problems with two－dimensional Eulerian hydrodynamic method［J］.Chinese Journal of Computational Physics，2003，20（2）：147－152.

［11］HAO Peng－cheng，F(xiàn)ENG Qi－jing，YAO Wen.Two－dimensional Euler adaptive mesh method on detonation［J］.Journal of Beijing Institute of Technology，2009，18（2）：141－145.

［12］Campbell A W.Diameter effect and failure diameter of a TATB－based explosive［J］.Propellants，Explosives，Pyrotechnics，1984，9（6）：183－187.

［13］Hill L G，Bdzil J B，Aslam T D.Front curvature rate stick measurements and detonation shock dynamics calibration for PBX 9502over a wide temperature range［C］∥Davis W C.Proceedings of the 11th International Detonation Symposium.Snowmass，CO：Office of Naval Research，1998：1029－1037.