唐 紅, 周俊輝, 呂珂臻, 陳學平, 豆育升

(1. 重慶郵電大學計算機科學與技術學院, 重慶 400065； 2. 中國工程物理研究院化工材料研究所, 四川 綿陽 621999)

1 引 言

高聚物粘結(jié)炸藥(Polymer Bonded Explosive,PBX),在壓制過程中會發(fā)生位移、變形甚至斷裂等各種力學行為, 但是這些力學行為很難直接觀察到,故人們難以明確其成型機理。為表征復合含能材料的壓制成型特性,許多學者通過數(shù)值模擬方法[1-3]計算炸藥不同組分中力學參量的變化,得到壓制過程的力學行為狀態(tài)及模擬數(shù)據(jù),對改進壓制工藝和控制炸藥質(zhì)量具有重要意義。

PBX的力學性能、熱力學性能和化學性質(zhì)都與細觀尺度上的物理和化學過程直接相關[4],因此在細觀尺度上研究PBX的壓制過程對于理解它的力學性能有重大意義。目前,工程上常用的研究PBX力學行為的Lagrangian有限元方法(Finite Element Method,FEM)是基于連續(xù)介質(zhì)的力學方法[5]。如蘭瓊等[2]利用有限元計算方法, 對奧克托今(HMX)基PBX炸藥溫壓時效處理過程的加熱加壓階段進行模擬,得到處理過程中炸藥件尺寸變化規(guī)律,并推導出炸藥件密度變化量,由此預測溫壓時效處理對PBX的作用效果; 張濤等[6]運用基于有限元方法對PBX粉末溫壓成型過程進行了數(shù)值模擬,獲得了粉末體幾何形變、應力場及相對密度分布等相關數(shù)據(jù)等。有限元方法的優(yōu)點是計算速度快,能夠在接近于工程的宏觀尺度上對材料的性質(zhì)進行研究,缺點是不能有效地在計算中考慮微觀結(jié)構(gòu)的性質(zhì),且在處理材料極大變形問題時,網(wǎng)格發(fā)生扭曲、畸變和纏繞,重新劃分網(wǎng)格非常困難,計算精度嚴重下降。近些年來,無網(wǎng)格法被提出并成為研究熱點,常用的無網(wǎng)格法有十幾種 ,其基本思想是將連續(xù)體離散為質(zhì)點單元的形式,在計算過程中, 所有信息都由這些質(zhì)點來表達, 避免了網(wǎng)格重新劃分的難題。作為無網(wǎng)格法之一,物質(zhì)點法(Material Point Method,MPM)本質(zhì)上是一個基于粒子的計算方法[7],相對其他基于網(wǎng)格的數(shù)值方法具有更高的效率和精度,同時在處理大變形問題以及一些帶有接觸的問題時比有限元法具有顯著的優(yōu)勢[8]。

因此,本研究主要應用MPM方法模擬細觀尺度下HMX基PBX壓制成型過程的力學行為,重點分析HMX基PBX壓制成型過程中炸藥顆粒變形、顆粒間的應力傳遞以及溫度變化等細觀力學行為。

2 MPM的基本理論

物質(zhì)點法在計算時滿足質(zhì)量守恒方程和動量守恒方程[7]:

(1)

ρa=·σ+ρb

(2)

式中,ρ為材料密度,kg·m-3;a是加速度,m·s-2;v是質(zhì)點速度,m·s-1;b是單位體積力,N;σ是柯西應力張量,Pa。

將基于背景網(wǎng)格的有限元形函數(shù),試函數(shù)ω帶入動量守恒方程(3)并在區(qū)域Ω上積分,可以得到動量方程的虛功方程[7]:

(3)

式中,dV和dS分別表示微分體積元和面積元;σs是比應力σs=σ/ρ;Γ為指定的應力邊界; 邊界應力為τ; 在指定位移邊界上的試函數(shù)ω=0。

由于把連續(xù)體離散成具有集中質(zhì)量的物質(zhì)點,在整個計算過程中每個物質(zhì)點的質(zhì)量保持不變,故自動滿足質(zhì)量守恒方程。因此, 密度可以表示為[7]:

