炸藥裝藥損傷行為數值模擬研究進展

張萌昭，屈可朋，沈 飛，吳翰林，周 濤

(西安近代化學研究所， 西安 710065)

炸藥是武器戰斗部的核心能源，是為武器提供殺傷性、破壞性、動力性的關鍵材料，其安全性也一直備受關注[1]。然而，在生產、運輸、發射過程中，由于加載沖擊[2]、應力狀態、環境溫度[3]、生產工藝等各方面復雜因素的影響，炸藥會出現各種形式的損傷，其力學性能、感度、破壞性等特性也會因此產生變化，甚至導致炸藥提前引爆，產生不可預估的后果[4]。為了保障戰斗部能夠安全穩定地實現打擊效果，必須對炸藥的損傷行為進行系統而深入的研究。

目前國內外對炸藥損傷行為的研究主要依靠兩種方法，一種主要依賴于實驗獲得相關數據，如目前常用的Hopkinson壓桿實驗、輕氣炮實驗、巴西實驗[5]、落錘試驗[6]等方式。實驗方式能夠真實、準確地得到不同炸藥在不同條件下的應力應變關系及損傷情況，但是由于炸藥的不穩定性，這種研究方法具有一定的危險性。此外試驗成本、實驗條件也往往會限制試驗過程的開展，得到的實驗結果一般不具有廣泛性和預測性，對于炸藥細微觀損傷的實時數據采集也很難實現。另一種則是數值模擬方法，國內外許多學者對炸藥的損傷特性進行了模擬研究，基于LS-DYNA等軟件，可以同時得到空間和時間尺度上的仿真數據，進而能夠觀測到材料細微觀尺度的內部變化過程。但是這種方法的準確性和適用范圍常受到材料的本構模型和模擬計算方法的影響。為了對目前炸藥損傷問題中常見的本構模型及算法進行詳細分析，列舉了國內外炸藥損傷本構模型及炸藥裝藥力學行為的數值計算模型，對比分析其優缺點及適用范圍，以期為研究人員選擇本構模型和仿真算法提供參考。

1 炸藥含損傷本構模型研究現狀

材料的本構關系能夠通過數學表達式來描述材料在外界因素作用下的力學行為及變形規律。通過給定相關參數，通過本構模型即可得到細致、準確的力學狀態。下面介紹幾種損傷本構模型。

1.1 粘彈性統計裂紋模型

(1)

其中:K為含裂紋擴展的等效力強度因子;K0為材料的斷裂韌性；vmax為最大擴展速度；m為速度系數。

(2)

其中：c為裂紋體的平均半徑；G為彈性元的剪切模量；τ為松弛時間；a為初始缺陷尺寸；n和N分別代表第n個彈性元、體元數量。

Vioco-SCRAM模型與SCRAM模型相比結構簡單、參數少，因此應用更加廣泛。除了可以計算裂紋的生長過程外，還可用于模擬炸藥的非沖擊點火過程，如孫寶平就基于該模型計算了侵徹過程中裝藥的整體升溫和點火機制[11]。但是該模型考慮的因素較多，各參量的確定雖然與SCRAM模型相比較為簡單，但仍然是一項繁瑣的工作。此外，該模型得到的計算結果與實驗結果相比，在某些情況下仍然有較大差距，而且對于炸藥細微觀的損傷現象表述仍不夠成熟和完善，需要進一步發展和改進。

針對Visco-SCRAM模型的不足，一些研究人員改進了模型參數的選取方式及參數取值范圍，提高計算精度的同時簡化了參數設置。如張延耿等[12]，基于幣形微裂紋面的對稱特性，利用半球體而非球體來描述裂紋的取向，將傳統Visco-SCRAM中微裂紋方向矢量的取值范圍縮小為原來的1/2。為了適應本構參數的變化，同時對初始缺陷尺寸表達式進行修正，得到了新的偏量本構關系率的表達形式，以不同的初始缺陷尺寸參數對應不同的應力狀態。此外，考慮到拉伸條件下線性的體量本構關系不能夠描述體積應變增量，因此結合文獻[10]中給出的微裂紋擴展速度方程對體量本構關系進行改進，最終得到改進后的本構關系表達式。對該模型利用LS-DYNA程序進行仿真，通過對比PBX炸藥的平板撞擊試驗數值模擬的計算結果，發現改進后的模型能夠更加合理、細致地反應炸藥在低速沖擊作用下的細微觀損傷破壞機制，但是該模型對一些參數有較強依賴性，數值模擬結果仍存在一定差異。

