王竟成, 羅景潤

(中國工程物理研究院總體工程研究所, 四川 綿陽 621999)

1 引 言

高聚物粘結炸藥(Polymer Bonded Explosive, PBX)是一種顆粒高度填充復合材料,具有復雜的微細觀結構。顆粒的性能、形貌、尺寸分布以及粘結劑的性質都對PBX的宏觀性能有著密切的影響。細觀力學逐漸發(fā)展了多種理論方法預測復合材料的有效性能[2-3],如上下界限法、稀疏法、Mori-Tanaka法、自洽法、微分法等,這些方法應用于PBX存在以下問題:

(1)PBX的顆粒含量很高,而傳統細觀力學方法預測顆粒含量較高的復合材料有效性能會出現很大偏差;

(2)PBX組分中顆粒與基體的彈性模量相差3～4個數量級,以致這些方法預測的結果具有量級上的差異;

(3)PBX炸藥顆粒是典型的不規(guī)則多面體,而細觀理論方法通常假設夾雜顆粒為球形或橢球形,難以很好地刻畫PBX的細觀結構特征。

由于上述細觀力學理論方法的局限性,構建合適的細觀數值模型,采用有限元法計算PBX的有效性能是一種有效途徑。目前,PBX的細觀數值模型主要有: 圓形顆粒規(guī)則排布和隨機排布[6-9]、六邊形顆粒規(guī)則排布、正方形顆粒模型?；谶@些模型,采用有限元方法,可以實現對PBX有效性能的預測,但是,這些數值模型并未反映PBX炸藥顆粒的不規(guī)則多面體特征,其預測結果與實測值差異較大。本研究采用Voronoi法以種子點為中心生成隨機多邊形顆粒,顆粒的平均尺寸通過種子數進行控制,建立了Voronoi細觀數值模型,力圖較真實反映PBX細觀顆粒特征,實現炸藥顆粒的高度填充。

2 數值模擬

2.1 細觀模型的建立

為考察顆粒形貌和排列分布方式對PBX有效模量的影響,首先建立四種構型的二維RVE模型: circle_r(圖1a)為圓形顆粒隨機排布模型,circle_1(圖1b)為圓形顆粒規(guī)則排布模型1,circle_2(圖1c)為圓形顆粒規(guī)則排布模型2,hexagon(圖1d)為六邊形顆粒規(guī)則排布模型。其中,圖1a,圖1b,圖1c為圓形顆粒的不同分布形式,顆粒含量均為49%,用于分析顆粒排列分布對有效模量的影響; 圖1d與圖1c的顆粒含量、分布形式相同,用于分析顆粒形貌的影響。

圖2a為超聲波下PBX的細觀形貌: 炸藥顆粒為不規(guī)則多面體,含量極高。為接近PBX的細觀形貌,實現炸藥顆粒的高度填充,更有效地模擬PBX的有效力學性能,建立不規(guī)則多邊形顆粒的Voronoi模型。在一定區(qū)域內播撒隨機種子點,采用Voronoi法[10]以種子點為中心生成隨機多邊形顆粒。顆粒的平均尺寸通過種子數進行控制。為保證顆粒生成的形態(tài)質量,有必要篩選種子點。圖2b為具有周期性結構的Voronoi圖,炸藥顆粒的平均粒徑為100 μm。進一步生成粘結劑層后,將節(jié)點信息導入ANSYS獲得有限元模型(圖2c),它與PBX真實的細觀形貌(圖2a)具有較好的相似性。將種子點固定,改變粘結劑的厚度構建不同顆粒含量的Voronoi模型,其顆粒形貌和分布形式完全相同,有利于研究顆粒含量對PBX有效性能的影響。

a. circle_rb. circle_1c. circle_2d. hexagon

圖1不同排布或形貌的代表性體積單元(RVE)模型

Fig.1RVE models with different particle distributions or shapes

a. digital image[11]b. Voronoi imagec. Voronoi model

圖2超聲波下PBX細觀形貌和Voronoi模型

Fig.2Digital image under ultrasonic and Voronoi model for PBX

計算中炸藥顆粒選用奧克托今(HMX)晶體,粘結劑選用聚酯型聚氨酯(Estane),均采用各向同性的彈性本構模型,其理想材料參數見表1[12]。從表1中可見基體Estane和炸藥顆粒HMX之間的彈性性能有巨大的差異,泊松比相差一倍,楊氏模量相差四個數量級。

表1各組分彈性參量

Table1The elastic properties of ingredients

ingredientE/MPaνmatrix：Estane0．770．499particles：HMX2．5325×1040．25

Note:Eis the Young′s modulus,νis the Poisson′s ratio.

