基于近場動力學非普通狀態理論的高聚物粘結炸藥損傷破壞模擬研究

李潘， 郝志明， 劉永平， 甄文強

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

0 引言

近場動力學(PD)方法是一種新的基于非局部思想的無網格方法，主要應用于損傷、破壞、失穩等問題研究[1-2]。最初Silling[3]提出的PD鍵理論采用對點力描述材料的本構關系，使得泊松比必須為定值，且不能準確地反映材料的黏性、塑性。為了克服上述缺點，Silling等[4]又進一步提出以狀態為基礎的PD理論，包括普通狀態理論和非普通狀態理論。其中PD非普通狀態理論引入了應力、應變等基本力學參量，能夠實現傳統本構關系和損傷破壞準則在PD框架下的應用。由于采用節點積分，PD非普通狀態理論引起零能模式問題，需要通過添加沙漏力進行控制[5-6]。與有限元、離散元或擴展有限元等方法相比，PD方法將結構離散為帶質量的物質點，根據空間積分求解運動方程，消除了網格依賴性和奇異性，損傷自然產生，在求解不連續問題時具有明顯優勢。

截止目前，學者們基于PD非普通狀態理論開展了金屬、混凝土、土體等材料的損傷破壞研究。Foster等[7-8]利用PD非普通狀態理論模擬了鋼板在沖擊載荷作用下的裂紋產生和擴展過程。Tupek等[9]將韌性損傷演化模型引入PD，模擬了鋁合金在沖擊載荷作用下的破壞過程。Fan等[10-12]將PD方法與光滑粒子流體動力學(SPH)方法耦合，模擬了土體在爆炸沖擊波作用下的響應。上述工作模擬了結構的損傷演化過程，但沒有對損傷過程中應力的分布情況進行分析。文獻[13-17]研究了混凝土、巖石等脆性材料中裂紋萌生、擴展以及成核過程，并給出了損傷過程中應力分布。但損傷區域的應力并沒有得到應有釋放，反而造成應力隨損傷演化逐漸增大、與實際不符的情況。

高聚物粘結炸藥(PBX)是一種含能材料，廣泛應用于國防和國民經濟各領域。由于力學行為較復雜，目前PBX損傷破壞研究主要采用有限元、離散元以及擴展有限元等方法[18-19]，PD方法在該方面的研究還比較少見。

本文介紹了PD非普通狀態基本理論，求解了運動方程的動態松弛方法，并引入沙漏力控制求解過程中的零能模式。在考慮損傷對PD力狀態影響的基礎上，提出通過影響函數引入損傷的方法，從而考慮了損傷對近似變形梯度和變形張量的影響。該方法應用于PBX巴西圓盤實驗的損傷破壞行為模擬，討論了沙漏力參數的取值。

1 PD非普通狀態理論

1.1 本構關系

PD非普通狀態理論的基本方程[4]為

(1)

(2)

式中：u為物質點x的位移；u′為物質點x′的位移；η=u′-u為兩物質點之間的相對位移。

物質點的近似變形梯度為

(3)

式中：K為形變張量；ω(|ξ|)為作用鍵的影響函數，它是兩點相對位置|ξ|的函數；dVξ為初始構形中物質點x所占的體積。

根據連續介質力學，應變張量可表示為

(4)

式中：I為單位矩陣。

令連續介質力學應變能和PD應變能相等，可得PD力狀態的表達式為

(5)

式中：P為第1類Piola-Kirchhoff應力；S為第2類Piola-Kirchhoff應力。

將(5)式代入(1)式，可得非普通狀態PD理論的運動方程為

(6)

1.2 運動方程求解

為了將PD方法應用于準靜態問題分析，采用自適應動態松弛法對(6)式進行求解，運動方程[20]可以寫成：

(7)

在自適應動態松弛法中，為了使運動方程的解盡快收斂到穩定解，每一個載荷步阻尼系數cn都需要重新計算，計算公式為

(8)

由于物質點與其近場范圍內其他物質點的弱耦合作用，使得結構出現剛體位移，造成PD非普通狀態理論中出現零能模式。控制零能模式需要在力狀態的基礎上添加沙漏力[5]

(9)

該方法是在物質點x與其近場范圍內的所有點之間增加線性彈簧，其彈性常數為cⅠ.

(10)

PD將結構離散為均勻分布的物質點，通過定義兩點間的對點力構建點的運動方程，求解運動方程從而得到物質點的速度和位移。本文采用Fortran語言開發了基于PD非普通狀態理論的損傷破壞模擬程序。

2 損傷定義

PD方法認為鍵的伸長率達到臨界伸長率時，兩點之間的作用鍵發生斷裂。鍵的伸長率s表達式[21]為

