一種高精度歐拉算法在成形裝藥金屬射流問題中的應用＊

陳千一

（中國商用飛機有限責任公司上海飛機設計研究院，上海201210）

空心炸藥能將炸藥爆炸產(chǎn)生的能量匯聚到一起，在局部產(chǎn)生極強的作用力；如果給空心炸藥的鏤空部分加上襯墊，炸藥爆炸產(chǎn)生的高壓氣體可以將襯墊壓潰，形成一股高速射流，這就是成形裝藥。

成形裝藥金屬射流屬于爆炸加載問題，涉及到炸藥的爆轟過程、不同相態(tài)物質(zhì)之間的相互作用以及固體的大變形甚至破碎等高度非線性問題，目前主要靠實驗手段研究此類問題，數(shù)值模擬暫時只能起到輔助作用。但是實驗成本昂貴，又具有一定的危險性，而數(shù)值模擬成本低廉、高效安全，且可以獲取一些實驗技術(shù)無法獲取的信息。近年來，一些研究人員嘗試用現(xiàn)代計算流體力學的高精度格式來模擬爆炸問題。劉凱欣課題組將一種改進的CE／SE格式與粒子水平集方法相結(jié)合，成功模擬了爆炸焊接［1］和金屬射流［2］等爆炸加載問題，計算結(jié)果與實驗及他人的計算結(jié)果吻合得較好，證實了計算方法的可靠性。

材料在高速沖擊載荷作用下的溫度計算一直是一個棘手的問題，溫度對材料的動力學行為有著重要的影響，從塑性流動到結(jié)構(gòu)性相變、甚至液化或者汽化，其中都有溫度的作用。大多數(shù)文獻只考慮由塑性變形引起的溫度變化，而忽略了沖擊溫升［3－6］。對于沖擊溫升的計算已有一些教科書給出了詳細的討論［7－8］，傳統(tǒng)的沖擊溫升計算方法在歐拉型數(shù)值程序中會遇到2方面的問題：首先是程序中難以判斷加載過程與卸載過程，且在沖擊動力學問題中會出現(xiàn)反復的加載與卸載的循環(huán)過程，使得實際操作起來更加困難；其次，由于歐拉程序記錄的是空間點的物理量，而不是跟蹤某一物質(zhì)點的始末，因此無法通過狀態(tài)量直接求出溫度的變化。

本文中，構(gòu)建一套高精度的求解成形裝藥金屬射流問題的數(shù)值方案，提出一種適用于歐拉程序的同時考慮塑性溫升和沖擊溫升的計算方法。首先使用固體炸藥爆炸的經(jīng)典算例對自行開發(fā)的數(shù)值程序進行驗證，然后對成形裝藥金屬射流問題進行模擬，給出炸藥的起爆過程、金屬襯墊的失效以及射流的形成過程。討論不同襯墊形狀對射流的形狀、速度和溫度的影響。

1 算法描述

1.1 控制方程

使用流體彈塑性模型求解金屬射流，一般將應力張量拆開表示，

式中：σij為應力張量的分量，sij為偏應力張量的分量，p為靜水壓。當固體物質(zhì)受到強沖擊載荷時，其力學行為接近于流體，此時可以忽略應力偏量。

將軸對稱問題簡化到柱坐標（r，z，θ）的r－z平面上，不考慮外力、熱傳導和熱擴散，方程組的守恒形式可以寫成

式中：Q為守恒量矢量，E和F分別為r和z方向上的流通量矢量，S為源項矢量，它們的表達式為

式中：ρ為物質(zhì)的密度，u和v分別為r和z方向上的速度分量，E＝ρe＋（1／2）ρ（u2＋v2）為物質(zhì)的體積能量，e為物質(zhì)的質(zhì)量內(nèi)能，σrr、σzz、σθθ和σrz為應力張量的分量，srr、szz、sθθ和srz為偏應力張量的分量，εp為等效塑性應變，T為物質(zhì)的溫度。對于固體炸藥的氣體產(chǎn)物，只需要求解前4個方程。襯墊材料失效后，開始塑性流動并形成射流。隨著溫度的上升，應力偏量趨近于0，金屬材料表現(xiàn)得更像是無黏流體。使用改進的時－空守恒元解元（CE／SE）格式［9－10］求解上述控制方程組。

1.2 溫度的計算

考慮塑性溫升和沖擊溫升對射流溫度的影響。塑性溫升的計算公式為

