基于CRKSPH的炸藥巴西壓裂仿真研究①

郭支明，秦棟澤，賈憲振

(1.中北大學，太原 030051；2.西安近代化學研究所，西安 710065)

0 引言

PBX9501炸藥(高聚物粘結炸藥)在外載荷作用下內部微缺陷會形核長大，并形成微裂紋進行演化，形成內部損傷，使其力學性能劣化。炸藥的損傷斷裂性能實驗主要有直接拉伸實驗、巴西實驗、三點彎曲實驗等。巴西實驗的試樣制備簡單，而且方便快捷[1]。由于PBX9501炸藥屬于粘彈性-準脆性材料，準確描述其力學響應及斷裂特性比較困難。國內一些學者對此方向進行了研究，趙四海[2]等利用VISCO-SCRAM 模型，計算了巴西壓裂實驗中PBX9501炸藥圓盤試樣的力學響應，并和試驗數據進行了對比。陳鵬萬等[3]采用巴西實驗研究了PBX9501炸藥的斷裂性能，并通過顯微鏡觀察了炸藥顆粒斷裂、界面脫粘、粘結劑撕裂等破壞形式。龐海燕[4]等采用標準巴西實驗、圓弧巴西實驗和橡膠墊巴西實驗測量得到PBX炸藥的拉伸強度。在數值計算方面，崔云霄[5]等基于EFG方法和線性內聚力模型，模擬了準靜態加載圓弧巴西實驗中PBX9501炸藥的斷裂行為。戴開達等[6]采用擴展有限元方法，對不同種類巴西實驗中PBX炸藥的裂紋起裂、擴展和斷裂行為進行了數值模擬。顏學堅[7]基于液壓致裂原理，自主搭建了測試平臺,發展了拉伸強度液壓致裂測試方法，實現了PBX炸藥的拉伸強度的準確測試。袁洪魏等[8]基于Boltzmann本構模型，建立了考慮靜水壓力影響的TATB基PBX準靜態Boltzmann-P非線性彈性本構模型。李尚昆等[9]從材料的力學行為特性、實驗方法、本構模型和強度理論四個方面，對高聚物粘結炸藥(PBX)的力學性能特征進行了歸納和評述。蒙君煚等[10]用分子動力學方法，模擬了功能助劑對界面結合能的影響規律，試驗結果與數值模擬結果吻合。原曾年等[11]采用CT技術和數值仿真的方式，研究了PBX材料樣品內的缺陷和損傷行為，模擬結果與試驗結果吻合較好。 性能測試方面，溫茂萍等[12]采用圓弧壓頭與差動變壓器引伸計相結合的巴西試驗方法，測試了脆性炸藥力學性能。陳科全等[13]以壓裝HMX基PBX炸藥為對象,采用改進后的圓弧巴西試驗,研究了不同預制缺陷條件下炸藥的破壞過程,并與無缺陷炸藥的試驗結果進行對比。劉晨等[14]設計了一種帶中心圓孔的平板試樣進行實驗，基于數字圖像相關方法(DICM)，對實驗進行應變場分析。

綜上所述，炸藥巴西壓裂過程的數值計算方法采用有限元或EFG方法(有背景網格)等基于控制方程弱形式的離散格式，這些基于網格的方法在計算固體斷裂等強間斷問題時，受到網格連通性的限制。本文使用基于控制方程強形式的離散格式——CRKSPH(守恒型再生核粒子方法)，對炸藥巴西壓裂過程進行模擬。這種純拉格朗日方法不再依賴于網格，在求解斷裂和裂紋擴展問題上有獨特的優勢。

1 基本理論

CRKSPH方法是SPH(Smoothed particle hydrodynamics，光滑粒子流體動力學)方法的一種改進。SPH是一種拉格朗日方法，它基于粒子近似對控制介質運動的積分或偏微分方程進行離散，連續場變量的值通過將離散粒子的性質(如質量，動量和能量)與內插核(通常用W表示)進行卷積來表示，是一種完全拉格朗日無網格數值方法，物質界面清晰，材料的變形不依賴于網格，而通過粒子的運動來描述。因此，SPH既具有拉氏計算中描述物質界面準確的優勢，又兼備無網格方法的長處，適宜計算帶有流體大變形及運動邊界的各類問題，如材料損傷斷裂、爆炸爆轟、沖擊侵徹等[15-22]。SPH方法最早由LUCY[23]和GINGOLD[24]等于1977年提出，并應用于天體物理學，SWEGLE 等[25]將SPH方法用于模擬爆炸問題。

