基于EFG方法的PBX炸藥圓弧巴西實驗的數值模擬

崔云霄，陳鵬萬，戴開達，王雷元，劉龑龍

(1. 北京理工大學爆炸科學與技術國家重點實驗室，北京100081； 2. 西北核技術研究所，陜西西安710024)

?

基于EFG方法的PBX炸藥圓弧巴西實驗的數值模擬

崔云霄1,2，陳鵬萬1，戴開達1，王雷元2，劉龑龍1

(1. 北京理工大學爆炸科學與技術國家重點實驗室，北京100081； 2. 西北核技術研究所，陜西西安710024)

摘要：為了對準靜態圓弧巴西實驗過程中PBX試樣的裂紋形成和擴展過程進行模擬，基于EFG方法和內聚力模型，利用LS-DYNA軟件的隱式分析模擬計算了實驗過程。結果表明，計算結果和文獻報道的實驗結果吻合較好，證明了該方法的可行性；采用線性內聚力模型可較好地描述PBX炸藥的斷裂過程，PBX試樣的起裂發生在壓縮載荷達到峰值后，隨著壓頭的進一步壓縮，試樣中心的裂紋沿直徑向壓頭方向擴展。

關鍵詞：固體力學；圓弧巴西實驗；EFG方法；內聚力模型；PBX；LS-DYNA軟件

引言

高聚物黏結炸藥(PBX)是由HMX炸藥、聚合物和添加劑組成的一種固體高能炸藥。研究其損傷斷裂性能，對于熱點點火預測、裝藥安全性評估具有重要意義[1]。在外部載荷作用下，PBX炸藥內部可能產生炸藥顆粒脫粘、微裂紋演化等宏觀或細觀損傷，使力學性能劣化，甚至會影響“熱點”的形成。

目前，PBX的損傷斷裂性能研究主要集中于拉伸斷裂，如直接拉伸實驗、巴西實驗、三點彎曲實驗等。相對而言，巴西實驗的試樣制備簡單、所需材料少、方便快捷[2]。傳統巴西實驗在實驗過程中加載接觸點容易產生應力集中，使得該位置可能提前發生破壞，影響測試結果的有效性[3]。Mellor等[4-7]采用圓弧巴西實驗，通過把壓頭接觸面改進為弧形，增大了接觸面，提高了實驗結果的可重復性。由于PBX炸藥兼有黏彈性和準脆性的力學特征，準確描述其力學響應及斷裂特性比較困難。黃濤[8]采用流形元法對PBX炸藥的損傷破壞過程進行了數值模擬，將建立在實驗測量上的彈性損傷模型引入到流形元程序，通過有預制裂紋和無預制裂紋兩種條件下的巴西實驗，探討了損傷對裂紋擴展過程的影響。趙四海[9]利用Visco-SCRAM模型計算了圓弧巴西實驗中PBX炸藥圓盤試樣的力學響應，結果顯示，在圓盤與圓弧墊塊接觸位置損傷變量較大，損傷呈啞鈴分布。

無網格迦遼金方法(EFG)是由Belytschko T提出的一種無網格數值計算方法[10]，具有較好的協調性和穩定性，在動力有限元軟件LS-DYNA中可進行基于EFG方法的模擬分析。本研究基于結合內聚力模型的EFG斷裂模擬方法，計算了準靜態加載圓弧巴西實驗中PBX炸藥的力學響應和斷裂行為，以期為含能材料損傷破壞行為的研究提供參考。

1裂紋模擬方法

1.1EFG方法

EFG方法的基本思想是：計算域Ω內的函數u(X)，可根據該域內的若干已知點，由最小二乘法構造其近似函數，即

(1)

式中：p(X)為基函數；a(X)為與X相關的系數，由二次方程式J(X)求極小值得到

(2)

式中：wi(X)為非負的權函數，其有一個影響范圍，超出影響范圍時，權函數為0。

新版LS-DYNA軟件中引入了結合內聚力模型的EFG斷裂擴展算法[11]，借助可視判據，程序自動將內聚力模型定義在單元邊緣。起裂判據采用最大拉應力判據，裂紋張開過程采用內聚力模型描述。使用該方法模擬計算材料或結構的斷裂問題，不需要在材料或結構中預置裂紋，對于三維模型，只要建立四面體單元作為背景網格。

1.2內聚力模型

內聚力模型本質上是一種唯象模型，其張力位移關系是對裂紋前端斷裂過程區的一種連續介質損傷力學描述。當達到最大應力時，內聚力區域開始起裂，而當達到臨界張開位移時，單元刪除，物理裂紋出現，此時，內聚力作功等于臨界能量釋放速率。基于EFG方法模擬材料或結構的斷裂破壞，內聚力模型采用的是Tvergaard-Hutchinson模型[12]，該模型是三線性(梯形)張力位移關系，如圖1所示。

