高聚物粘結炸藥動態損傷破壞的數值刻畫

魏強， 黃西成， 陳剛， 陳鵬萬

(1.北京理工大學 爆炸科學與技術國家重點實驗室, 北京 100081; 2.中國工程物理研究院 總體工程研究所, 四川 綿陽 621999)

0 引言

顆粒粘結炸藥在低速撞擊等環境下可能發生點火，一般情況下，不會發展成為爆轟，但仍會對人員和器材設備形成威脅[1]。高聚物粘結炸藥(PBX)的這種機理尚不明確的點火行為與其內部損傷的形成息息相關[2-4]，刻畫PBX炸藥的損傷及其演化過程存在著很多困難。從試驗的角度來講，炸藥內部損傷過程難于捕捉，且存在隨機性(由于初始缺陷的隨機性分布)；此外，在描述炸藥的破壞時也缺乏合適的表征參量，來建立材料細觀破壞形式與材料宏觀力學響應之間的聯系。從數值仿真的角度來講，首先需要解決的就是所用數值方法本身在描述材料宏、細觀破壞時存在的局限性以及伴隨引入的虛假數值結果圖像(如應變局部化問題)；其次需要解決的就是建立能正確刻畫材料力學行為的本構模型。由于存在上述困難，且不易解決，因而現今對PBX炸藥損傷破壞過程的研究主要是捕捉一些破壞圖像，增進對相關破壞過程的認識[5-7]。

相比于試驗手段，數值仿真能提供更多的過程信息，這對炸藥安全性的評估是有利的；商業有限元軟件的廣泛應用也為數值仿真研究提供了很大的便利。但PBX炸藥具有明顯的黏彈性，拉伸、壓縮不對稱性，壓力相關性，應變率相關性，溫度相關性等[8]，此外，其變形破壞過程中可能還會伴有明顯的化學反應過程[9]，這些復雜的特性使得PBX炸藥的力學行為難于刻畫，更沒有現成的本構模型可用。由于PBX炸藥的宏觀力學響應形態與混凝土材料類似，而混凝土材料模型經過多年的發展已經形成了一套較為完備的描述體系，因此借鑒混凝土材料模型來描述PBX炸藥的力學行為成為最常用的方法之一[10-12]。相較而言，PBX炸藥的力學行為比混凝土材料更為復雜，借鑒混凝土材料的本構模型主要是用來刻畫PBX炸藥在復雜應力及應變率影響下的力學行為；為了更真實地描述PBX炸藥的力學行為，仍需要對現有的混凝土本構模型進行進一步的發展。

本文在有限元分析軟件AUTODYN中，以Karagozian & Case (K&C)模型[13]框架為基礎，依據PBX炸藥的力學特性對模型進行了修改和擴展。在修改的模型基礎上，分析了材料的應變率效應及損傷破壞的描述方法等，以期較好地描述低速撞擊等環境下PBX炸藥的損傷破壞，為后續炸藥安全性評估研究提供基礎。本文的研究對象為法國原子能委員會研究人員研究的某奧克托今(HMX)基PBX炸藥，其力學行為類似于PBX9501[11].

1 本構模型

1.1 壓力依賴性

本文采用的模型考慮了PBX炸藥力學行為的壓力依賴性。模型以K&C模型的框架為基礎，修改了原K&C模型中的損傷演化模式，計及了材料強度、模量、損傷演化模式等的壓力依賴性。

與很多地質類材料本構模型的構造方法相同，模型采用了不同強度面插值的方法來描述材料當前的屈服強度，強度面間的插值函數是用戶輸入的損傷演化曲線。此處采用一條應力σ- 應變ε曲線來說明模型的描述方法。用于插值的不同強度面在應力- 應變曲線中體現為不同的特征點，如圖1所示。模型中采用3個特征點來刻畫材料的應力- 應變曲線，分別是初始屈服點(材料線性彈性段和非線性硬化段的臨界點)、最大應力點(材料硬化過程和軟化過程的臨界點)、殘余應力點(表征材料破壞后的殘余承載能力)。

圖1 材料應力- 應變曲線示意圖Fig.1 Stress-strain curve