(4)

式中,mp是物質(zhì)點p的質(zhì)量,g;δ是狄拉克函數(shù);xp是物質(zhì)點p的坐標;Np是物質(zhì)點總數(shù)。

物質(zhì)點數(shù)量和質(zhì)量在整個計算過程中始終不變,質(zhì)量守恒方程自動滿足。在求解動量方程時,質(zhì)點和背景完全固連,隨著背景網(wǎng)格一起運動,因此可通過建立在背景網(wǎng)格節(jié)點上的有限元形函數(shù)Ni(x)實現(xiàn)物質(zhì)點和背景網(wǎng)格節(jié)之間信息的映射。網(wǎng)格節(jié)點和物質(zhì)點這件變量的映射關系為[9]:

(5)

式中,ψp表示物質(zhì)點攜帶的變量;ψi表示網(wǎng)格節(jié)點上的變量;nu表示單元節(jié)點數(shù);xp表示物質(zhì)點的坐標。

最終(3)式可寫成下面的節(jié)點離散形式:

(6)

式中節(jié)點質(zhì)量為:

(7)

內(nèi)力:

(8)

外力:

(9)

(6)式為動量方程在網(wǎng)格節(jié)點上的離散形式,這里采用顯式積分進行求解,具體步驟可以參考文獻[7,9]。

物質(zhì)點法實現(xiàn)了物質(zhì)點到網(wǎng)格和網(wǎng)格到物質(zhì)點的2次映射,物質(zhì)點的應力可利用本構(gòu)方程由式(5)和材料方程更新,得到下一時刻物質(zhì)點所攜帶的變量。物質(zhì)點法的背景網(wǎng)格只在每個時間步內(nèi)和物體固連,在每個時間步結(jié)束時,丟棄已經(jīng)變形的背景網(wǎng)格,這樣避免了傳統(tǒng)的有限元法中背景網(wǎng)格固定不動,造成網(wǎng)格的畸變和纏結(jié)。

3 PBX壓制數(shù)值模擬

3.1 PBX細觀結(jié)構(gòu)模型

圖1a為HMX基PBX炸藥在細觀結(jié)構(gòu)顯微鏡照片[10],圖1a中顯示PBX炸藥顆粒的等效直徑最大為0.2 mm,粘結(jié)劑厚度約為0.01～0.05 mm。其中HMX顆粒大小、形狀各不相同,分布也極其不規(guī)則,HMX顆粒之間還存在粘結(jié)劑。在如此復雜的PBX炸藥細觀結(jié)構(gòu)條件下,建立一種直接反應炸藥細觀結(jié)構(gòu)的模型存在一定的難度。PBX炸藥是由炸藥顆粒和粘結(jié)劑按一定的比例混合壓制成型的,在細觀計算模型的構(gòu)建中需要把PBX炸藥壓制成型過程考慮進去。

a. microstructure ofPBX explosives b. distribution of explosiveparticle

圖1 PBX炸藥細觀結(jié)構(gòu)顯微鏡照片和炸藥顆粒分布

Fig.1 Micrographs of the microstructure of PBX explosives and the distribution of explosive particles

再利用圖像處理技術,將文獻[10]中細觀結(jié)構(gòu)顯微鏡照片(圖1a)構(gòu)建成較合理的PBX炸藥細觀結(jié)構(gòu)模型(圖1b)。在細觀結(jié)構(gòu)計算模型的構(gòu)建中,將雜亂的晶體顆粒近似成圓形,同時用粘結(jié)劑對顆粒外表層進行包裹。在圓形近似的過程中,密集的粒徑分布也被近似成稀疏的分布。這樣的抽象近似對實際力學性能的影響還需進一步研究。

3.2 計算模型及方法