1.2 含裂紋擴展損傷的粘塑性本構方程

現在國內外的微裂紋統計損傷模型，研究的重點較多基于炸藥的粘彈性行為，但一些力學性能試驗表明，炸藥在外力的作用下也會表現出明顯的塑性特性?；谶@種現象，成麗蓉等[13]建立了能夠反映PBX材料的粘塑性特性的含損傷本構方程。該模型同樣假設初始狀態的微裂紋各向同性分布，并將一個阻尼器和塑性單元并聯以后與彈性體串聯以構建彈粘塑性模型，如圖1。

圖1 彈粘塑性模型示意圖

得到的偏應力張量表達式如下：

(3)

為了驗證該模型的準確性，成麗蓉等人對PBX炸藥試件進行了靜動態單軸壓縮力學實驗及斷裂性能試驗，得到的試驗數據與仿真結果進行比較，發現雖然該模型能夠描述PBX炸藥的粘塑性及脆性損傷力學的響應過程，但是模擬曲線在應力強度超過最大值以后裂紋擴展速率迅速增大，模擬曲線與實驗數據并未完全吻合，尤其在應變率較高的情況下擬合效果有較大偏差。

1.3 含各向異性損傷的粘彈性統計微裂紋本構模型

除了上文提出的基于初始材料各向同性假設的本構模型外，也有一些模型考慮了微裂紋法向造成的影響，提出了含各向異性損傷的粘彈性統計微裂紋本構模型(Aniso-Visco-SCRAM)。如圖2。

圖2 Aniso-Visco-SCRAM模型示意圖

張延耿等[14]改變了原本Visco-SCRAM模型中利用積分形式計算微裂紋擴展產生的應變的方法，利用疊加法將不同法向上的微裂紋應變分量求和，得到宏觀上的微裂紋擴展產生的附加應變。與傳統Visco-SCRAM模型得到的裂紋在各個角度擴展情況相同相比，Aniso-Visco-SCRAM更好地描述了PBX炸藥在低速沖擊載荷下的損傷機理和變形特征。但是該模型所需的部分參數無法通過試驗確定，而且這些參數的取值對模型的準確性有較大的影響，所以需要更進一步的定量研究以更好地描述損傷破壞過程。Marwen Chatt等[15]也對含各向異性損傷的本構模型進行了研究。他所提出的模型是一種新的粘彈性微平面模型，使用Maxwell鏈模型和微平面方法并與含各向異性損傷結合來模擬準脆性PBX材料的行為。該模型基于Carol等人提出的V-D-T分解方法，在微平面上將應變張量沿微平面方向投影，之后采用內部微平面法，在每個方向都定義一個Maxwell鏈模型，該模型由10個粘彈性體和1個彈性體并聯而成。V-D-T分解方法如圖3所示。

圖3 V-D-T分解示意圖

定義微平面內部的行為需要確定每個微平面方向和每個元素的彈性應變和彈性應力。并且該模型認為損傷的演化僅與彈性元件的相互作用力有關而與粘彈性元件無關。與實驗數據相比，當考慮包絡曲線、橫向應變和松弛步階時，數值預測與實驗數據擬合較好。但在描述壓縮的模擬卸載過程時，所提出的模型仍然缺乏準確性，沒有考慮材料的塑性和不可逆性的損傷。同時該模型所需的參數較多，需要高性能的計算機提供支持。

1.4 內聚裂紋模型