在初始屈服點之前材料處于彈性變形階段，當應力超過初始屈服點后，模型即認為材料開始了損傷演化；當前屈服點由最大應力點與初始屈服點或殘余應力點通過非線性插值來確定。損傷演化曲線作為非線性插值函數參與了當前屈服點的確定，其主要作用是建立應力和應變之間的聯系，它是由材料單軸壓縮應力- 塑性應變曲線(初始屈服點至殘余應力點段)做無量綱處理后獲得的。為了考慮復雜的應力狀態，3個特征點的應力信息在應力空間內形成了3個強度面，而應變信息則是通過損傷演化參量來考慮應力狀態的影響的。通過損傷演化曲線應力和應變信息之間建立起了聯系，這樣，模型就能完整地描述材料在不同應力狀態下的應力- 應變響應。

3個特征點在空間形成的最大強度、殘余強度、初始強度3個強度面可以分別表示為

(1)

(2)

(3)

式中：a0、a1、a2、a1f、a2f、a0y、a1y、a2y均為待定參數；p為靜水壓力(取壓縮為正，拉伸為負)。

材料當前的屈服強度可通過上述3個強度面插值得到，用方程表述為

(4)

式中：λ為損傷演化參量；η(λ)為損傷演化曲線(用于非線性插值)；λ1對應于η(λ)剛達到1時的λ值。損傷演化參量λ的表達式為

(5)

式中：εp為塑性應變；r為應變率增強因子；εm為不同約束壓力pa下材料的最大失效應變(對應于不同約束壓力下材料的最大強度)，與壓力相關。

在約束壓力較低時，上述描述方法已經可以較好地描述材料的應力- 應變響應。但當約束壓力大于脆延轉化壓力后，試驗結果顯示，材料的應力- 塑性應變曲線不再出現軟化現象(見圖2)。這時，當λ>λ1后，材料的損傷演化過程就無法再采用(4)式來描述，需要重新給出描述方法。為了解決這個問題，需要對K&C模型進行修改，具體為在高約束壓力下，取消原模型中的軟化插值過程，引入一個簡單的線性硬化過程，即認為在較高約束壓力下，材料后期會進入一個線性硬化階段，模型也直接變為塑性線性硬化模型(硬化模量依賴于約束壓力)。這樣，以材料的脆延轉化壓力為界，模型中相當于構造了兩種類型的損傷演化曲線，用來描述材料在不同約束壓力下的應力- 應變響應(存在軟化現象和不存在軟化現象)。

圖2 不同約束壓力下PBX炸藥的應力- 塑性應變曲線[9]Fig.2 Stress-plastic strain curves of PBX explosives under different confined pressures[9]

在選擇采用的損傷演化曲線時，模型中采用了約束壓力的概念，與靜水壓力不同，約束壓力不會隨著單向壓縮加載而發生變化，如單軸壓縮過程中，約束壓力始終為0 MPa. 約束壓力概念的引入對描述材料模量的壓力相關性至關重要，顯而易見，材料模量的變化不可能與靜水壓力形成一一對應的關系。模型中，約束壓力定義為

pa=p-Δσa/3，

(6)

式中：Δσa為單向壓縮加載路徑上(單軸壓縮與靜水壓力組合形成的加載路徑)的von Mises等效應力值。在單向壓縮情況下，這種定義無疑是正確的；對于其他復雜應力狀態，約束壓力的定義是一種等效。上述模型成功地捕捉到了PBX炸藥在不同約束壓力下的應力- 應變響應。

1.2 應變率效應

圖3 不同應變率下材料強度隨壓力的變化情況[14]Fig.3 Variation of material strength with pressure under different strain rates[14]

本文中采用徑向應變率增強的方式來描述材料的應變率效應，通過圖3中的字母標識，對這種應變率增強方式進行說明：假設點C當前的壓力為G，模型中認為這是應變率增強后的壓力，首先通過當前的應變率增強因子(通過當前應變率得到)把壓力縮小，得到對應準靜態時的壓力D，進而得到準靜態時的強度E；再通過應變率增強因子放大，得到應變率增強后的強度F，作為點C的強度值。由幾何相似性不難證明，當材料強度隨壓力的變化為直線時(如圖3中所示的情況)，徑向應變率增強方式得到的結果與圖3中的試驗結果相一致；但是從本文采用的試驗數據[9,11]來看，在更寬的壓力范圍內，PBX炸藥的強度隨壓力近似呈雙線性變化，Wiegand等[15]就曾建議直接采用雙線性方程來描述PBX炸藥強度隨壓力的變化情況。為了方便介紹，圖4中比照圖3給出了較大壓力范圍內PBX炸藥強度隨壓力變化的示意圖，線段na和ab對應于準靜態時的情形，點n對應的強度為材料的單軸壓縮強度。