2.2 RVE尺寸的選擇

代表性體積單元(Representative Volume Element, RVE)是細觀力學中的一個基本概念,需要滿足尺度的二重性: 宏觀上可以看成一個物質點,RVE中的宏觀應力、應變場可視為均勻的; 細觀上包含了足量的細觀結構信息,可以代表局部連續(xù)介質的統計平均。細觀應力應變場只通過它們的體積平均值對材料的宏觀性能產生影響[13]。不同材料選取的RVE尺寸有很大的差別,一般來說需要滿足:Lc?LRVE?LM,Lc為材料組分特征尺寸,LRVE為RVE尺寸,LM為結構宏觀特征尺寸。只有當LRVE/Lc足夠大時,非均勻材料組分結構的統計平均性質才能趨于穩(wěn)定,細觀結構的不確定影響才能忽略不計。

適當尺寸的RVE,其有效彈性性質不受邊界條件類型的影響[14-15]。Lemaitre[16]曾對各種典型材料代表性單元的最小尺寸進行了粗略估計,金屬: 0.1 mm×0.1 mm×0.1 mm,混凝土: 100 mm×100 mm×100 mm,木材: 10 mm×10 mm×10 mm,高分子和顆粒復合材料: 1 mm×1 mm×1 mm。Ostoja-Starzewski[17]研究片狀或針狀顆粒復合材料的有效剛度后發(fā)現,在最大許可誤差為5%和顆粒剛度是基體剛度100倍的情況下,RVE最小尺寸是顆粒特征尺寸的10～20倍。真實炸藥顆粒的平均尺寸約100 μm,因此選擇PBX二維RVE模型的尺寸為1000 μm×1000 μm。

2.3 有限元計算方法

采用周期性的RVE模型(圖3)來模擬均勻化介質的力學行為,平行于邊界面的邊界節(jié)點位移通過以下關系耦合起來:

ui-uj=0

(1)

式中，ui,uj分別為邊界節(jié)點i和j平行于邊界面的位移,μm。采用易于實現的均勻應變邊界條件進行有效模量的計算,即左邊界和下邊界分別約束住x和y向位移,上邊界和右邊界進行耦合約束,上邊界通過施加位移載荷來模擬單軸壓縮試驗,見圖4。有限元模擬采用平面應變假設和平面八節(jié)點單元plane183。

圖3周期性的Voronoi模型

Fig.3Periodic Voronoi models

圖4邊界條件

Fig.4Boundary conditions

根據平面應變假設和邊界條件,由:

(2)

可得:

(3)

(4)

式中,Ux為右邊界節(jié)點x向位移,μm;Uy為上邊界節(jié)點y向位移,μm。

根據文獻[18],平均應力既可以使用體積分得到,也可以通過面積分得到:

(5)

式中，F代表作用反力,N為RVE邊界上節(jié)點的個數。對于當前二維問題,可簡化為:

(6)

式中，Fre為ANSYS后處理中提取的上邊界節(jié)點反力之和,N。利用式(3)、(4)、(6)有:

(7)

進一步求解出體積模量K和剪切模量G:

(8)

式中,Uy為采用均勻應變邊界條件時,在上邊界節(jié)點施加的位移載荷,μm;Ux為后處理中提取的右邊界節(jié)點位移Ux,μm;Fre為上邊界節(jié)點反力之和,N;ν為泊松比,E為楊氏模量,Pa;K為體積模量,Pa;G為剪切模量,Pa。

3 結果與討論

表2列舉了圓形顆粒不同分布形式之間的有限元結果,用于分析顆粒分布形式對有效模量的影響。circle_r與circle_1各項差異都較小; 在相同顆粒含量下,不同分布形式之間泊松比ν和體積模量K的變化相對較小,而楊氏模量E和剪切模量G的差異卻較大,說明E和G對顆粒分布形式更為敏感。以circle_2為標準,c=49%時,circle_1中E高出76.97%,G高出78.22%;c=64%時,E和G分別高出158.96%和160.56%?？梢?這兩種排布方式隨著顆粒含量的增加,差異在進一步增大。以上分析說明,顆粒的排布方式對PBX的有效模量有顯著影響,且其影響隨著顆粒含量的上升有增大的趨勢; 顆粒含量很高時,使用炸藥顆粒規(guī)則排布的簡單模型可能引起PBX有效性能預測的極大偏差; 因此,在構建細觀模型時應當考慮炸藥顆粒空間分布特征。

表2不同分布形式的有限元結果

Table2Finite element results of different particle distributions

modelc/%νE/MPaK/MPaG/MPacircle＿r49．220．49665．29260．711．77circle＿149．000．49645．38249．911．80circle＿249．000．49803．04254．431．01circle＿164．000．490819．50351．866．54circle＿264．000．49667．53370．812．51

