基于Voronoi細觀數(shù)值模型預(yù)測PBX的有效彈性模量

王竟成, 羅景潤

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

1 引 言

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

(1)PBX的顆粒含量很高,而傳統(tǒng)細觀力學(xué)方法預(yù)測顆粒含量較高的復(fù)合材料有效性能會出現(xiàn)很大偏差;

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

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

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

2 數(shù)值模擬

2.1 細觀模型的建立

為考察顆粒形貌和排列分布方式對PBX有效模量的影響,首先建立四種構(gòu)型的二維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的細觀形貌,實現(xiàn)炸藥顆粒的高度填充,更有效地模擬PBX的有效力學(xué)性能,建立不規(guī)則多邊形顆粒的Voronoi模型。在一定區(qū)域內(nèi)播撒隨機種子點,采用Voronoi法[10]以種子點為中心生成隨機多邊形顆粒。顆粒的平均尺寸通過種子數(shù)進行控制。為保證顆粒生成的形態(tài)質(zhì)量,有必要篩選種子點。圖2b為具有周期性結(jié)構(gòu)的Voronoi圖,炸藥顆粒的平均粒徑為100 μm。進一步生成粘結(jié)劑層后,將節(jié)點信息導(dǎo)入ANSYS獲得有限元模型(圖2c),它與PBX真實的細觀形貌(圖2a)具有較好的相似性。將種子點固定,改變粘結(jié)劑的厚度構(gòu)建不同顆粒含量的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)晶體,粘結(jié)劑選用聚酯型聚氨酯(Estane),均采用各向同性的彈性本構(gòu)模型,其理想材料參數(shù)見表1[12]。從表1中可見基體Estane和炸藥顆粒HMX之間的彈性性能有巨大的差異,泊松比相差一倍,楊氏模量相差四個數(shù)量級。

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

適當(dāng)尺寸的RVE,其有效彈性性質(zhì)不受邊界條件類型的影響[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,高分子和顆粒復(fù)合材料: 1 mm×1 mm×1 mm。Ostoja-Starzewski[17]研究片狀或針狀顆粒復(fù)合材料的有效剛度后發(fā)現(xiàn),在最大許可誤差為5%和顆粒剛度是基體剛度100倍的情況下,RVE最小尺寸是顆粒特征尺寸的10～20倍。真實炸藥顆粒的平均尺寸約100 μm,因此選擇PBX二維RVE模型的尺寸為1000 μm×1000 μm。

2.3 有限元計算方法

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

ui-uj=0

(1)

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

圖3周期性的Voronoi模型

Fig.3Periodic Voronoi models

圖4邊界條件

Fig.4Boundary conditions

根據(jù)平面應(yīng)變假設(shè)和邊界條件,由:

(2)

可得:

(3)

(4)

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

根據(jù)文獻[18],平均應(yīng)力既可以使用體積分得到,也可以通過面積分得到:

(5)

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

(6)

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

(7)

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

(8)

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

3 結(jié)果與討論

表2列舉了圓形顆粒不同分布形式之間的有限元結(jié)果,用于分析顆粒分布形式對有效模量的影響。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有效性能預(yù)測的極大偏差; 因此,在構(gòu)建細觀模型時應(yīng)當(dāng)考慮炸藥顆?？臻g分布特征。

表2不同分布形式的有限元結(jié)果

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

表3不同顆粒形貌的有限元結(jié)果

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

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

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

a. Poisson′s ratio

b. Young′s modulus

c. bulk modulus

d. shear modulus

圖5各模型有效性能計算結(jié)果對比

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

4 結(jié) 論

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

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

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

參考文獻:

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

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

[2] 黃克智, 黃永剛. 固體本構(gòu)關(guān)系. 北京: 清華大學(xué)出版社, 1999: 120-172.

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

[3] 胡更開, 鄭泉水, 黃筑平. 復(fù)合材料有效彈性性質(zhì)分析方法. 力學(xué)進展, 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] 賈憲振, 王浩, 王建靈. 炸藥有效彈性性能的細觀尺度仿真預(yù)估. 含能材料, 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炸藥有效彈性模量的有限元模擬. 北京理工大學(xué)學(xué)報, 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圖算法. 計算機應(yīng)用, 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] 張妍, 韓林. 細觀力學(xué)基礎(chǔ). 北京: 科學(xué)出版社, 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.