圖4 應變率效應修正說明示意圖Fig.4 Schematic diagram of correction of strain rate effect

如圖4所示，在較高約束壓力下(壓力達到ab段時)，PBX炸藥材料強度隨壓力變化的斜率會減小，這時，徑向應變率增強的方式會高估材料的應變率增強效應(徑向應變率增強的方式等價于斜率不變的情形)。本文對較高壓力下材料的應變率增強效應進行了修正，徑向應變率增強的方式形式簡潔，實現方便，本文不作改動，修正的主要是應變率增強因子。在圖4中，延長ab段與cd段(單軸壓縮路徑)相交于點m，對mb段而言，這就構造出了與圖3類似的情形，不過此時交點m對應的強度不再是材料的單軸壓縮強度。對應于ab段新的應變率增強因子可表示為

rn=1+(ro-1)fc/σm，

(7)

式中：fc為材料的單軸壓縮強度；ro為未修正前的應變率增強因子；σm為點m對應的強度值。本文中，壓縮時的動態增強因子rc采用了文獻[14]中的試驗數據：

(8)

由于沒有在公開文獻中查到法國原子能委員會研究人員在動態拉伸方面的研究工作，因而拉伸時的動態增強因子rt采用了文獻[16]中與之力學行為類似的PBX9501試驗數據：

(9)

除了材料強度與應變率相關外，本文中還假設材料的剪切模量與應變率增強因子成正比；而對于材料的最大失效應變，本文中認為其與應變率無關[17]。圖5給出了不同應變率增強因子下，材料的單軸壓縮應力- 應變曲線(由單個六面體單元計算輸出)，計算結果可以直觀地展示出上述假定。

圖5 單軸壓縮應力- 應變曲線隨動態增強因子的變化Fig.5 Variation of stress-strain curves in uniaxial compression with dynamic enhancement factor

1.3 偏平面形狀的描述

材料偏平面形狀由類三角形向圓形轉變的過程，也是材料由低約束壓力下的脆性行為向高約束壓力下的塑性行為轉變的一種體現。K&C模型中經驗性地給出了混凝土材料拉伸、壓縮子午線之比(偏平面形狀變化的依賴量)隨壓力變化的公式；把PBX炸藥的材料參數帶入到公式中發現，在單軸壓縮時，材料的拉伸、壓縮子午線之比已經達到0.90，即偏平面形狀已接近圓形，顯示出了明顯的塑性材料的特性，這與PBX炸藥的力學特性不相符(在單軸壓縮情況下，PBX炸藥材料體現出明顯的脆性特性)。

對于PBX炸藥在復雜應力下力學特性的試驗研究還很缺乏。唐維[18]對某PBX代用材料開展了相關試驗研究，根據試驗結果，在低約束壓力下，材料拉伸、壓縮子午線之比近似為2/3[19]，這與Riedel-Hiermaier-Thoma模型中建議的0.680 5相近[20]。本文中取0.68作為約束壓力較低時，材料拉伸、壓縮子午線的比值。隨著約束壓力變大，材料偏平面的形狀逐漸變為圓形，即拉伸、壓縮子午線之比逐漸變為1.00，本文中把炸藥晶粒屈服強度的值(HMX晶粒的屈服強度為130 MPa[21])作為此轉變壓力。這里隱含一假設，即炸藥晶粒本身的屈服強度沒有明顯的壓力依賴性。PBX炸藥在受拉時體現出明顯的脆性特性，因而當壓力為負時，本文沿用K&C模型中的建議，把拉伸、壓縮子午線之比取為0.50. 本文中采用的拉伸、壓縮子午線之比可歸納為

(10)

式中：σY為炸藥晶粒的屈服強度。各點間采用線性插值過渡。

2 結構破壞中材料損傷的描述

2.1 變形局部化