由于金屬襯墊為固體材料，而固體材料的沖擊絕熱線與卸載等熵線十分接近［7］，假設加載過程和卸載過程都是等熵的，因此有

式中：dS為材料的熵增，cV為材料的比定容熱容。由γ／V＝γ0／V0以及dV＝－1／ρ（d）2ρ，可以得到

這種近似的方法比較適合于歐拉型的數(shù)值程序，計算時不需要判斷加載還是卸載。將塑性溫升和沖擊溫升結(jié)合并寫成歐拉形式

方程（3）中的ST就可以寫為

2 算法驗證

一端固壁一端自由的一維炸藥棒是一個模擬固體炸藥的經(jīng)典算例，這個問題的爆轟過程的數(shù)值模擬已經(jīng)十分成熟，本文中將它作為一個基準算例來驗證我們的程序。一根長0.1m的TNT炸藥棒從固壁端起爆，用燃燒函數(shù)［11］和JWL狀態(tài)方程［12］描述炸藥的爆轟過程和爆轟產(chǎn)物。TNT炸藥的材料參數(shù)如表1所示，其中ρ0為炸藥的密度，pCJ為炸藥的C－J爆壓，vd為炸藥的爆速，e為炸藥的內(nèi)能，A、B、R1、R2和ω為JWL方程中的系數(shù)。

表1 TNT炸藥和B炸藥的材料參數(shù)Table 1 Material parameters of TNT and composition B

圖1為起爆10μs后的壓力分布曲線，分別使用了5種不同的網(wǎng)格尺寸：在0.1m的長度上劃分200、500、1 000、2 000、4 000個網(wǎng)格。從圖中可以看到，隨著網(wǎng)格數(shù)n的增加，數(shù)值結(jié)果逐漸收斂。當網(wǎng)格數(shù)大于或等于1 000時，計算結(jié)果不受網(wǎng)格尺寸的影響。圖2為炸藥起爆后14μs內(nèi)不同時刻的壓力分布曲線，時間間隔為2μs。圖中的虛線是由實驗測得的TNT炸藥的C－J爆壓（21GPa），可以看到，隨著爆轟過程的推進，壓力峰值越來越接近C－J爆壓。

圖1 不同網(wǎng)格尺寸下起爆后10μs時的壓力分布Fig.1 Pressure distribution at 10μs after detonation with different resolutions

圖2 起爆后14μs內(nèi)不同時刻的壓力分布圖，時間間隔2μsFig.2 Pressure profiles along the TNT slab at an interval of 2μs until 14μs after detonation

3 數(shù)值結(jié)果與分析

對4種不同形狀的金屬襯墊進行模擬，考察金屬襯墊的形狀對射流的影響。圖3為初始時刻4種不同形狀金屬襯墊的構(gòu)型圖。圖3（a）～（c）是圓錐形金屬襯墊，它們的頂角θ依次為43.6°、77.3°、120°；圖3（d）的襯墊是球面的一部分，球面半徑為0.06m。4種構(gòu)型的金屬襯墊的裝藥半徑均為0.04m。裝藥采用B炸藥，其材料參數(shù)如表1所示。金屬襯墊全部為OFHC銅，使用Johnson－Cook本構(gòu)模型和 Mie－Gruneisen狀態(tài)方程模擬這種材料。炸藥的起爆方式可分為點起爆和面起爆，本文中全部采用面起爆。

圖3 初始時刻4種不同形狀襯墊的構(gòu)型圖Fig.3 Configurations of the shaped charges with different－shaped liners

圖4為不同形狀的金屬襯墊在裝藥起爆后50μs時的變形圖和密度分布等值云圖，圖4（a）～（c）是圓錐形襯墊的計算結(jié)果，圖4（d）是球面襯墊的計算結(jié)果?？紤]到球面襯墊的頂角為180°，可以看到，隨著金屬襯墊頂角的增大，射流的長度逐漸變短。這是因為頂角為43.6°的襯墊產(chǎn)生的射流速度最大，在同樣的時間里射流的變形最大，被拉伸的長度最長。此外，還可以看到，隨著襯墊頂角的增大，射流回流部分（即翅膀以上部分）所占的比重逐漸減小，也就是說，隨著襯墊頂角的增大，金屬襯墊的利用率有所增加。從密度云圖可以看到，射流后部的表面，也就是與炸藥的接觸面附近密度最大，是由于受到炸藥的爆炸載荷作用后，材料急劇壓縮；而前部的射流密度較低，這是由于射流被拉長導致局部的體積膨脹。