由于炸藥在某些情況下與混凝土等有相似的力學特性，也有一些由基于混凝土材料的本構模型發展而來的用于炸藥的本構模型，例如內聚裂紋模型。內聚裂紋模型可以計算材料任意位置的裂紋的產生或擴展，該模型所需的材料參數少，物理意義明確。該模型的特點在于能夠避免處理裂紋尖端的奇異的特點，通過單元函數的擴展而直接作用于網格內部。內聚裂紋模型的基本方程如下：

(4)

崔云霄[16]等人基于內聚裂紋模型對PBX模擬材料進行了動態巴西試驗的數值模擬，利用LS-DYNA軟件得到變形場及應變測試結果，得到拉伸應變的峰值誤差小于5%，對于裂紋張開角度的模擬誤差也能控制在1°之內，因此內聚裂紋模型能夠較好地模擬PBX-M材料在動態巴西實驗過程中的裂紋擴展過程。吳艷青[17]等人也采用了內聚力模型來研究PBX的損傷，該模型的特點在于分別考慮了粘彈性和脆性損傷的來源，認為炸藥的粘彈性主要來自粘結劑，因此基于HMX的彈塑性模型將粘結劑的功能通過內聚力模型表現，進而直接置于HMX晶體的外界面。

1.5 Z-W-T模型

我國國內使用較多的本構模型由朱兆祥、王禮立和唐志平等人，基于Green-Rivlin本構理論提出的Z-W-T模型。該模型將一個非線性彈簧和兩個Maxwell體并聯進而描述高聚物的非線性粘彈性力學行為。其方程形式為：

(5)

其中:E0、α、β、E1、E2、θ1、θ2均為材料常數。但是這種模型難以用于描述PBX炸藥的損傷情況。因此孫文旭等[18]將損傷行為引入Z-W-T模型，建立了含損傷的非線性粘彈性本構模型。他們對一種新型的抗過載澆鑄PBX炸藥進行準靜態力學實驗和分離式霍普金森壓桿實驗，通過遺傳算法擬合，對高應變率下PBX炸藥在到達破壞前的動態力學行為進行了較好的描述。屈可朋等[19]采用含顆粒相結構應變率效應的Z-W-T模型對RDX基含鋁炸藥的沖擊特性進行研究，與通過分離式霍普金森壓桿實驗得到的實驗數據擬合較好，表明該模型可以較好地描述RDX基含鋁炸藥在高應變率下的損傷行為。

2 數值模擬方法研究現狀

隨著計算機技術的發展，炸藥損傷的數值模擬方法得到了極大的發展，一定程度上彌補了傳統試驗方法的限制和不足，得到了廣泛應用。數值模擬方法能夠全方位且直觀精確地反應炸藥在不同尺度上的力學響應特性。下面將對一些數值模擬方法進行介紹。

2.1 有限元方法

有限元方法是一種常用的研究PBX力學行為的方法，這種方法基于連續介質力學，能夠在宏觀尺度上研究材料的力學性能，其主要思想在于通過變分方法求解出偏微分方程的近似解，求解時將物體離散化，并對其中每個微小單元進行求解，并通過單元與單元之間的關系進行約束，具有計算速度快，而且相關的計算軟件較為成熟的優點。石嘯海等[20]就利用有限元方法，針對戰斗部在侵徹過程中產生的裝藥損傷問題進行研究，分別對兩組模型進行了仿真模擬，計算了PBX的裝藥損傷演化過程，并且分析了裂紋寬度、過載、軸向應力和損傷程度等因素，比較了有機玻璃、聚四氟乙烯、酚醛樹脂、尼龍四種材料的緩沖性能，驗證了內聚力本構模型的準確性。魏強等[21]則基于 Karagozian & Case模型，利用有限元方法對PBX炸藥在撞擊作用下的損傷情況，引入了特征長度以保證計算精度。但是有限元方法也存在不足，在有限元方法中，網格的邊界與模型的幾何邊界一致，因而在變形過程中會出現網格的畸變、纏繞等影響計算精度的問題，同時對于裂紋的生長過程的研究也有一定限制，很難對結構的細微觀變化進行準確描述。

2.2 擴展有限元方法