采用有限元方法描述材料損傷和破壞過程時存在困難，面臨的困難之一是對材料軟化過程地描述。材料在應變軟化過程中會出現變形局部化的問題，這在巖土材料本構模型的研究中有過很多的討論[22]，其本質是軟化過程中材料變形地高度不均勻化。試驗中得到的材料應力- 應變響應的軟化段既包含了材料本身的力學特性，同時也包含了變形局部化后引入的結構效應。Wang等[23]就曾采用局部化結構特征與材料力學特性相結合的方法對某PBX代用材料的單軸壓縮應力- 應變曲線進行了描述，與試驗結果符合得很好。遺憾的是，大多數試驗研究中，一般只給出了材料的應力- 應變曲線，很少對試件破壞后的形貌進行描述，這對試驗數據的理解造成了一定的困難。

網格劃分的尺寸、計算中的數值誤差等同樣會導致數值計算中出現結構變形的局部化，這類似于試驗中出現的變形局部化，但產生的原因不同。為了減小這種不真實的變形局部化，除了要保證合適的計算網格尺寸外，還需要在單元中引入特征長度來解決。常見的引入方法有兩種[24]：方法1是采用非局部化模型的形式，其基本思想是，在給定的特征長度內(不受單元尺寸的約束)基于單元的積分點對塑性應變等參量進行平均化處理，從而引入特征長度的影響，避免數值上導致不真實的變形局部化；方法2是采用斷裂能加特征長度的形式，就其實際效果來看，相當于調整應力- 應變曲線的軟化段(下降段)，人為減緩單元的軟化速度，以此來避免有限元計算中出現的變形局部化問題。方法1需要計算程序中相關算法的支持，而方法2僅需調節材料的參數就可實現。

特征長度本意是用來定義材料變形局部化的尺度，如材料在破壞過程中形成的裂縫尺寸等。但在實際破壞過程中，裂縫尺寸受到應變率、應力狀態等很多因素的影響，導致很難給出合適的特征長度。本文中沒有直接引入特征長度的概念，而是通過調整應力- 塑性應變曲線的軟化段來近似考慮。結合模型的特點，調整時保持了應力不變，只對軟化段對應的塑性應變部分進行調節。采用單軸壓縮試驗進行校驗，壓縮試件分為6 mm×6 mm×12 mm的長方體和φ6 mm×12 mm的圓柱體兩種，壓頭和試件的接觸面不考慮摩擦。計算中，試件網格尺寸考慮了1.5 mm、1.0 mm和0.5 mm 3種。計算結果顯示，當下降段塑性應變放大為原來的4倍時，輸出結果與試驗結果取得了較好的一致性，如圖6所示。

圖6 單軸壓縮應力- 應變曲線計算與試驗結果的對比Fig.6 Comparison of calculated and experimental stress-strain curves of uniaxial compression

計算結果同時顯示出，即使網格的大小和壓縮試件的形狀有所不同，材料應力- 應變曲線的硬化段依然保持了很好的一致性，而軟化段就出現了明顯的分散性。相比而言，采用圓柱型試件得到的軟化段一致性要更好一些。

2.2 拉伸失效后的處理

對失效單元的處理是采用有限元方法描述材料損傷和破壞過程時面臨的另一個困難。單元的失效在程序實現時等同于單元的刪除，這會在結構中形成一個尺度虛假的裂紋，從而影響結構的后續響應。PBX炸藥易于產生拉伸斷裂，如果程序中不考慮拉伸破壞的話，就無法表征材料的破壞形貌；但如果僅把拉伸破壞的單元作常規失效處理的話，又會對結構后續的壓縮響應造成明顯的影響。為了更好地表征PBX炸藥材料，需要對拉伸失效后單元的處理方式進行修正。

AUTODYN軟件中有一種失效后可再次恢復承載能力的單元失效處理方式，遺憾的是在進行程序二次開發時這種失效方式無法直接調用。借鑒這個思路，本文對拉伸失效作如下處理：單元中的應力達到拉伸強度時，把單元的剪切模量和屈服強度置為小量(原數值的十萬分之一)，并把此單元設置為拉伸失效狀態進行輸出顯示；當拉伸失效單元重新承受壓縮載荷時，單元恢復壓縮承載能力，即認為拉伸失效不會對單元的壓縮過程產生影響。本文在描述材料的拉伸破壞時，沒有考慮拉伸過程中的損傷演化，認為當單元中的應力達到拉伸強度后，材料直接發生脆性斷裂。