Note:νis the Poisson′s ratio;Eis the Young′s modulus;Kis the bulk modulus;Gis the shear modulus;cis content.

圖1c(圓形顆粒circle_2)和圖1d(六邊形顆粒)擁有完全相同的排布形式,顆粒尺寸也大小相當,用于分析顆粒形貌對PBX有效模量的影響。有限元結果見表3,從表3中可以看出不同顆粒形貌間ν和K的偏差均低于2%,E和G對顆粒形貌也更為敏感; 以circle_2為標準,三種顆粒含量下,六邊形顆粒的楊氏模量E分別高出36.18%,114.08%,125.48%,剪切模量G分別高出36.63%,115.14%,126.53%。隨著顆粒含量的增加,顆粒形貌對有效模量的影響在逐漸增大; 顆粒含量很高時,粗略地將顆粒近似為圓形或六邊形也會給有效模量的預測結果帶來較大的偏差; 因此,構建細觀模型時有必要更精準地刻畫PBX的顆粒形貌,而Voronoi模型提供了一種較好的方法。

表3不同顆粒形貌的有限元結果

Table3Finite element results of different particle shapes

modelc/%νE/MPaK/MPaG/MPacircle＿249．000．49803．04254．431．01hexagon49．000．49734．14255．571．38circle＿264．000．49667．53370．812．51hexagon64．000．492816．12371．885．40circle＿272．250．494516．56502．995．54hexagon72．250．487737．34507．4212．55

PBX顆粒含量約90%,顆粒的分布形式和形貌對有效模量的影響非常大。因此與PBX真實細觀形貌更接近的Voronoi模型能更好地預測其有效模量。此外,Voronoi模型能夠實現炸藥顆粒的超高填充度,這是圓形顆粒很難滿足的。表4列舉了顆粒含量為49%～94.9%的Voronoi模型有效模量的有限元結果。隨著顆粒含量的增加,泊松比不斷減小,E、K、G均迅速增大。c=94.9%時有限元結果E=1.410 GPa,而準靜態(tài)下相同含量PBX楊氏模量實驗值為(1.015±0.11) GPa[19]。有限元計算中使用理想的材料參數,炸藥晶體與粘結劑之間的界面處理為共節(jié)點的完美粘接。而實際上組分材料內部存在眾多微裂紋、微孔洞等細觀損傷形式,組分界面間也存在各種缺陷,這都對PBX有效性能有弱化作用,因此有限元結果較實驗值偏大。此外,目前有限元模型中未考慮顆粒粒徑的分布,這可能也是造成差異的原因。

表4不同顆粒含量下Voronoi模型的有限元結果

Table4Finite element results of Voronoi models with different particle contents

c/%νE/MPaK/MPaG/MPa49．000．49675．12258．881．7164．000．492715．51352．425．2072．250．486636．66454．8712．3381．000．478984．24664．1328．4890．250．4341483．431222．28168．5594．900．39401410．842217．03506．06

圖5展示了不同顆粒含量下Voronoi模型計算結果的變化趨勢,同時給出了circle_2模型、hexagon模型的計算結果作為對比。由于組分Estane與HMX的模量相差巨大,顆粒含量對有效模量的影響十分明顯。泊松比ν隨顆粒含量的增加快速下降(圖5a),楊氏模量E(圖5b)、體積模量K(圖5c)、剪切模量G(圖5d)三者隨顆粒含量的增加近似指數上升。各模型之間體積模量差異相對最小。顆粒含量低于75%時,hexagon模型與Voronoi模型計算結果較為一致; 顆粒含量高于75%時,兩者出現明顯差異,hexagon模量更明顯地高估了PBX的有效模量。

a. Poisson′s ratio

b. Young′s modulus

c. bulk modulus

d. shear modulus

圖5各模型有效性能計算結果對比

Fig.5Comparison of the effective property results calculated by different models

4 結 論

(1) 由于PBX顆粒含量極高,組分彈性模量相差3～4個數量級,楊氏模量和剪切模量對炸藥顆粒的分布形式和形貌都很敏感。

(2) 以隨機種子點為中心,采用Voronoi法生成隨機多邊形顆粒,建立Voronoi細觀數值模型,既較真實地反映了PBX的細觀顆粒特征,又實現了炸藥顆粒的高度填充,其楊氏模量預測結果與實測值比較接近。

(3) 基于Voronoi模型的結果表明,PBX的楊氏模量、體積模量、剪切模量均隨著顆粒含量的增加近似指數上升,而泊松比隨顆粒含量的增加快速下降。

參考文獻:

[1] 羅景潤. PBX的損傷、斷裂及本構關系研究. 綿陽: 中國工程物理研究院, 2001.

LUO Jing-run. Study on damage, fracture and constitutive relation of PBX. Mianyang: China Academy of Engineering Physics, 2001.

[2] 黃克智, 黃永剛. 固體本構關系. 北京: 清華大學出版社, 1999: 120-172.

HUANG Ke-zhi, HUANG Yong-gang. Solid constitutive relations. Beijing: Tsinghua University Press, 1999: 120-172.

[3] 胡更開, 鄭泉水, 黃筑平. 復合材料有效彈性性質分析方法. 力學進展, 2001, 31(3): 361-393.

HU Geng-kai, ZHEN Quan-shui, HUANG Zhu-ping. Micromechanics methods for effective elastic properties of composite materials.AdvancesinMechanics, 2001, 31(3): 361-393.

[4] 敬仕明, 李明, 龍新平. PBX有效彈性性能研究進展. 含能材料, 2009, 17(1): 119-124.

JING Shi-ming, LI Ming, LONG Xin-ping. Progress in predicting the effective elastic properties of PBX.ChineseJournalofEnergeticMaterials(HannengCailiao), 2009, 17(1): 119-124.

[5] Rae P J, Goldrein H T, Palmer S J, et al. Quasi-static studies of the deformation and failure ofβ-HMX based polymer bonded explosives.MathematicalPhysical&EngineeringSciences, 2002, 458: 743-762.

[6] Annapragada S R, Sun D W, Garimella S V. Prediction of effective thermo-mechanical properties of particulate composites.ComputationalMaterialsScience, 2007, 40: 255-266.

[7] 賈憲振, 王浩, 王建靈. 炸藥有效彈性性能的細觀尺度仿真預估. 含能材料, 2013, 21(4): 469-472.

JIA Xian-zhen, WANG Hao, WANG Jian-ling. Mesoscale simulation of effective elastic properties of explosive.ChineseJournalofEnergeticMaterials(HannengCailiao), 2013, 21(4): 469-472.

[8] 戴開達, 劉龑龍, 陳鵬萬. PBX炸藥有效彈性模量的有限元模擬. 北京理工大學學報, 2012, 32(11): 1154-1158.

DAI Kai-da, LIU Yan-long, CHEN Peng-wan. Finite element simulation on effective elastic modulus of PBX explosives.TransactionsofBeijingInstituteofTechnology, 2012, 32(11): 1154-1158.

[9] Biswajit B. Effective elastic moduli of PBX from finite element simulations[EB/OL]. arXiv: cond-mat/0510367.

[10] 孫繼忠, 胡艷, 馬永強. 基于Delaunay三角剖分生成Voronoi圖算法. 計算機應用, 2010, 30(1): 75-77.

SUN Ji-zhong, HUYan, MA Yong-qiang. Voronoi diagram generation algorithm based on Delaunay triangulation.JournalofComputerApplication, 2010, 30(1): 75-77.

[11] CHEN Peng-wan, HUANG Feng-lei, DING Yan-sheng. Microstructure, deformation and failure of polymer bonded explosives.JournalofMaterialsScience, 2007, 42: 5272-5280.

[12] Ananda B, Zhou M. A lagrangian framework for analyzing microstructure level response of polymer-bonded explosives.ModelingandSimulationinMaterialsScienceandEngineer, 2011, 19: 1-24.

[13] 張妍, 韓林. 細觀力學基礎. 北京: 科學出版社, 2014: 4-7.

ZHANG Yan, HAN Lin, Foundation of mesomechanics. Beijing: Science Press, 2014: 4-7.

[14] Kanit T, Forest S, Galliet I, et al. Determination of the size of the representative volume element for random composites statistical and numerical approach.InternationalJournalofSolidsandStructure, 2003, 40: 3647-3679.

[15] Sab K. On the homogenization and the simulation of random materials.EuropeanJournalofMechanicsSolids, 1992, 11: 585-607.

[16] Lemaitre J. Formulation and identification of damage kinetic constitutive equations.ContinuumDamageMechanicsTheoryandApplications, 1987, 295: 37-89.

[17] Ostoja S M. Random field models of heterogeneous materials.SolidsandStructures, 1998, 35: 2429-2455.

[18] Nguyen V D, Bechet E, Geuzaine C, et al. Imposing periodic boundary condition on arbitrary meshes by polynomial interpolation.ComputationalMaterialScience, 2012, 55(5): 390-406.

[19] Gray G T, Idar D J, et al. High- and low-strain rate compression properties of several energetic materail composites as a function of strain rate and temperature∥11thinternational detonation symposium, Snowmass, Colorado, 1998.