由PBX炸藥顆粒分布模型圖(圖1b),設計基于物質(zhì)點法計算模型圖,如圖2所示。模型由剛性壓頭、剛性模具、炸藥顆粒HMX和粘接劑Estane組成。壓頭的壓制面設為剛性加壓面,模具內(nèi)壁設為剛性約束面。紅色的圓形粒子代表炸藥顆粒,炸藥顆粒外表覆蓋的一層藍色物質(zhì)為Estane粘結(jié)劑。炸藥顆粒之間存在著少許孔隙,孔隙中為空物質(zhì)。模型中共有90個直徑范圍為0.02～0.2 mm的炸藥顆粒,顆粒之間緊密排布,外層為厚度為0.005～0.010 mm的Estane粘結(jié)劑,HMX的體積分數(shù)為85%,粘結(jié)劑體積分數(shù)為6%,其余體積為孔隙。模型中,剛性壓頭以0.3 m·s-1的速度向下壓縮炸藥,炸藥和周圍環(huán)境的初始溫度均為353.15 K,模擬時間設為1 ms。

圖2 計算模型圖

Fig.2 Calculation model

模型中HMX顆粒采用彈塑性材料模型和Grüneisen狀態(tài)方程描述。壓頭和模具材料為剛性。彈塑性材料模型的應力表示為[11]:

(10)

γ0+αμE

(11)

式中,p為壓力,MPa;μ=ρ/ρ0-1=1/Vr-1;c代表常溫常壓無擾動狀態(tài)聲速,m·s-1;ρ0為初始密度,kg·m-3;Vr為相對體積,m3;γ0為Grüneisen常數(shù) ;E為內(nèi)能,J;s1、s2和s3為沖擊波波速-波后粒子速度曲線的斜率系數(shù);α為一階體積修正因數(shù)。HMX顆粒材料模型參數(shù)見表1[12]。

表1 HMX材料模型參數(shù)

Table 1 The model parameters of HMX

density/g·cm-3shearmodulusG/GPaPoisson′sratioyieldstress/MPaspecificheatcapacity/J·kg-1·K-1thermalconductivity/W·K-1·m-11．902．700．3210011500．25

模型中粘結(jié)劑為Estane,采用Prony 級數(shù)形式的粘彈性本構(gòu)模型和Tait狀態(tài)方程 ,Prony 級數(shù)形式粘彈性模型的應力表示為[13]:

(12)

式中,σ為應力,MPa;G(t)為剪切松弛核函數(shù),e為應變偏量部分(剪切變形),Δ為應變體積部分(體積變形),t為當前時間,s;τ為過去時間,s;I為單位張量。用Prony級數(shù)表示粘彈性屬性的基本形式為[11]:

G(t)=G

(13)

K(t)=K

(14)

式中,G和Gi是剪切模量,MPa;K和Ki是體積模量,和是各Prony級數(shù)分量的松弛時間(Relative time),s。Tait狀態(tài)方程為[13]:

(15)

B(T)=CK(0,T)

(16)

式中,p為壓力,MPa;T為溫度,K;C=0.0894是通用的Tait常數(shù);V(0,T) 是在0壓下由溫度決定的體積函數(shù);K(0,T) 是在0壓下溫度決定的體積模量函數(shù)。粘接劑Estane材料模型參數(shù)見表2[12]。

表2 Estane材料模型參數(shù)

Table 2 The model parameters of Estane

bulkmodulus/GPashearmodulus/GPadensity/g·cm-3specificheatcapacity/J·kg-1·K-1thermalconductivity/W·K-1·m-13．60．271．1011500．25

4 模擬結(jié)果及分析

4.1 PBX壓制過程中的變形

剛性壓頭以0.3 m·s-1的速度從上往下壓縮時,PBX炸藥體系在不同時刻的變形情況如圖3所示。數(shù)值模擬結(jié)果表明,PBX在壓縮的過程中,其顆粒形狀是不斷變化的。圖3a和圖3b為壓縮過程中的整合階段,炸藥顆粒重新排列,形成接觸擠壓。在圖3a中,位于接觸面的炸藥顆粒首先受到剛性壓頭的擠壓,炸藥顆粒向下的運動,顆粒間的縫隙不斷減小。在0.50 ms時,炸藥顆粒已趨于壓實的狀態(tài),各顆粒之間接觸擠壓已經(jīng)基本形成。圖3c和圖3d表示炸藥顆粒的鞏固階段發(fā)生塑性形變的過程。在圖3c中炸藥體系壓縮到0.75 ms時,炸藥顆粒之間均受到不同程度力的擠壓,發(fā)生塑性變形,顆粒輪廓明顯改變。一直壓縮到圖3d的塑性變形的狀態(tài)。炸藥顆粒的特征基本消失,形成密實整體,壓制過程變?yōu)閷φㄋ幷w的壓縮作用。