采用單個六面體單元對程序的執行情況進行檢查。單元一端加載，另一端限制加載方向的自由度；加載端進行拉伸- 壓縮- 拉伸- 壓縮的反復加載過程，輸出的應力- 時程曲線如圖7所示。第1次拉伸過程，讓單元進入了拉伸失效，隨后的壓縮加載過程，單元完全恢復了承載能力，響應與材料的單軸壓縮過程完全一致；第2次進入拉伸過程后，單元中的應力首先進行了卸載，當應力卸載為0 MPa時，單元中的應力不會隨著拉伸載荷的增大而增長，即單元已經失去了拉伸承載能力；當單元再次進入壓縮過程時，單元先是經歷彈性加載，當應力到達當前的屈服強度后，應力響應沿著上次的壓縮損傷路徑繼續演化。單元是通過速度邊界條件進行加載的，加載過程沒有考慮應變率效應。圖7中標示的虛線位置為拉伸、壓縮加載發生轉化的時刻。

圖7 拉伸- 壓縮- 拉伸- 壓縮加載過程的應力- 時程曲線Fig.7 Stress-time history curve of tension-compression-tension-compression loading process

在實際的應用中發現，這種拉伸失效單元的處理方法易受單元變形情況及應力狀態的影響，在考慮拉伸裂紋擴展的算例中，該方法并不適用。

3 Steven試驗仿真計算

3.1 計算模型介紹

法國原子能委員會研究人員通過試驗手段研究了某PBX炸藥(力學行為類似于PBX9501炸藥)的力學特性及其安全性[9，11，14,25-26]，研究給出了不同約束壓力下PBX炸藥的應力- 應變響應，模量的壓力相關性，材料的應變率效應等一系列的試驗結果。在采用Steven試驗研究炸藥的點火問題時，還對未發生點火的炸藥試樣進行了解剖觀察，給出了試樣在撞擊后的破壞變形形貌。本文就采用這一系列的試驗結果，來考察模型對PBX炸藥在撞擊作用下損傷破壞情況的刻畫能力。

下面對Steven試驗的計算模型進行簡要的介紹，模型中：撞擊彈頭質量1.2 kg，撞擊速度為58 m/s，材料為鋼；靶體為φ200 mm×35 mm的圓柱型鋼盒，撞擊端面鋼層的厚度為3 mm，鋼盒內含有φ100 mm×13 mm的圓柱型PBX炸藥試件及內徑為100 mm、外徑為120 mm、厚度同樣為13 mm的有機玻璃環。計算模型如圖8所示，炸藥試件的網格尺寸控制在1 mm左右。計算中，鋼采用了理想彈性塑性模型來描述，有機玻璃采用了彈性模型來描述，具體參數參見文獻[11]。

圖8 Steven試驗計算模型Fig.8 Numerical simulation model of Steven test

為了更好地描述試驗數據，本文對PBX炸藥的強度面進行了分段擬合，炸藥本構模型中強度面的參數如表1、表2所示(在高約束壓力下材料不會發生軟化，因而表2中沒有殘余強度面的參數)，試驗數據參考文獻[9，11]。本文中采用的損傷演化曲線根據文獻[9]中的試驗數據獲得，得到的曲線如圖9所示，圖中η為采用材料單軸壓縮強度歸一化后的應力(硬化段部分取初始屈服強度時η=0，取最大強度時η=1)。計算中，炸藥的剪切模量和體積模量考慮了材料的壓力相關性，試驗數據取自于文獻[25]，如圖10所示。

表1 約束壓力小于100 MPa時強度面的參數

表2 約束壓力大于100 MPa時強度面的參數

圖9 損傷演化曲線Fig.9 Damage evolution curve

圖10 模量隨約束壓力的變化Fig.10 Variation of modulus with confined pressure

3.2 計算與試驗結果對比