圖1　Tvergaard-Hutchinson模型的梯形張力位移關系Fig.1　Trapezoid relation between the traction anddisplacement of Tvergaard-Hutchinson model

Tvergaard-Hutchinson模型使用一個無量綱張開位移λ來度量法向斷裂(λ3)和切向斷裂(λ1，λ2)相對位移之間的相互作用

(3)

式中：NLS和TLS分別為界面法向和切向的最大位移張開值；<>為McCauley括號，用來區分拉伸(λ3≥0)和壓縮(λ3<0)。

張力的控制方程為：

(4)

(5)

在數值模擬中，除了最大應力值及臨界張開位移作為參數外，還要選擇給出δ1和δ2。研究表明[13]，對于脆性材料，線性內聚力模型比較合適，因此將δ1和δ2設為0，Tvergaard-Hutchinson模型的張力位移關系退化為線性關系。

2圓弧巴西實驗及計算模型

實驗材料為PBX9501，其配方(質量分數)為：HMX 95%、Estane 2.5%、NP 2.5%，圓盤尺寸為Φ25.4mm×6.35mm，常曲率墊塊與圓盤的半徑比為1.25，MTS材料試驗機常速率加載，頂部圓弧壓頭的移動速率為0.5mm/min，實驗裝置如圖2所示[6]。根據實驗裝置，建立三維全尺寸計算模型。

圖2　實驗裝置示意圖Fig.2　Schematic diagram of the experimental apparatus

采用四面體單元對巴西圓盤試樣進行網格劃分，單元總數為97344個，節點總數為18729個，劃分后的模型如圖3所示。采用六面體單元對兩個圓弧壓頭進行網格劃分，單元總數為9680個，節點總數為12222個；上壓頭施加位移邊界，模擬壓頭運動，下壓頭施加固定邊界，壓頭與試樣之間采用自動面面接觸。

圖3　巴西圓盤試樣的網格劃分模型Fig.3　FEM mesh model of Brazilian disk specimen

為了模擬實驗采用的常速率加載，顯式分析需要較長的時間，因而采用隱式分析。壓頭為鋼材料，采用彈性本構，密度為7850kg/m3，彈性模量為211GPa，泊松比為0.3。對于PBX炸藥，斷裂前將其力學行為近似描述為彈性，材料參數為：密度1.829g/cm3，彈性模量2.8GPa，泊松比0.3。達到起裂判據后，結合內聚力模型描述其內部裂紋的形成。參考PBX9501的抗拉強度，內聚拉伸強度取1.8MPa。選定內聚拉伸強度后，根據斷裂能采用線性內聚力模型估算臨界張開位移。根據文獻[14]，其斷裂韌性約0.29MPa·m1/2，平面應變下斷裂能與斷裂韌性的關系為

(6)

式中：KIC為斷裂韌性；v為泊松比；E為彈性模量。

因此，由式(6)換算得到斷裂能，并根據式(5)可以確定出臨界張開位移。在此基礎上經過試算，采用的PBX材料的臨界張開位移為0.032mm。選擇過大的臨界張開位移，試樣在2%的直徑變化下都不發生斷裂；而選擇過小的臨界張開位移，試樣會提前破壞，且網格發生畸變。

3計算結果與討論

模擬得到壓頭施加的壓縮載荷隨試樣直徑變化的曲線見圖4，載荷(p)利用試樣的厚度(w)和半徑(R)歸一化。圖4中還給出了試樣中心的裂紋張開位移曲線，該曲線是根據試樣中心位置單元兩側節點的x方向位移曲線相減得到的，裂紋張開位移(δ)以試樣直徑(D)歸一化。

圖4中A、B、C和D四個圓圈分別對應40.65、44.40、45.25和46.45s四個時刻，其中A點為壓縮載荷達到峰值前的一個時刻，試樣尚處于彈性變形，直徑變化約1.3%，內部微裂紋開始演化，裂紋張開位移曲線的斜率已發生變化；B點為壓縮載荷達到峰值11.06MPa的時刻；C點為壓縮載荷開始下降而整體變形仍在增加的時刻，材料發生軟化。由圖4可見，裂紋張開位移曲線的斜率在B點和C點之間發生突變，如圖4中的實心原點，對應無量綱位移為1.48%，即45.1s時，試樣起裂。實際實驗中，采用1N的載荷預壓，計算中未考慮預壓過程，而在無量綱位移達到0.1%前(約3s)，模擬得到壓頭載荷無較大變化，扣除該部分時間，則起裂時間為42.1s，與基于DIC方法測試得到的起裂時間37.4s較為接近[6]。D點為壓縮載荷下降到局部最低值的時刻，隨著壓頭位移的增加，載荷將開始回升。起裂后，裂紋張開位移快速增加，并開始沿直徑向壓頭方向擴展。