a. 0.25 ms b. 0.50 ms

c. 0.75 ms d. 1 ms

圖3 PBX炸藥壓制過程中不同時刻的變形模擬結(jié)果

Fig.3 Deformation simulation results of PBX at different times during compression

4.2 PBX壓制過程中的應力變化

圖4為PBX炸藥壓制過程中不同時刻應力分布模擬結(jié)果。從圖4a和圖4b看出,在0.25 ms時,應力沿炸藥顆粒間的接觸面向剛性模具底板傳播,并形成了多條的應力鏈。同一個PBX炸藥內(nèi)部應力分布不均勻,炸藥體系頂部處于高應力狀態(tài), 底部部分應力較小。應力集中點主要集中在顆粒與剛性壓頭的接觸部分,PBX炸藥內(nèi)部存在應力梯度,各炸藥顆粒之間應力差較大。圖4c和圖4d可以看出炸藥體系應力梯度在減小,這是由于大部分炸藥顆粒已發(fā)生滑移和變形, 并且處于塑性變形狀態(tài)。在壓縮至1 ms時,炸藥體系已基本處于密實狀態(tài),但各部分應力分布仍然呈不均勻狀態(tài)。由于在炸藥壓制成型過程中,炸藥顆粒間所受的壓力和摩擦力略大于其他部位,而導致應力在顆粒接觸面較大。從圖4可以看出,早期階段炸藥體系呈現(xiàn)出一個比較松散的狀態(tài),在不斷壓縮過程中,各顆粒之間受到壓力開始接觸擠壓而導致應力變化。這與劉群[14]等利用非線性有限元計算方法模擬JO-9159炸藥顆粒壓制成型過程得到的結(jié)果“應力集中主要出現(xiàn)在藥粒與約束面的接觸部分;當藥粒之間空隙被填滿, 藥粒進入塑性變形后, 藥粒內(nèi)部壓力迅速升高且壓力趨于一致”類似。

a. 0.05 ms b. 0.25 ms

c. 0.75 ms d. 1 ms

圖4 PBX炸藥壓制過程中不同時刻應力分布模擬結(jié)果

Fig.4 Stress distribution simulation results of PBX at different times during compression

圖5為模擬壓制過程中PBX炸藥體系整體應力-應變率曲線。圖5表明,在壓縮初始階段,應力上升速度較快,應力和應變呈線性關系,PBX 材料表現(xiàn)為彈性性質(zhì)。超過材料的屈服力后,受到粘結(jié)劑性質(zhì)影響,PBX材料開始表現(xiàn)為粘性性質(zhì),應力上升速度減緩。上述結(jié)論與蔡宣明[15]等關于動態(tài)力學響應及細觀損傷破壞模式的研究結(jié)果“彈性階段,應力-應變曲線開始階段具有較好的線性關系; 屈服階段,隨著應變的增加,應力開始減小; 緊接著進入強化階段,PBX 模擬材料抗壓能力又開始增強,應力隨著應變的增大而增大”相符。

圖5 PBX炸藥體系壓縮的應力-應變率曲線

Fig.5 Compressive stress-strain rate curve for PBX system

4.3 PBX壓制過程中的溫度變化