撞擊后炸藥的變形形貌如圖11所示，在撞擊成坑部位，模型很好地捕捉到了坑體的形貌。同樣，模型也成功捕捉到了炸藥試件的隆起變形，但計算結果在數值上與試驗結果存在差距。隆起過程涉及到材料的破壞、滑移及由此引起的體積膨脹等問題，主要是材料的局部化變形過程，采用有限元方法進行描述時存在困難。雖然沒有很好地描述到材料的隆起變形，但模型仍捕捉到了試件中的一些破壞形態，如隆起過程中形成的剪切破壞、炸藥試件靠近有機玻璃部分底部的翹起現象等。

圖11 撞擊后炸藥試樣變形情況Fig.11 Deformation of explosive specimen after impact

文獻[15]給出了撞擊后試樣(同類PBX炸藥的代用材料)背面典型的破壞形態，數值計算結果與其對比如圖12所示。計算得到的破壞形態與試驗結果整體類似，顯示出輻射狀與環狀相疊加的破壞形態。從整體效果來看，本文采用的模型已較好地刻畫了炸藥試樣被撞擊后的破壞形貌。當然，炸藥試樣變形破壞成因復雜，尤其是撞擊中心部位(經歷過高約束壓力等復雜應力狀態)，本文提出的模型還不能正確描述，其中的變形破壞機理還需要進一步地認識和研究。

圖12 撞擊后炸藥試件背面的破壞情況Fig.12 Damage in the rear of the explosive specimen after impact

除上述試驗結果外，文獻[11]中還給出了一組PBX炸藥與鋼盒接觸界面處應力隨時間的變化曲線。本文對此工況也進行了計算，計算結果與試驗結果對比如圖13所示。試驗中，PBX炸藥發生了起爆，因而應力響應在曲線尾部會有一個突躍，而本文中沒有考慮炸藥的化學反應，計算結果中沒有顯現出這一現象。從對比結果來看，數值計算結果的峰值比試驗結果要高，應力上升速度快于試驗結果，而應力下降速度慢于試驗結果。但從應力響應持續時間及應力響應曲線中出現的平臺特征來看，本文中的模型也顯示出了很好的刻畫能力。造成數值計算結果與試驗結果差異的可能原因有很多，如界面接觸情況等，但由于缺乏更多的試驗細節，本文就不再展開分析。

圖13 應力時程曲線計算與試驗結果的對比Fig.13 Comparison of calculated and experimental stress-time-history curves

4 結論

本文針對復雜應力狀態下PBX炸藥動態響應描述的問題，在K&C模型的基礎上，修改了模型中的損傷演化模式，考慮了材料模量、強度等的壓力相關性，計及了材料的應變率效應，給出了材料拉伸、壓縮子午線之比的取值，處理了結構在損傷破壞過程中變形局部化的問題，提出了拉伸失效單元的處理方法，最終得到的模型較好地描述了撞擊作用下PBX炸藥的損傷破壞形貌及力學響應。

在模型的實現過程中，得出以下結論：

1)材料屈服強度隨壓力呈線性變化時，如果變化的斜率不隨應變率發生改變，就相當于采用徑向應變率增強的方式來描述材料的應變率效應，二者是等價的；PBX炸藥強度隨壓力的變化并不是簡單的線性變化，因而在采用徑向應變率增強方式時，需要做出相應的修正。

2)從拉伸、壓縮子午線之比的數值來看，在較低約束壓力下(壓縮情況下)，K&C模型中經驗方程給出的數值偏大，在單軸壓縮時，拉伸、壓縮子午線之比已經達到0.90，顯然不合理，需要修正；PBX炸藥拉伸、壓縮子午線之比隨壓力的變化還需要進一步的試驗研究。

3)在采用有限元方法刻畫PBX炸藥的變形破壞過程時，需要解決材料變形局部化的問題；數值上解決這個問題的思路主要是抑制計算中不真實的變形局部化的出現，使結構整體響應與試驗結果保持一致；由于只在加載方向進行了校準，雖然可以較好地捕捉到試件加載方向的變形情況，但對于其他未校準方向，就不能很好地描述，當然，變形局部化涉及到對空間不連續問題的描述，這需要進一步的研究。

4)在有限元程序中，單元的失效等同于單元的刪除，而PBX炸藥易于產生拉伸破壞，如果把拉伸破壞簡單的處理為單元失效的話，將會不真實地削弱材料的壓縮承載能力，在程序實現時，要盡量地減小拉伸失效對壓縮過程的影響。