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

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

（1.北京應用物理與計算數學研究所，北京 100094；2.中國工程物理研究院研究生部，北京 100088）

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

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

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

1 計算方法

1.1 控制方程

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

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

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

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

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

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

1.2 數值計算方法

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

1.3 狀態方程

對炸藥反應物、產物均采用JWL狀態方程

式中：i為s或g，分別代表反應物和產物；Ai、Bi、、、ωi和為參數；V＝ρ0／ρ，E＝ρ0e，e為比內能，T為溫度。混合物狀態方程采用體積可加、內能可加、溫度平衡和壓力平衡等假設得到。若定義

由內能可加、等溫假設和狀態方程，可得

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

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

2 鈍感炸藥的爆轟數值模擬

2.1 計算參數

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

2.2 鈍感炸藥的沖擊點火

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

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

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

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

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

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

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

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

2.4 鈍感炸藥直徑效應

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

2.5 鈍感炸藥死區數值模擬

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

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

3 結 論

在彈塑性流體力學程序中引入點火增長的反應率模型和炸藥減敏模型，借助網格自適應加密技術，對鈍感炸藥的沖擊點火、炸藥直徑效應和帶拐角炸藥的死區形成等爆轟現象進行了數值模擬。數值結果表明，標準反應率模型能夠正確模擬平面爆轟波的特征狀態和鈍感炸藥的直徑效應，在引入減敏模型后能夠較好地模擬鈍感炸藥的死區形成過程。鈍感炸藥的數值研究表明，自主研發的彈塑性流體力學程序應用于炸藥研究是可行的。

［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，Forest 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］馮其京，何長江，張敏，等.二維高速碰撞問題歐拉數值模擬的混合網格計算［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，FENG 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.