圖6為PBX炸藥體系內(nèi)不同時刻溫度分布模擬結(jié)果。圖6顯示,炸藥體系受壓后,剛性壓頭與PBX炸藥體系接觸面溫度首先升高,隨著壓制的進行,炸藥內(nèi)部溫度逐漸升高,而達到平衡后,炸藥內(nèi)部存在溫差。從圖6a和圖6b看出,炸藥體系在0.10 ms時,炸藥體系頂部在剛性壓頭擠壓和顆粒間摩擦的作用下溫度逐漸升高; 0.25 ms時,體系內(nèi)部的局部溫度最大上升了10 K。炸藥體系底部由于受到壓力影響較小,溫度明顯低于頂部。圖6c和圖6d可以看出溫度變化已經(jīng)趨于平衡,表明體系內(nèi)部溫差逐漸減小。在1 ms時,炸藥體系已接近密實狀態(tài),除中心處附近溫度較高以外,其余顆粒溫度基本相同,體系內(nèi)部最大溫差在20 K左右。圖7為模擬過程中,炸藥體系整體溫度變化曲線。從圖6、圖7得出,模擬結(jié)束時,體系整體溫度上升到了359.12 K,局部最高溫度為385.74 K。溫度變化的模擬結(jié)果與文獻[16]的結(jié)果“PBX炸藥在升壓過程中均出現(xiàn)內(nèi)部溫度瞬時增高的現(xiàn)象,達到平衡后,炸藥內(nèi)部存在溫差,藥柱中心處藥粒溫度最高的結(jié)果,壓藥柱內(nèi)部最大溫升為29 K”基本一致。

a. 0.10 ms b. 0.20 ms

c. 0.75 ms d. 1 ms

圖6 PBX炸藥壓制過程中不同時刻溫度分布模擬結(jié)果

Fig.6 Temperature distribution simulation results of PBX at different times during compression

圖7 PBX炸藥體系壓縮的溫度時間變化曲線

Fig.7 Compressive temperature-time curve for PBX system

5 結(jié) 論

采用物質(zhì)點法,對HMX基PBX炸藥壓制過程中的力學行為進行了細觀尺度的模擬。在模擬壓縮過程中,剛性壓頭以0.3 m·s-1的速度從上往下壓縮,模擬時間持續(xù)1 ms。模擬結(jié)果表明,PBX在壓縮的過程中,位于接觸面的炸藥顆粒首先受到剛性壓頭的擠壓,炸藥顆粒向下的運動,顆粒間的縫隙不斷減小; 在0.50 ms時,炸藥顆粒已趨于壓實的狀態(tài); 在0.75 ms時,炸藥顆粒之間均受到不同程度力的擠壓,發(fā)生塑性變形,顆粒輪廓明顯改變; 至1 ms時PBX無顯著的顆粒特征,形成密實整體。

結(jié)果表明,壓縮中炸藥顆粒存在整合和鞏固兩個階段。在整合階段(即0～0.5 ms),PBX顆粒受到壓力而重新排列,應力和溫度集中主要出現(xiàn)在炸藥顆粒與壓縮面的接觸部分,并形成了多條應力鏈。應力鏈沿炸藥顆粒間的接觸面向上和向下傳播。炸藥體系頂部應力、溫度變化均大于底部部分; 當炸藥顆粒之間空隙被填滿,炸藥體系進入鞏固階段(即0.5～1 ms)。該階段發(fā)生塑性形變,炸藥體系內(nèi)部炸藥體系應力梯度在減小,炸藥顆粒間所受的壓力和摩擦力略大于其他部位,而導致應力在顆粒接觸面較大; 在壓制過程中,炸藥體系溫度升高,局部溫度最大上升了10 K。最后,炸藥體系達到平衡,內(nèi)部存在溫差,體系內(nèi)部最大溫差在20 K左右。除中心處附近溫度較高以外,其余顆粒溫度基本相同。

本研究中的計算模型和模擬方法同樣適用于其他組成的PBX炸藥,在對其它組分的PBX炸藥進行模擬時,只需要輸入其相關的物理、化學性質(zhì)參數(shù)。由于計算程序的限制,本研究中只對HMX基PBX壓制過程進行了二維下的數(shù)值模擬,下一步工作是建立HMX基PBX炸藥細觀結(jié)構(gòu)的三維計算模型,對炸藥壓制過程進行三維數(shù)值模擬計算。不同壓力下、不同顆粒尺徑下等多種情況下對應力的傳播、變形的影響,以及炸藥顆粒的損傷本構(gòu)模型的研究有待進一步研究。

參考文獻:

[1] Banerjee B, Cady C M, Adams D O. Micromechanics simulations of glass-estane mock polymer bonded explosives[J].ModellingandSimulationinMaterialsScienceandEngineering, 2003, 11(4): 457.