盡管SPH方法已成功應用于許多領域，但傳統的SPH方法仍存在嚴重的缺陷。其中，最嚴重的是SPH缺乏零階一致性，即所謂的“E0級誤差”[26-27]。SPH的另一個常見問題是精確捕捉沖擊波陣面所需的人工粘性公式。 MONAGHAN和GINGOLD[28]的標準粘性計算公式因阻尼過大，使結果嚴重失真。目前為止，人工粘性計算仍是一個尚未解決的問題。

解決零階精度問題一個有意義的嘗試是再生核(RK)方法[18,29]， RK方法以附加的自由度對SPH插值內核進行了改進，可以精確地再生任意階的多項式。這完全消除了SPH的零階誤差，但這種方法使粒子之間的內核不對稱。本文使用CRKSPH(守恒型再生核粒子方法)[29]，其使用線性再生核，基于RK插值重新推導了動量守恒方程。這種方法可以精確地保持線性動量守恒，以達到一定精度，可有效消除零階插值誤差和過度使用人工粘度。

CRKSPH保留了傳統SPH方法的許多優點，同時對SPH的缺點進行了改進，特別是人工粘性過高和零階精度誤差。為此，CRKSPH方法對傳統SPH方法進行三個變化：線性再生核函數、守恒型的控制方程、精確的人工粘性算子。

1.1 線性再生核函數

為解決SPH無法重構所需階數的多項式，在傳統SPH插值核函數中添加一些修正項，以精確再生常量、線性或更高階數的項。這種插值方法，稱為再生核方法(RPKM)，再生核公式可以推廣到任意階。 RK插值公式和變量場的梯度表示為

(1)

(2)

1.2 守恒型控制方程

關于再生核插值和梯度，在建立守恒型控制方程時，CRKSPH方法利用了MLSPH構造中的一些派生方法。

假設插值使用通用的核函數ψ(x)，其中守恒密度U和通量F定義為

(3)

通過插值函數和體積積分：

(4)

可得到動力動量和能量方程：

(5)

其中的力對是反對稱的。因此，這樣推導出來的動量方程，再加上再生核形式，可以保證精確的線性動量守恒。

1.3 人工粘性算子

這里使用經典的van Leer限制器[30-31]作為新算子的基礎，這個算子是基于網格方法的計算流體力學理論。修正后的算子如下：

(6)

(7)

(8)

這里使用普通的RK梯度算子來計算該速度梯度：

(9)

等式(6)的第一部分，就是普通的van Leer限制器；當ηij<ηcrit時激活，當點靠近時將限制器強制為零。

2 損傷演化模型

PBX9501材料在拉伸過程中，隨著變形的增加,損傷演化可表示為應變的函數。對于PBX9501材料，可建立如下損傷演化方程[32]，損傷參數D(ε)演化根據公式：

(10)

式中B和n為材料常數；εc為材料的斷裂應變；ε0為損傷開始發展時的臨界應變。

對式(10)積分，有

(11)

式中D0為材料的初始損傷。

PBX9501材料的初始損傷D0可通過材料的實際密度與理論密度的偏差來初步估計。根據文獻[30-31]實驗結果，可確定損傷演化方程中的材料常數，如表1所示。

表1 PBX9501材料參數

從而可建立材料的損傷本構關系：

σ=(1-D)Eε

(12)

式中σ為材料應力張量；E為彈性模量矩陣；ε為應變張量。

3 數值模型和仿真結果

計算時，將三維巴西實驗問題簡化為平面應力問題，建立的CRKSPH模型如圖1左圖所示，圓盤試樣直徑為20 mm，共劃分14 400個粒子。為了模擬真實的實驗加載，圓盤下部放置在剛性面上，上部施加位移載荷0.05 mm/min。圖1右圖為標準巴西壓裂試驗基本破壞形式示意圖。

圖1 CRKSPH計算粒子模型(左)和標準巴西壓裂試驗基本破壞形式 (右)

3.1 仿真結果

以上文所述加載工況(加載速率0.05 mm/min)和粒子模型進行巴西圓盤仿真實驗，位移演化過程如表2所示，損傷及裂紋演化情況如表3所示。

表2 PBX9501試件不同時刻位移云圖

