李潘， 郝志明， 劉永平， 甄文強(qiáng)

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

0 引言

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

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

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

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

1 PD非普通狀態(tài)理論

1.1 本構(gòu)關(guān)系

PD非普通狀態(tài)理論的基本方程[4]為

(1)

(2)

式中：u為物質(zhì)點(diǎn)x的位移；u′為物質(zhì)點(diǎn)x′的位移；η=u′-u為兩物質(zhì)點(diǎn)之間的相對位移。

物質(zhì)點(diǎn)的近似變形梯度為

(3)

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

根據(jù)連續(xù)介質(zhì)力學(xué)，應(yīng)變張量可表示為

(4)

式中：I為單位矩陣。

令連續(xù)介質(zhì)力學(xué)應(yīng)變能和PD應(yīng)變能相等，可得PD力狀態(tài)的表達(dá)式為

(5)

式中：P為第1類Piola-Kirchhoff應(yīng)力；S為第2類Piola-Kirchhoff應(yīng)力。

將(5)式代入(1)式，可得非普通狀態(tài)PD理論的運(yùn)動(dòng)方程為

(6)

1.2 運(yùn)動(dòng)方程求解

為了將PD方法應(yīng)用于準(zhǔn)靜態(tài)問題分析，采用自適應(yīng)動(dòng)態(tài)松弛法對(6)式進(jìn)行求解，運(yùn)動(dòng)方程[20]可以寫成：

(7)

在自適應(yīng)動(dòng)態(tài)松弛法中，為了使運(yùn)動(dòng)方程的解盡快收斂到穩(wěn)定解，每一個(gè)載荷步阻尼系數(shù)cn都需要重新計(jì)算，計(jì)算公式為

(8)

由于物質(zhì)點(diǎn)與其近場范圍內(nèi)其他物質(zhì)點(diǎn)的弱耦合作用，使得結(jié)構(gòu)出現(xiàn)剛體位移，造成PD非普通狀態(tài)理論中出現(xiàn)零能模式?？刂屏隳苣Ｊ叫枰诹顟B(tài)的基礎(chǔ)上添加沙漏力[5]

(9)

該方法是在物質(zhì)點(diǎn)x與其近場范圍內(nèi)的所有點(diǎn)之間增加線性彈簧，其彈性常數(shù)為cⅠ.

(10)

PD將結(jié)構(gòu)離散為均勻分布的物質(zhì)點(diǎn)，通過定義兩點(diǎn)間的對點(diǎn)力構(gòu)建點(diǎn)的運(yùn)動(dòng)方程，求解運(yùn)動(dòng)方程從而得到物質(zhì)點(diǎn)的速度和位移。本文采用Fortran語言開發(fā)了基于PD非普通狀態(tài)理論的損傷破壞模擬程序。

2 損傷定義

PD方法認(rèn)為鍵的伸長率達(dá)到臨界伸長率時(shí)，兩點(diǎn)之間的作用鍵發(fā)生斷裂。鍵的伸長率s表達(dá)式[21]為

(11)

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

(12)

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

下面基于應(yīng)力狀態(tài)定義損傷[13]。首先定義物質(zhì)點(diǎn)x和x′之間的平均第一主應(yīng)力

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

(13)

定義各物質(zhì)點(diǎn)的局部損傷值

(14)

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

(15)

ft為材料拉伸強(qiáng)度。

2.1 損傷對力狀態(tài)的影響

物質(zhì)點(diǎn)之間的力狀態(tài)方程(5)式可以改寫為

(16)

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

(17)

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

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

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

(18)

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

3 巴西圓盤實(shí)驗(yàn)?zāi)M

3.1 計(jì)算模型

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

表1 PBX9501材料參數(shù)[22]

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

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

3.2 結(jié)果與討論

3.2.1 僅考慮損傷對力狀態(tài)的影響模擬結(jié)果

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

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

Tab.2 Damage evolution of Brazilian disc

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

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

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

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

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

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

3.2.3 損傷隨加載位移的變化

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

4 結(jié)論

本文研究了PD非普通狀態(tài)理論求解線彈性損傷破壞問題的方法，并應(yīng)用于PBX巴西圓盤實(shí)驗(yàn)的損傷破壞模擬，獲得了以下結(jié)論：

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

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

參考文獻(xiàn)(References)

[1] 黃丹, 章青, 喬丕忠, 等. 近場動(dòng)力學(xué)方法以及應(yīng)用[J].力學(xué)進(jìn)展, 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òng)力學(xué)理論及其在巖石類材料變形過程的數(shù)值模擬[D].重慶: 重慶大學(xué), 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] 黃西成, 李尚昆, 魏強(qiáng), 等. 基于XFEM與Cohesive模型分析PBX裂紋產(chǎn)生與擴(kuò)展[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] 王竟成, 羅景潤. 不同細(xì)觀力學(xué)方法預(yù)測高聚物粘結(jié)炸藥有效彈性模量的比較[J]. 兵工學(xué)報(bào), 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.