郭支明，秦棟澤，賈憲振

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

0 引言

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

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

1 基本理論

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

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

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

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

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

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

(1)

(2)

1.2 守恒型控制方程

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

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

(3)

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

(4)

可得到動(dòng)力動(dòng)量和能量方程：

(5)

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

1.3 人工粘性算子

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

(6)

(7)

(8)

這里使用普通的RK梯度算子來(lái)計(jì)算該速度梯度：

(9)

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

2 損傷演化模型

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

(10)

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

對(duì)式(10)積分，有

(11)

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

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

表1 PBX9501材料參數(shù)

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

σ=(1-D)Eε

(12)

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

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

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

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

3.1 仿真結(jié)果

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

表2 PBX9501試件不同時(shí)刻位移云圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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