為了彌補有限元方法的不足，發展出了擴展有限元方法(XFEM)。擴展有限元方法由Belytschko和Black[22]首先提出，該方法改變了傳統有限元方法的劃分方法，并不對結構內部的裂紋和不連續面等問題進行網格劃分，而是將增強形函數加至傳統有限元的位移模式中來代表計算區域內的間斷，圖4表示裂紋平板有限元模型，其中淺黑色部分為裂紋，圖中方框節點表示被完全切斷的集合，而圓圈節點則表示圍繞裂尖的節點。當網格密度足夠時，可以更加方便地模擬材料中的裂紋、空洞等情況[23]，具有計算精度高、無需重新劃分網格等優點。

圖4 裂紋平板有限元模型示意圖

戴開達等[24]用擴展有限元方法對PBX炸藥在多種巴西實驗中的變形破壞過程進行了研究，當單元被裂紋完全穿過時使用的增強型函數為

ψJ(x)=NJ(x)H(f(x))

(6)

其中NJ(x)為通用增強型函數，H(x)在(x-x*)·n≥0取值為1，否則為-1。f(x)的計算式如下：

(7)

當節點處于裂尖周圍時，則增強形函數為

(8)

基于上述2個增強形函數，可用下式來表示含裂紋二維板的位移場：

(9)

其中:S為所有節點的集合；Sh為被裂紋完全切斷的單元的節點集合；Sc為圍繞裂尖的單元的節點；aJ、bK為節點的附加自由度。

利用上述公式，基于擴展有限元法模擬了巴西實驗中裂紋起裂、擴展和斷裂等行為，并且得到了在不同加載形式下圓盤式樣的不同應力應變分布。

Huang X C等[25]則基于擴展有限元模型和內聚模型，對含孔洞平板的力學性能進行研究，發現對該平板施加整體壓力時，孔洞周圍的局部拉伸力會導致局部裂紋產生。并且通過數值模擬得到了裂紋在整個試驗過程中的整體走勢、拐點，包括起裂時刻、裂紋初期展速等數值，證明了有限元方法和內聚模型的準確性。但是擴展有限元方法在解決三微裂紋擴展和群裂紋等方面仍然存在缺陷。

2.3 物質點法

上面介紹的方法計算過程依賴于網格的劃分，存在一定的局限性，因此近些年來無網格法被提出并且逐漸發展。物質點法(MPM)就是常用的無網格法之一。MPM與網格法本質上的區別是物質點法是一種基于粒子的計算方法[26]，MPM首先將連續體通過某種方法離散成一個個具有集中質量的物質點，在計算過程中物質點的質量保持不變，之后通過建立背景網格節點上的有限元函數實現質點和背景的固連關系。他們之間的映射關系表示為：

(10)

其中：ψp為物質點的變量；ψi為網格節點變量；nu為單元節點數；xp為物質點坐標。在此映射關系下，動量方程的虛功形式可以表示為離散形式：

(11)

(12)

通過上式和材料本身的本構關系方程即可計算得到物質點的應力，并且物質點法的背景網格隨時間不斷更新，避免了有限元法的缺陷，同時在計算過程中MPM也滿足了質量守恒定律和動量守恒定律[27]。唐紅[28]等人基于物質點法，對奧克托今基高聚物粘結炸藥的力學性質進行模擬，并且建立了炸藥顆粒隨機分布的力學模型，對PBX炸藥壓制過程中的二維數值模擬進行模擬，得到了炸藥在炸藥過程中顆粒和縫隙的變化過程，同時還檢測了炸藥體系溫度的變化情況。MPM與其他方法相比，有著較高的效率、精度和穩定性。但是由于質點積分存在誤差，對于小變形問題，MPM的精度要低于拉格朗日有限元方法，因此該方法更加適用于大變形，

2.4 無網格伽遼金方法

無網格伽遼金方法(EFG)是另外一種使用較為廣泛的無網格方法。這種方法是由Belytschko T[29]提出，對形函數的導數和邊界條件的施加方法進行改進，提高了計算的收斂速度和精度，彌補了傳統有限元方法中容易出現的體積鎖死現象的問題。其具體方法是通過最小二乘法構建基于一定已知點構造出區域內的近似函數，即：

