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

郭支明，秦棟澤，賈憲振

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

0 引言

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

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

1 基本理論

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

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

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

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

1.1 線性再生核函數(shù)

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

(1)

(2)

1.2 守恒型控制方程

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

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

(3)

通過插值函數(shù)和體積積分：

(4)

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

(5)

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

1.3 人工粘性算子

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

(6)

(7)

(8)

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

(9)

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

2 損傷演化模型

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

(10)

式中B和n為材料常數(shù)；εc為材料的斷裂應(yīng)變；ε0為損傷開始發(fā)展時的臨界應(yīng)變。

對式(10)積分，有

(11)

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

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

表1 PBX9501材料參數(shù)

從而可建立材料的損傷本構(gòu)關(guān)系：

σ=(1-D)Eε

(12)

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

3 數(shù)值模型和仿真結(jié)果

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

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

3.1 仿真結(jié)果

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

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

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

從表2和表3中的結(jié)果可知：

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

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

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

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

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

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

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

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

圖4 計算得到的應(yīng)力位移曲線

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

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

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

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

表4 三種加載速率下相同中心位置位移對應(yīng)的斷裂形態(tài)

從計算結(jié)果可看出：

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

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

4 結(jié)論

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

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

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

本文主要以探索新算法在炸藥材料破壞仿真方面的應(yīng)用，試驗數(shù)據(jù)和圖片均取自相關(guān)文獻(xiàn)，下一步研究中，將會結(jié)合試驗數(shù)據(jù)對計算模型進(jìn)行進(jìn)一步修正。