(11)

二維各向同性線彈性材料的鍵臨界伸長率sc表達式為

(12)

式中：K為材料的體積模量；G為剪切模量；Gc為斷裂能。

下面基于應力狀態定義損傷[13]。首先定義物質點x和x′之間的平均第一主應力

σ1(x,x′)=[σ1(x)+σ1(x′)]/2.

(13)

定義各物質點的局部損傷值

(14)

式中：φ(x,t)取值范圍在[0,1]之間，0表示該物質點鄰域內無作用鍵斷裂，1表示鄰域內的鍵全部斷裂；μ(x,x′,t)為物質點對是否斷裂的函數，其表達式為

(15)

ft為材料拉伸強度。

2.1 損傷對力狀態的影響

物質點之間的力狀態方程(5)式可以改寫為

(16)

為了表征研究系統的宏觀損傷情況，定義損傷因子d為

(17)

式中：nf為整體離散空間中斷裂鍵總數；ntot為整體離散空間中初始作用鍵總數。

2.2 損傷對近似變形梯度和變形張量的影響

在結構含初始缺陷或在模擬過程中產生損傷時，采用(16)式只是排除了兩物質點之間相互作用力的傳遞，實際上作用鍵的斷裂還將影響到(3)式中近似變形梯度F和變形張量K的計算。目前PD計算中，文獻[13-17]仍考慮已斷裂鍵的作用，造成損傷處應力沒有釋放。為了克服這一困難，通過影響函數引入損傷

(18)

這樣近似變形梯度F和變形張量K中就不再考慮已斷裂鍵的影響，損傷處應力得到相應的釋放。對于應變梯度較大的區域，若物質點近場范圍內的作用鍵完全斷裂，則形變張量K和近似變形梯度F產生奇異，因此有必要對周圍作用鍵完全斷裂的物質點進行特殊處理，即當物質點x與其他物質點不發生作用時，不再對該點的K和F進行計算，直接強制其應變、應力值為0.

3 巴西圓盤實驗模擬

3.1 計算模型

PBX9501的巴西圓盤試件如圖2(a)所示，圓盤半徑R=25.40 mm，厚度W=6.35 mm，將其簡化為平面應力問題進行模擬。將巴西圓盤均勻離散為8 273個粒子，相鄰粒子間的間距Δx=0.25 mm，近場半徑δ=3Δx，離散模型如圖2(b)所示。采用方塊對巴西圓盤施加位移載荷，2 000個加載步，每步施加位移為1.5×10-4mm. PBX9501的材料參數見表1.

表1 PBX9501材料參數[22]

采用1.2節中添加沙漏力的方法進行零能模式控制，通過零能模式對位移場的影響，確定沙漏力參數取值。圖3給出了不同沙漏力參數對x方向位移的影響，從中可以看出：當cI=0時，位移沿x軸的分布曲線產生劇烈振蕩；當cI=0.005×1023時，PD方法計算結果偏小，且位移在x軸的端部產生輕微振蕩；當cI分別為0.1×1023、0.25×1023、0.5×1023時，PD方法模擬結果與有限元方法結果吻合較好；當cI取值繼續增大時，PD方法計算結果越來越偏離有限元方法結果。因此，采用添加沙漏力的方法進行零能模式控制時，參數取值應該控制在一個較精確的取值附近。若取值太小，則零能模式引起的數值不穩定性無法得到有效抑制；若取值太大，則計算結果同樣偏離真實值。

圖4對cI=0和cI=0.5×1023兩種情況下的x軸方向位移云圖進行了對比，結果表明不進行零能模式控制時，位移場產生數值振蕩，分布云圖虛化。選取合適的零能模式控制方法和控制參數，能夠使數值振蕩得到抑制，從而獲得穩定的位移場。

3.2 結果與討論

3.2.1 僅考慮損傷對力狀態的影響模擬結果

采用(16)式引入損傷對力狀態的影響，表2給出了巴西圓盤的損傷演化過程。從表2中可以看出，當加載步達到1 212步時，試件在中心開始起裂，形成豎向裂紋，此時中心處表征物質點局部損傷的參數φ約為0.2，即物質點周圍的作用鍵發生部分斷裂。隨著加載位移的進一步增大，裂紋上下尖端產生應力集中，使得裂紋尖端應力大于材料拉伸強度，裂紋向上下兩端快速擴展。在裂紋向兩端擴展的同時，已開裂處物質點的損傷進一步加劇，當加載步達到12 16步時，中心處物質點發生斷裂。隨著越來越多物質點的斷裂，宏觀裂紋形成。當加載步達到1 230步時，裂紋沿著載荷施加方向形成一條長長的裂紋。從x軸方向位移云圖可以看出，損傷的初始階段，位移基本上是連續的，隨著裂紋的擴展以及宏觀裂紋的形成，位移表現出明顯的不連續性。