(13)

其中：u(X)為計算域內函數；p(X)為基函數；a(X)為J(X)的極小值，與X相關，其中J(X)為：

(14)

其中權函數wi(X)有一個影響范圍，其值非負，當超出該范圍時其值取0。

崔云霄[30]等人基于EFG對PBX炸藥的圓弧巴西實驗進行了數值模擬，得到了x方向的應變峰值約為6%，與實驗結果較一致，較為準確地得到了準靜態壓縮下裂紋的產生即擴展過程。

2.5 離散元方法

離散元方法(DEM)是一種較為成熟的數值模擬方法。這種方法是由Cundall[31]最早提出，基于牛頓第二定律，將塊體單元設定為理想剛體，并且基于塊體單元的幾何形狀和變形特征得到不同的形狀、性質離散元，用顯式差分的計算方法，先得到各個元素的運動方程，之后通過進一步計算得到物質的整體運動形態，適用于非連續截止問題即非連續體破壞問題。在炸藥的沖擊過程中，對細觀模擬能夠取得不錯的計算結果和精度。金博[32]在研究炸藥的反應速率與裂紋的擴展和生長之間關系時，采用DEM方法進行數據模擬，得到了不同發射速度下的裂紋分布、成長狀況，進而利用裂紋成核理論描述反應過程，進一步得到了較為完整的反應速率模型。但是這種方法的精確度對主要參數的依賴性較高，計算效率較低。

2.6 近場動力學方法

此外近些年也有學者基于近場動力學來研究炸藥損傷問題。近場動力學方法(PD)最早由美國的Silling[33]提出，這種方法和MPM相似，首先將固體離散為一系列帶質量的物質點，但是PD考慮了物質點和其相鄰的一定范圍內物質點之間的相互作用力，如圖5，這種相互作用力包含材料的物性信息，從而避免了傳統的應力-應變形式，進而不需要連續性假設。

圖5 物質點的運動及相互之間作用力示意圖

但是以鍵為基礎的動力學過于簡化，在使用過程中受到限制，故PD理論又進一步發展成為以狀態為基礎的近場動力學理論Silling[34]。李潘[35]基于狀態PD理論對裂紋在PBX炸藥裂紋擴展中的應用進行了研究，使用了顯式的計算方法，通過巴西圓盤實驗、帶長方孔平板拉伸實驗、含裂紋板動態斷裂試驗等多種試驗研究PBX炸藥的力學行為，建立了近場動力學蠕變模擬方法。秦洪遠等[36]則通過非局部鍵形近場動力學建模，對含初始裂紋的巴西圓盤劈裂過程中裂紋的擴展及破壞形式進行研究，形成了較為完整的近場動力學數值體系，得到了不同初始裂紋長度、不同初始中心裂紋傾角下的圓盤裂紋生長情況。PD方法不需要借助外部準則，避免了裂紋尖端奇異性問題，同時提高了計算效率。

對以上幾種數值模擬方法進行總結，具體情況如表1所示。

表1 幾種數值模擬方法

3 結論

針對炸藥的損傷，列舉了一些較為常用的含損傷本構模型，介紹了一些數值模擬方法。由于炸藥較為復雜的力學特征，很難得到一個廣泛適用的本構模型，雖然所介紹的本構模型都可以用來描述炸藥的損傷行為，但是其側重點和適用范圍及驗證方法都不盡相同。含損傷的本構模型在細微觀層面的發展還有欠缺。

一般受限于實驗條件，對炸藥的損傷進行全面系統的受力情況、力學性能分析是比較困難的，因此數值模擬方法在實際工程中相當重要。但是現有文獻介紹的數值模擬方法在多尺度分析方面有很大不足，而且對于數值模擬的理論分析較為欠缺。

為了更好地研究炸藥裝藥損傷，需要在細微觀層面進行更加深入的探討和研究，構建更加合理、準確的含損傷細微觀本構模型，提升數值模擬精度和計算效率。