圖4　施加載荷隨PBX試樣直徑變化的曲線Fig.4　Curves of applied load vs. variation of diameterof PBX sample from simulation

試樣在不同時刻沿x方向的位移云圖見圖5，分別對應圖4的A、B、C、D時刻。

圖5　不同時刻PBX試樣x方向的位移云圖Fig.5　Selected contour plots of x-displacement of PBXsample from simulation at different time

從圖5可以看出，C時刻在試樣的中心附近出現宏觀裂紋，到D時刻后，裂紋已向兩端擴展，且裂紋寬度增加。采用DIC方法得到試樣在C時刻和D時刻沿x方向位移云圖見圖6[6]。

圖6　DIC方法測得PBX試樣x方向的位移云圖Fig.6　Selected contour plots of x-displacement ofPBX sample from DIC measurement

從圖6可以看出，數值模擬得到的位移云圖與采用DIC方法獲取的位移云圖具有相似的分布，試樣沿直徑方向存在一條宏觀裂紋。在C時刻，模擬得到x方向的位移峰值比DIC方法得到的結果稍偏大。

試樣在不同時刻x方向的應變云圖見圖7，分別對應圖4的A、B、C、D時刻。

圖7　不同時刻PBX試樣x方向的應變云圖Fig.8　Selected contour plots of x-strain of PBXsample from simulation at different time

由圖7可以看出，在圓盤的中心周圍有一個應變等值線圈，位于等值線圈內的單元處于拉伸狀態。B時刻試樣中心出現局部高應變區；C時刻試樣中心出現宏觀裂紋，裂紋周圍還存在高應變區；D時刻裂紋向兩端擴展，試樣內部的應變能釋放，裂紋附近的拉應變減小，裂紋尖端附近仍存在應變集中區。

采用DIC方法計算得到起裂時刻和裂紋擴展過程x方向的應變云圖見圖8[6]。試樣在不同時刻y方向的位移云圖見圖9。

對比圖7和圖8可知，模擬得到的x方向的應變結果和實驗結果較一致，應變峰值約為6%，說明基于EFG的斷裂模擬方法對PBX材料斷裂破壞過程進行模擬分析是可行的。

圖8　DIC方法測得PBX試樣x方向的應變云圖Fig.8　Selected contour plots of x-strain of PBXsample from DIC measurement

圖9　不同時刻試樣y方向的位移云圖Fig.9　Selected contour plots of y-strain of samplefrom simulation at different time

由圖9可以看出，試樣中與壓頭接觸的部位位移峰值約0.3mm，沒有出現過早破壞，說明采用圓弧壓頭后改善了標準巴西實驗中試樣端頭的應力集中問題。

從數值模擬的結果看，在施加的壓縮載荷達到峰值時，PBX炸藥并未發生宏觀斷裂，而是在緊跟載荷峰值之后，圓盤試樣的中心出現裂紋，與實驗結果一致。隨著壓頭位移的增加，裂紋沿試樣直徑向兩端擴展。需要指出的是，目前該計算還存在不穩定的情況，其方法還在發展之中。

4結論

(1)利用LS-DYNA程序的EFG方法，結合內聚力模型，對PBX圓弧巴西實驗進行了數值模擬，再現了準靜態壓縮下PBX圓盤的裂紋萌生及擴展的過程，計算得到的變形場分布與實驗結果符合較好。

(2)采用線性演化內聚力模型及合適的參數可以較好地描述PBX炸藥斷裂破壞過程中的能量耗散行為，根據模擬計算結果可知，試樣起裂發生于壓縮載荷達到峰值之后。

參考文獻：

[1]陳鵬萬, 黃風雷. 含能材料損傷理論及應用[M]. 北京: 北京理工大學出版社, 2006.

CHEN Peng-wan, HUANG Feng-lei. Damage Theory and Application of Energetic Materials[M]. Beijing: Beijing Institute of Technology Press, 2006.

[2]Johnson H D. Mechanical properties of high explosives, MHSMP-74-19[R]. Amarillo: Mason and Hanger-silas Mason Co., Inc., 1974．

[3]龐海燕, 李明, 溫茂萍, 等. PBX巴西試驗與直接拉伸試驗的比較[J]. 火炸藥學報, 2011, 34(1):42-44.