表2 巴西圓盤的損傷演化過程

Tab.2 Damage evolution of Brazilian disc

將局部損傷值φ(x,t)≥0.99的物質點刪除，從而給出PD方法模擬得到的巴西圓盤試件破壞模式，并與實驗結果進行對比。如圖5所示，PD數值計算在最后形成一條豎向裂紋，損傷破壞過程與實驗基本吻合。但在損傷過程中由于應力得不到有效釋放，導致損傷處應力越來越大，影響了損傷過程中的位移和應力分布模擬。造成PD方法模擬結果刪除了中間三排物質點，形成較寬的裂紋。

3.2.2 考慮損傷對近似變形梯度和變形張量影響的模擬結果

從表2中的應力云圖可以看出，損傷處應力沒有得到釋放，且數值隨損傷演化增大，與實際情況不符。為了解決上述問題，在3.2.1節考慮損傷對PD力狀態影響的基礎上，再引入損傷對近似變形梯度和變形張量的影響，模擬巴西圓盤的損傷破壞過程。表3給出了改進損傷引入方法之后巴西圓盤的損傷演化過程以及應力分布云圖。從表3中可以看出，損傷處應力得到釋放，隨損傷演化在裂尖處重新集中，從而促使裂紋繼續向前擴展，最后形成宏觀裂紋，位移呈現明顯不連續性。

將局部損傷值φ(x,t)≥0.99的物質點刪除得到巴西圓盤的破壞模式，如圖6所示。結合圖5對比分析發現，除少數位置外巴西圓盤只刪除了中間一排點，比之前刪除三排點能夠更準確地表征真實裂紋開裂情況。但是由于PD方法是一種非局部方法，且非普通狀態理論中物質點之間的關聯性更強，使得損傷演化過程中有些作用鍵無法及時斷裂，造成損傷跨點傳播。

表3 通過影響函數引入損傷后巴西圓盤的損傷演化過程

Tab.3 Damage evolution of Brazilian disc after introducing damage through the influence function

3.2.3 損傷隨加載位移的變化

圖7將PD模擬結果和實驗結果載荷- 位移曲線進行比較，二者吻合較好。PD方法計算得到的試件強度比實驗結果稍高，這是由于PD方法采用線彈性損傷模型，沒有反映材料的塑性損傷特征。圖8給出了載荷和損傷因子隨位移的變化，當位移達到0.134 5 mm時，巴西圓盤開始產生裂紋，載荷- 位移曲線呈現非線性特征；隨后裂紋迅速擴展，導致載荷劇烈下降；裂紋繼續擴展直至貫穿，試件發生完全斷裂，與Zhou等[24]測得實驗現象一致性很好。

4 結論

本文研究了PD非普通狀態理論求解線彈性損傷破壞問題的方法，并應用于PBX巴西圓盤實驗的損傷破壞模擬，獲得了以下結論：

1) PD非普通狀態理論在不進行零能模式控制的情況下，數值振蕩明顯，得不到合理結果。采用添加沙漏力的方式進行零能模式控制時：若沙漏力參數太小，則零能模式無法得到有效抑制；若沙漏力太大，則結果仍將偏離真實值。因此，需將沙漏力控制在一個合理范圍內。

2) 提出通過影響函數的損傷引入方法，從而考慮損傷對近似變形梯度和變形張量的影響，解決了發生損傷處的應力釋放問題。模擬了PBX巴西圓盤實驗中的裂紋萌生、擴展直至貫穿破壞全過程，試件載荷- 位移曲線、破壞模式與實驗結果一致。

參考文獻(References)

[1] 黃丹, 章青, 喬丕忠, 等. 近場動力學方法以及應用[J].力學進展, 2010, 40(4): 448-459.

HUANG Dan, ZHANG Qing, QIAO Pi-zhong, et al. A review on peridynamics method and its applications[J]. Advances in Mechanics, 2010, 40(4):448-459. (in Chinese)

[2] Gerstle W, Sau N, Silling S A. Peridynamic modeling of concrete structures[J]. Nuclear Engineering and Design, 2007, 237(12/13):1250-1258.

[3] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1):175-209.

[4] Silling S A, Epton M, Weckner O, et al. Peridynamic states and constitutive modeling[J]. Journal of Elasticity, 2007, 88(2):151-184.