在整個仿真試驗過程中，試件沿載荷方向處于壓縮狀態，沿水平方向處于拉伸狀態；在載荷作用點(接觸點)附件一定區域內會發生剪切破壞，最終在試件內部沿著載荷作用方向會形成一條貫穿裂紋，將試件一分為二，計算結果可充分地展現這一過程。

從表2和表3中的結果可知：

(1)t=0.36 s圓盤中心位置微缺陷發展為局部損傷，在隨后的加載過程中逐步形成裂紋，并沿著縱向向上和向下擴展；

(2)在圓盤的加載點兩側，由于局部剪切作用產生局部破壞；

(3)在試件的底部附件區域由于剪切和拉伸綜合作用也發展出了多處不同程度的損傷裂紋，損傷情況基本以圓盤縱向對稱軸呈對稱分布。

將表2中位移云圖和圖1右圖中基本破壞形式示意圖做對比，可看出仿真結果和巴西壓裂基本破壞模式保持一致。圖2所示為0.5 s時微觀缺陷分布、損傷參數D和文獻中試驗結果[6]。

圖2 0.5 s時微觀缺陷分布(左)、損傷參數D和文獻中試驗結果(中和右)

從圖2斷裂情況和文獻中試驗結果[6]進行對比可知，試件損傷、裂紋的產生和擴展方向均和試驗基本吻合，而且基于CRKSPH方法仿真，除了可追蹤明顯的裂紋外，還可對產生損傷的位置進行預測和追蹤，并且可有效處理多裂紋的產生和合并。

從圖3所示試件頂部裂紋分布和試驗數據對比情況可看出，仿真結果中頂部3條主裂紋分布和試驗結果保持一致[4]。圖4所示為基于本構關系(8)計算得到的試件中心位置應力-應變曲線。

圖3 試件頂部裂紋分布情況對比 (左為試驗結果和右為仿真結果)

圖4 計算得到的應力位移曲線

從圖4中可看出，應變在0～0.002 6這個階段，應力基本呈線性增大，材料處于彈性變形階段；當應變大于0.002 6時，應力達到最大值5.87 MPa，此時試件中心的拉應力超過抗拉強度，從而產生初始裂紋；裂紋繼續擴展，導致試件中應力減小，當應力減小到2.857 MPa時，形成貫穿裂紋。

3.2 不同加載速率對試件斷裂形式的影響

由于炸藥在運輸和使用過程中會遇到各種載荷情況，有必要考慮加載速率對斷裂形式的影響。本節基于上文的計算方法和模型，考慮不同加載速率下試件的斷裂形態和規律。

工況和計算結果如表4所示。這里主要考慮的加載速率為0.03、0.05、0.07 mm/min。

表4 三種加載速率下相同中心位置位移對應的斷裂形態

從計算結果可看出：

(1)加載速率不同，試件中裂紋起始位置不同(如表4中第一行各圖中的黃色圓圈標記處)。加載速率為0.03 mm/min時，裂紋起始位置位于試件中心點下方2.4 mm處，加載速率為0.05 mm/min時，裂紋起始位置位于中心點，加載速率為0.07 mm/min時，裂紋起始位置有2個，一個位于中心點上方2.5 mm處，另一個位于中心點下方2.5 mm處。

(2)試件中心點位移相同的情況下，各加載速率對應的斷裂形態不同。加載速率越大，壓裂所產生的裂紋數量越多，而且在試件外表面，裂紋數量從底部往上沿外邊緣快速增加。

4 結論

(1)采用 CRKSPH方法，對PBX9501炸藥標準巴西圓盤壓裂過程計算結果中，試件破壞方式和裂紋走向與文獻中實驗結果基本保持一致。

(2)不同加載速率對試件中裂紋的產生和擴展都有一定的影響，隨著加載速率的增加，試件中裂紋數量迅速增加，碎裂程度增大。

(3)計算結果表明，CRKSPH和材料的損傷本構相結合，能夠有效求解固體炸藥損傷、斷裂和裂紋擴展問題，可有效處理彈脆性材料中多裂紋產生、合并等物理現象。計算方法和結果對炸藥材料工程設計有一定的指導意義。

本文主要以探索新算法在炸藥材料破壞仿真方面的應用，試驗數據和圖片均取自相關文獻，下一步研究中，將會結合試驗數據對計算模型進行進一步修正。