PANG Hai-yan, LI Ming, WEN Mao-ping, et al. Comparison on the Brazilian test and tension test of the PBX[J]. Chinese Journal of Explosives and Propellants, 2011, 34(1):42-44.

[4]Mellor M, Hawkes I. Measurement of tensile strength by diametric compression of discs and annuli[J]. Eng Geol, 1971, 5:173-225.

[5]陳鵬萬, 黃風雷, 張瑜. 用巴西實驗評價炸藥的力學性能[J].兵工學報,2001,22(4): 533-537.

CHEN Peng-wan, HUANG Feng-lei, ZHANG Yu. Brazilian test and its application in the study of the mechanical properties of explosives[J]. Acta Armamentarii, 2001,22(4): 533-537.

[6]Liu C, Thompson D G, Lovatoy M L, et al. Macroscopic crack formation and extension in pristine and artificially aged PBX 9501[C]∥Proceedings of the 14th International Detonation Symposium. Idaho: Office of Naval Research, 2010.

[7]陳科全, 藍林鋼, 路中華, 等. 含預制缺陷PBX炸藥的力學性能及破壞形式[J]. 火炸藥學報, 2015, 38(5):51-53.

CHEN Ke-quan, LAN Lin-gang, LU Zhong-hua, et al. Mechanical properties and failure modes of PBX explosives with different prefabricated defects[J]. Chinese Journal of Explosives and Propellants, 2015, 38(5):51-53.

[8]黃濤. 高聚物粘結炸藥損傷破壞的流形元法模擬研究[D]. 北京：北京理工大學，2006.

HUANG Tao. Numerical simulation of damage and fracture of polymer bonded explosives by manifold method[D]. Beijing: Beijing Institute of Technology, 2006.

[9]趙四海. 用粘彈性統計裂紋模型模擬高能炸藥的力學響應和非沖擊點火[D]. 長沙: 國防科技大學, 2011.

ZHAO Si-hai. A viscoelastic statistic crack constitutive model for mechanical response and the Non-shock ignition of high explosives[D]. Changsha: National University of Defense Technology, 2011.

[10]Belytschko T, Lu Y Y, Gu L. Crack propagation by element-free galerkin methods[J]. Engineering Fracture Mechanics, 1995, 51: 295-315.

[11]GUO Yong, Wu C T. XFEM and EFG cohesive fracture analysis for brittle and semi-brittle materials[C]∥11th International LS-DYNA Users Conference. Livermore: Livermore Software Technology Corp, 2010, 421-432.

[12]Tvergaard V, Hutchinson J W. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids[J]. Journal of the Mechanics and Physics of Solids, 1992, 40(6): 1377-1397.

[13]Xu X P, Needleman A. Analysis of ductile crack growth by means of a cohesive damage model[J]. International Journal of Fracture, 1996, 81(2): 99-112.

[14]LIU Cheng, Cady C M, Rae P J, et al. On the quantitative measurement of fracture toughness in high explosive and mock materials[C]∥Proceedings of the 14th International Detonation Symposium. Idaho: Office of Naval Research, 2010.

Numerical Simulation of PBX under Circular Anvil Brazilian Test Based on EFG Method

CUI Yun-xiao1,2，CHEN Peng-wan1，DAI Kai-da1，WANG Lei-yuan2，LIU Yuan-long1

(1. State Key Laboratory of Explosion Science and Technology，Beijing Institute of Technology, Beijing 100081， China;2.Northwest Institute of Nuclear Technology，Xi′an 710024，China)

Abstract:To simulate the crack formation and extension process of PBX sample in the quasi-static circular anvil Brazilian test, based on EFG method and cohesive model, the experimental process was calculated by implicit analysis and simulation of LS-DYNA software. The results show that the calculated results are in good agreement with reported experimental ones in literature, proving the feasibility of the method. The fracture process of PBX can be well described by linear cohesive model. The crack initiation of PBX sample occurs after the compression load reaches the peak value, and with the further compression of compression anvil, the crack of the specimen center extends along the diameter direction to the anvil.

Keywords:solid mechanics；circular anvil Brazilian test；EFG method；cohesive model；PBX；LS-DYNA software

中圖分類號：TJ55; O346

文獻標志碼：A

文章編號：1007-7812(2016)01-0034-05

作者簡介:崔云霄(1980-)，男，博士研究生，從事裝藥安全性和點火機理研究。E-mail:yunxiaocui@163.com

基金項目:國家自然科學基金資助(No.11202027，No.11221202，No.U1330202)

收稿日期:2015-08-29；修回日期:2016-01-04

DOI：10.14077/j.issn.1007-7812.2016.01.005