[5] Breitenfeld M S, Geubelle P H, Weckner O, et al. Non-ordinary state-based peridynamic analysis of stationary crack problems[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 272(11):233-250.

[6] Silling S A. Stability of peridynamic correspondence material models and their particle discretizations[J]. Computer Methods in Applied Mechanics and Engineering,2017, 322(15):42-57.

[7] Foster J T, Silling S A, Chen W W. Viscoplasticity using peridynamics[J]. International Journal for Numerical Methods in Engineering, 2010, 81(10):1242-1258.

[8] Foster J T. Dynamic crack initiation toughness: experiments and peridynamic modeling[D]. Lafayette, IN, US: Puedue University, 2009.

[9] Tupek M R, Rimoli J J, Radovitzky R. An approach for incorporating classical continuum damage models in state-based peridynamics[J]. Computer Methods in Applied Mechanics and Engineering,2013, 263(24):20-26.

[10] Fan H, Bergel G L, Li S. A hybrid peridynamics-SPH simulation of soil fragmentation by blast loads of buried explosive[J]. International Journal of Impact Engineering, 2016, 87(1):14-27.

[11] Lai X, Ren B, Fan H F, et al. Peridynamics simulations of geomaterial fragmentation by impulse loads[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(12):1304-1330.

[12] Ren B, Fan H F, Bergel G L, et al. A peridynamics-SPH coupling approach to simulate soil fragmentation induced by shock waves[J]. Computational Mechanics, 2015, 55(2):287-302.

[13] Zhou X P, Wang Y T, Xu X M. Numerical simulation of initiation, propagation and coalescence of cracks using the non-ordinary state-based peridynamics[J]. International Journal of Fracture, 2016, 201(2):213-234.

[14] Zhou X P, Wang Y T, Qian Q H. Numerical simulation of crack curving and branching in brittle materials under dynamic loads using the extended non-ordinary state-based peridynamics[J]. European Journal of Mechanics A/Solids, 2016, 60(5):277-299.

[15] Zhou X P, Wang Y T. Numerical simulation of crack propagation and coalescence in pre-cracked rock-like Brazilian disks using the non-ordinary state-based peridynamics[J]. International Journal of Rock Mechanics and Mining Sciences, 2016, 89(5):235-249.

[16] Wang Y T, Zhou X P, Xu X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics[J]. Engineering Fracture Mechanics, 2016, 163(20):248-273.

[17] 谷新保. 近場動力學理論及其在巖石類材料變形過程的數值模擬[D].重慶: 重慶大學, 2015.

GU Xin-bao. Peridynamic theory and its numerical simulation in the deformation and damage process of the rock-like materials[D]. Chongqing: Chongqing University, 2015. (in Chinese)

[18] 黃西成, 李尚昆, 魏強, 等. 基于XFEM與Cohesive模型分析PBX裂紋產生與擴展[J]. 含能材料, 2017, 25(8):694-700.

HUANG Xi-cheng, LI Shang-kun, WEI Qiang, et al.Analysis of crack initiation and growth in PBX energetic material using XFEM-based cohesive method[J]. Chinese Journal of Energetic Materials, 2017, 25(8):694-700. (in Chinese)

[19] 王竟成, 羅景潤. 不同細觀力學方法預測高聚物粘結炸藥有效彈性模量的比較[J]. 兵工學報, 2017, 38(6):1106-1112.

WANG Jing-cheng, LUO Jing-run. Comparison of effective moduli of polymer bonded explosive predicted by different micromechanical methods[J]. Acta Armamentarii, 2017, 38(6):1106-1112. (in Chinese)

[20] Kilic B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials[D]. Ann Arbor, MI, US: The University of Arizona, 2008.

[21] Kilic B, Agwai A, Madenci E. Peridynamic theory for progressive damage prediction in center-cracked composite laminates[J]. Composite Structures, 2009, 90(2): 141-151.

[22] Tan H, Liu C, Huang Y, et al. The cohesive law for the particle/matrix interfaces in high explosives[J]. Journal of the Mechanics and Physics of Solids, 2005, 53(8):1892-1917.

[23] Liu C, Brownign R. Fracture in PBX 9501 at low rates[C]∥Proceedings of the 12th International Detonation Symposium. San Diego, CA, US: Office of Naval Research, 2002: 11-16.

[24] Zhou Z B, Chen P W, Guo B Q, et al. Quasi-static tensile deformation and fracture behavior of a highly particle-filled compo-site using digital image correlation method[J]. Theoretical and Applied Mechanics Letters, 2011(5):10-13.

[25] Liu C, Thompson D G. Macroscopic crack formation and extension in pristine and artificially aged PBX 9501[R]. Santa Fe, NM, US: Los Alamos National Laboratory, 2010.