圖4 t＝50μs時不同形狀的金屬襯墊形成的射流的變形圖與密度分布Fig.4 Deformation and density contour plots of the jets at 50μs from liners with different shapes

在程序中監(jiān)測射流物質(zhì)的最大速度和最大溫度并每隔1μs記錄1個值。圖5（a）、（b）分別為4種不同形狀的金屬襯墊形成的射流的速度最大值和溫度最大值的時間歷程曲線?？梢钥吹?，隨著襯墊頂角的增大，射流速度逐漸減小。溫度的時間歷程曲線略有不同，在受到爆炸載荷作用的初期，襯墊急劇壓縮，塑性變形產(chǎn)生塑性溫升，密度增大產(chǎn)生沖擊溫升，材料內(nèi)部的溫度迅速達到最大值。隨著射流的形成，盡管塑性變形仍然使得溫度上升，但是由于射流的伸長，材料密度快速下降，沖擊溫升部分開始貢獻負溫升。因此，在爆炸載荷作用的初期，材料的溫度隨著襯墊頂角的減小而增大，而當射流形成并趨于穩(wěn)定之后，4種襯墊頂角的射流溫度十分接近。

圖5 不同形狀襯墊下射流內(nèi)部速度最大值和溫度最大值的歷史曲線Fig.5 Maximum－velocity histories and maximum－temperature histories of the jet from liners with different shapes

4 結(jié) 論

模擬了由4種不同形狀襯墊構(gòu)成的成形裝藥，通過分析發(fā)現(xiàn)，在爆炸載荷作用的初期，材料的溫度隨著襯墊頂角的減小而增大，而當射流形成并趨于穩(wěn)定之后，4種構(gòu)型的襯墊產(chǎn)生的射流溫度十分接近。隨著襯墊頂角的增大，射流的長度逐漸變短，速度逐漸減小，對金屬襯墊的利用率有所增加。

本文中目前只完成了射流的形成過程，射流的穿甲問題是接下去需要開展的工作，這涉及到另一個內(nèi)容十分豐富的領域——穿甲力學。此外，還可以考慮將現(xiàn)有的計算平臺向三維擴展。

［1］Liu Wei－dong，Liu Kai－xin，Chen Qian－yi，et al.Metallic glass coating on metals plate by adjusted explosive welding technique［J］.Applied Surface Science，2009，255：9343－9347.

［2］Chen Qian－yi，Liu Kai－xin.A high－resolution Eulerian method for numerical simulation of shaped charge jet including solid－fluid coexistence and interaction［J］.Computers ＆ Fluids，2012，56：92－101.

［3］Udaykumar H S，Tran L，Belk D M，et al.An Eulerian method for computation of multimaterial impact with ENO shock－capturing and sharp interfaces［J］.Journal of Computational Physics，2003，186（1）：136－177.

［4］Tran L B，Udaykumar H S.A particle－level set－based sharp interface Cartesian grid method for impact，penetration，and void collapse［J］.Journal of Computational Physics，2004，193（2）：469－510.

［5］Ma S，Zhang X，Lian Y，et al.Simulation of high explosive explosion using adaptive material point method［J］.Computer Modeling in Engineering ＆ Sciences，2009，39（2）：101－123.

［6］Teng X，Wierzbicki T.Evaluation of six fracture models in high velocity perforation［J］.Engineering Fracture Mechanics，2006，73（12）：1653－1678.

［7］Meyers M A.Dynamic behavior of materials［M］.New York：Wiley，1994.

［8］湯文輝，張若棋.物態(tài)方程理論及計算概論［M］.長沙：國防科技大學出版社，1999.

［9］Wang Jing－tao，Liu Kai－xin，Zhang De－liang.An improved CE／SE scheme for multi－material elastic－plastic flows and its applications［J］.Computers ＆ Fluids，2009，38（3）：544－551.

［10］Chen Qian－yi，Wang Jing－tao，Liu Kai－xin.Improved CE／SE scheme with particle level set method for numerical simulation of spall fracture due to high－velocity impact［J］.Journal of Computational Physics，2010，229：7503－7519.

［11］李德元，徐國榮，水鴻壽，等.二維非定常流體力學數(shù)值方法［M］.北京：科學出版社，1998.

［12］Lee W H.Computer simulation of shaped charge problems［M］.World Scientific：Singapore，2006.