[2] 蘭瓊, 唐維 ,賀建華, 等. 某HMX基PBX溫壓時效處理過程變形規(guī)律數(shù)值模擬[J]. 含能材料, 2015, 23(6): 543-547.

LAN Qiong, TANG Wei, HE Jian-hua, et al. Simulation on deformation of HMX based PBX by thermal-pressure treatment[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2015, 23(6): 543-547.

[3] Lusti H R, Gusev A A. Finite element predictions for the thermoelastic properties of nanotube reinforced polymers[J].ModellingandSimulationinMaterialsScienceandEngineering, 2004, 12(3): S107.

[4] Menikoff R, Shaw M S. Review of the forest fire model for high explosives[J].CombustionTheoryandModelling, 2008, 12(3): 569-604.

[5] SaeedMoaveni. Finite element analysis theory and application with ANSYS[M]. Beijing: Publishing House of Electronics Industry, 2012.

[6] 張濤, 趙北君, 朱世富, 等. PBX 粉末成形的數(shù)值模擬研究[J]. 材料工程, 2009, 35(5): 68-72.

ZHANG Tao, ZHAO Bei-jun, ZHU Shi-fu, et al. Numerical simulation study of PBX powder-forming[J].JournalofMaterialsEngineering, 2009, 35(5): 68-72.

[7] Sulsky D, Chen Z, Schreyer H L. A particle method for history-dependent materials[J].ComputerMethodsinAppliedMechanicsandEngineering, 1994, 118(1): 179-196.

[8] Hornemann U, Holzwarth A. Shaped charge penetration in alumina targets[J].InternationalJournalofImpactEngineering, 1997, 20(1): 375-386.

[9] Andersen S, Andersen L. Analysis of spatial interpolation in the material-point method[J].Computers&Structures, 2010, 88(7): 506-518.

[10] Peterson J R, Wight C A,Berzins M. Applying high-performance computing to petascale explosive simulations[J].ProcediaComputerScience, 2013, 18: 2259-2268.

[11] 王晨, 陳朗, 劉群, 等. 多組分PBX炸藥細觀結(jié)構(gòu)沖擊點火數(shù)值模擬[J]. 爆炸與沖擊, 2014, 34(2): 167-173.

WANG Chen, CHEN Lang, LIU Qun, et al. Numerical simulation for analyzing shock to ignition of PBXs with different compositions in meso-structural level[J].Explosion&ShockWaves, 2014, 34(2):167-173.

[12] 董海山. 周芬芬. 高能炸藥及相關物性能[M]. 北京: 科學出版社,1989: 334-340.

DONG Hai-shan, ZHOU Fen-fen. High energy explosives and correlative physical properties[M]. Beijing: Science Press, 1989: 334-340.

[13] Xue L, Borodin O, Smith G D, et al. Micromechanics simulations of the viscoelastic properties of highly filled composites by the material point method (MPM)[J].ModellingandSimulationinMaterialsScienceandEngineering, 2006, 14(4): 703.

[14] 劉群, 陳朗, 魯建英, 等. 炸藥顆粒壓制成型數(shù)值模擬[J]. 高壓物理學報, 2009, 23(6): 421-426.

LIU Qun, CHEN Lang, LU Jian-ying, et al. Numerical simulation of explosive particles compaction[J].ChineseJournalofHighPressurePhysics, 2009, 23(6): 421-426.

[15] 蔡宣明, 張偉, 魏剛, 等. PBX 模擬材料動態(tài)力學響應及細觀損傷模式[J]. 含能材料, 2014, 22(5): 658-663.

CAI Xuan-ming, ZHANG Wei, WEI Gang, et al. Dynamic mechanics response and mesoscopic damage of a PBX simulant[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2014, 22(5): 658-663.

[16] 賈憲振, 王曉峰, 陳松, 等. 炸藥等靜壓工藝安全性的數(shù)值模擬研究[J]. 兵工自動化, 2014,33(7): 65-67

JIA Xian-zhen, WANG Xiao-feng, CHEN Song, et al. Numerical study of safety of explosive isostatic pressing[J].OrdnanceIndustryAutomation, 2014, 33(7): 65-67.