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

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

1 引 言

近場動力學(PD)是一種新的基于非局部思想的無網(wǎng)格方法。在分析破壞問題時有限元、有限差分等傳統(tǒng)方法會產生裂紋尖端奇異性,擴展有限元雖然已經(jīng)解決了很多裂紋擴展和連接問題,但在分析三維裂紋擴展和群裂紋等復雜破壞問題時面臨挑戰(zhàn)。無網(wǎng)格法消除了網(wǎng)格依賴性,但在模擬斷裂時遇到了張力不穩(wěn)定問題。分子動力學方法也被用來模擬裂紋擴展和連接問題,但存在計算時間長、計算效率低等問題。PD方法兼有分子動力學方法和無網(wǎng)格方法的優(yōu)點,避免了傳統(tǒng)方法在面臨不連續(xù)問題時的奇異性,又突破了分子動力學方法在計算尺度上的缺陷。因此,該方法在研究損傷、斷裂、失穩(wěn)等問題時具有明顯優(yōu)勢。PD方法還有待發(fā)展,傳統(tǒng)的粘彈性、塑性以及彈、塑、粘性耦合的材料性質在PD本構模型中的表述尚待深入研究[1-3]。

瀝青、混凝土、高聚物、固體推進劑、高聚物粘結炸藥等廣泛應用于國民經(jīng)濟建設和國防工業(yè)中,在溫度和機械載荷的作用下呈現(xiàn)出明顯的非線性粘彈性特征,這類材料力學行為的研究越來越受到重視[4-6]。孟紅磊[7]提出了一種含累積損傷的非線性粘彈性本構方程來描述推進劑的拉伸應力-應變關系,并將其引入有限元分析中獲得了較好的計算結果。馮震宇[8]將非線性粘彈性朱-王-唐本構模型應用于飛機風擋的數(shù)值模擬,結果表明本構模型較好地模擬了風擋材料的力學行為。

由于傳統(tǒng)數(shù)值模擬方法易在損傷處產生奇異,且具有網(wǎng)格依賴性等缺點,PD方法已開始應用于粘彈性材料的力學行為模擬。Mitchell[9]基于并聯(lián)形式的Maxwell模型建立了PD線性粘彈性本構模型。Azizi[10]將Burgers模型引入PD鍵理論,得到PD線性蠕變本構力函數(shù),模擬了高分子材料在低應力水平下的線性蠕變,數(shù)值模擬與實驗結果較吻合,但模型本構力函數(shù)中的橫截面積對于點對相互作用沒有明確的物理意義。由于高分子、混凝土、高聚物粘結炸藥等材料在高應力作用下呈現(xiàn)出非線性蠕變特性,模型不能較好地模擬這類材料在溫度和應力共同作用下的非線性蠕變行為。

本研究利用Burgers粘彈性模型表征蠕變柔度主曲線,結合非線性粘彈體的時間-應力等效原理,得到不同應力水平下蠕變柔度的表達式。將其代入PD本構力函數(shù),推導出非線性粘彈體的PD蠕變本構力函數(shù),從而建立起一種可應用于高聚物粘結炸藥的近場動力學蠕變模擬方法。利用該方法對PBX9502在不同溫度和不同應力作用下的蠕變行為進行模擬,獲得與實驗一致的結果。

2 PD方法的基本理論

2.1 運動方程

PD方法將物體所占區(qū)域離散成具有一定質量的物質點,域內任一物質點xk與其周圍一定范圍(近場范圍H=H(xk,δ)={xj∈R:‖xj-xk‖≤δ})內的其它物質點xj之間存在相互作用f,如圖1所示,也稱為本構力函數(shù)[11]：

f=f(xk,xj,u(xk,t),u(xj,t),t)

(1)

根據(jù)牛頓第二定律可得物質點xk的運動方程為

b(xk,t)

(2)

式中,ρ為物質點密度,u為物質點的位移,b為單位體積物質所受的外載荷。

圖1物質點與其近場范圍內其它物質點的相互作用

Fig.1Pairwise interaction of a material point with its neighboring points

2.2 本構力函數(shù)的基本形式

對于線彈性各向同性材料,Madenci[12]給出了本構力函數(shù)f的基本形式

(3)

式中,yj-yk=xj-xk+uj-uk;s=(|yj-yk|- |xj-xk|)/(|xj-xk|),為作用鍵伸長率;c為與結構尺寸、近場半徑和材料柔度相關的參數(shù)。令經(jīng)典理論的應變能密度與PD理論的應變能密度相等,可分別得到一維、二維和三維情況下c的表達式。

(4)

式中,J為材料柔度,δ為近場半徑,h為二維板厚度,A為一維桿的橫截面積。

通過改變式(4)中材料柔度的取值,可以將線彈性材料的本構力函數(shù)擴展到非線性分析。

2.3 方程求解方法

PD方法的運動方程式(2)可離散為

(5)

其中,n為時間步數(shù)。

可見,近場動力學方法給出的運動方程是動力學形式的,用它來計算靜力學問題時還需做一定的處理,如可采用引入人工阻尼的動態(tài)松弛方法來求解。Underwood[13]提出了一種自適應動態(tài)松弛法求解近場動力學運動方程,其中阻尼隨物質點的位移變化,能較快地使結果收斂。

將式(5)引入阻尼cn,并改寫成矩陣形式

(6)

由中心差分法得到

(7)

初始條件為

(8)

式中,F為物質點所受合力,D為密度矩陣,其對角元素滿足

(9)

式中,e為x,y,z方向的單位矢量,cn為阻尼系數(shù)

(10)

(11)

動態(tài)松弛法是通過添加人工阻尼,從而求得函數(shù)的靜態(tài)解的一種方法。阻尼越大,收斂也就越快,但是人工阻尼的大小不能超過臨界阻尼,否則會造成計算時間過長。

3 非線性粘彈體的PD蠕變本構力函數(shù)

3.1 非線性粘彈體的時間-應力等效原理

瀝青、混凝土、固體推進劑、高聚物以及高聚物粘結炸藥等在溫度相同的條件下,應力水平越高,材料的蠕變應變就越大,材料呈現(xiàn)出非線性粘彈性特征。此時,不能只考慮時間和溫度,還需要考慮應力水平對蠕變行為的影響。非線性粘彈體的時間-應力等效原理[14]認為,材料受載應力水平對蠕變柔度的影響與溫度相似,也具有等效性。依據(jù)自由體積理論,推導出了時間-應力等效原理的表達式[15]：

J(σ,t)=bσJ(σ0,t/aσ)

(12)

式中,aσ和bσ分別為應力水平和豎直移位因子,具有與溫度移位因子類似的形式,J(σ0,t)為參考應力σ0下的蠕變柔度主曲線。

3.2 蠕變柔度主曲線的表征

粘彈性材料同時具有彈性和粘性特征,根據(jù)流變學理論采用彈性和粘性元件組合描述其粘彈性行為。粘性元件與彈性元件常見的組合模型包括Maxwell模型、Kelvin模型、Burgers模型以及其它復雜模型。一般而言,材料模型的選擇和確定應該遵循以下原則: (1)模型能夠很好地反映材料的力學性能; (2)模型應盡可能簡單、直觀。Burgers粘彈性模型是由Maxwell單元和Kelvin單元串聯(lián)組成的四參數(shù)模型,可以表示高聚物粘結炸藥粘彈行為的主要特征[16]。由實驗測得的高聚物粘結炸藥材料的蠕變曲線與Burgers模型一致性很好?；谝陨蠋c,本研究選擇了Burgers粘彈性模型來描述高聚物粘結炸藥的粘彈性行為。

參考應力σ0下的蠕變柔度主曲線可以通過Burgers粘彈性模型來描述,如圖2所示。

圖2Burgers粘彈性模型

Fig.2Burgers viscoelastic model

根據(jù)胡克定律和牛頓流體定律可得到,當應力為常數(shù)時Burgers模型的應變表達式為

(13)

式中,E1,E2為彈簧的彈性模量,MPa;η1,η2為粘壺的粘滯系數(shù),MPa;t為時間,s。

令蠕變柔度J(σ0,t)=ε(σ0,t)/σ0,則:

(14)

由此,畫出Burgers模型的蠕變曲線,如圖3所示。

圖3Burgers模型的蠕變曲線

Fig.3Creep curve of Burgers model

3.3 PD蠕變本構力函數(shù)

結合式(12)和式(14)可以得到不同應力水平下,非線性粘彈體的蠕變柔度表達式

(15)

當時間t為定值時,材料的蠕變柔度保持不變。因此,可以將式(15)代入式(4),再代入式(3),得到蠕變各個時刻對應的各參數(shù)均具有明確物理意義的PD方法本構力函數(shù)

(16)

利用Fortran自編程序實現(xiàn)算法,首先將結構離散為均勻分布的物質點,各個時刻物質點間的相互作用力用式(16)表示,采用動態(tài)松弛法求得物質點的位移。最后將該方法應用于PBX9502材料的蠕變行為模擬。

4 應用實例

4.1 應力移位因子的確定

Gagliardi[17]針對PBX9502開展了不同溫度、不同應力下的圓柱體單軸壓縮蠕變實驗,實驗結果見圖4。通過公式J=ε/σ,并取對數(shù),得到對數(shù)蠕變柔度隨時間的變化曲線,如圖5所示?？梢钥闯?24 ℃時,1.7237,3.4475 MPa低應力水平下的蠕變柔度曲線幾乎完全重合,PBX9502呈現(xiàn)出了線性粘彈性特性。隨著溫度的上升,應力水平越高,材料的蠕變應變就越大,PBX9502呈現(xiàn)出非線性粘彈性。同時,不同應力下的蠕變柔度曲線具有相似性,這與溫度變化時的蠕變行為類似,即時間和應力對PBX9502蠕變行為的影響也具有等效性。

圖4PBX9502的蠕變實驗結果

Fig.4Creep test results of PBX9502

a. 24 ℃

b. 50 ℃

c. 70 ℃

圖5不同應力下PBX9502的對數(shù)蠕變柔度曲線

Fig.5Logarithmic creep compliance curves of PBX9502 under different stresses

以σ0為參考應力,將其對應的對數(shù)蠕變柔度主曲線向其他應力水平下的蠕變曲線做相應的豎直移位,使兩曲線重合,得到的蠕變曲線如圖6所示。圖6中顯示,經(jīng)過移位的對數(shù)蠕變主曲線與其他應力水平下的蠕變曲線幾乎完全重合,這是由于不同應力下的蠕變柔度曲線具有相似性。因此,應力水平移位因子aσ=1。

a. 24 ℃

b. 50 ℃

c. 70 ℃

圖6豎直移位后PBX9502的對數(shù)蠕變柔度曲線

Fig.6Logarithmic creep compliance curves of PBX9502 after vertical translation

另外,從式(12),可以看出豎直移位因子bσ等于不同應力下初始柔度的比值。

根據(jù)如圖4所示的PBX9502的蠕變試驗數(shù)據(jù),擬合得到Burgers模型參數(shù),結果見圖7和表1。

圖7Burgers模型參數(shù)擬合

Fig.7Parameter fitting of Burgers model

表1PBX9502的Burgers模型參數(shù)

Table1Parameters of Burgers model for PBX9502

T/℃σ0/MPaE1/GPaη1/GPaE2/GPaη2/GPa245．32294．83903．9394×1076．40377．4478×105501．72372．85382．6254×1074．71644．1235×105703．44752．15471．7862×1071．23570．9212×105

Note：σ0is reference stress,E1,E2are elastic modulas of the spring,η1,η2are cofficients of viscosity of the dashpot.

通過蠕變柔度主曲線與其他應力水平下蠕變柔度曲線初始柔度的比值,得到的應力豎直移位因子bσ的取值列于表2。

表2應力豎直移位因子bσ的取值

Table2Values of stress vertical translation factorbσ

T/℃σ0stress/MPaσ=1．7237σ=3．4475σ=5．3229245．32290．80830．7697501．72371．09561．3352703．44750．74881．2600

4.2 PBX圓柱體受壓蠕變行為的PD模擬

如圖8所示,蠕變實驗采用圓柱體試樣[17],直徑12.7 mm,高25.4 mm,上下受均勻分布的壓應力。模型被離散為均勻的物質點,計算參數(shù)如下: dt=1.0 s,密度ρ=1900 kg/m3,泊松比ν=0.25,物質點間距dx=0.6 mm,近場半徑δ=3dx。

圖9給出初始應變和應力隨載荷步的變化趨勢,可以看出,PBX炸藥柱的應力和應變約在150載荷步達到穩(wěn)定。

a. geometric modelb. discrete model

圖8數(shù)值計算模型

Fig.8Numerical simulation model

a. change of initial strain with time step

b. change of stress with time step

圖9應力和應變隨載荷步的變化

Fig.9The change of stress and strain with time step

將參考應力σ0對應蠕變柔度主曲線的Burgers模型參數(shù)以及應力豎直移位因子bσ代入式(16)的PD蠕變本構力函數(shù),得到的模擬蠕變曲線與實驗曲線的對比情況見圖10,可以看出,模擬曲線與實驗曲線吻合很好。

a. 24 ℃

b. 50 ℃

c. 70 ℃

圖10PBX9502模擬蠕變曲線與試驗曲線的對比

Fig.10Comparison of simulated and experimental creep curves of PBX9502

5 結 論

近場動力學方法是一種新的無網(wǎng)格方法,它在分析損傷、斷裂和失穩(wěn)等不連續(xù)問題時具有優(yōu)勢。但是,作為一種新的無網(wǎng)格方法,PD方法還在發(fā)展中,在粘彈性、塑性、損傷斷裂與破壞等問題分析上還有待發(fā)展,需要做進一步的研究。

本研究通過非線性粘彈體的時間-應力等效原理,與Burgers粘彈性模型得到了不同應力作用下材料蠕變柔度的表達式,推導出非線性粘彈體的PD蠕變本構力函數(shù),從而建立起可應用于高聚物粘結炸藥的近場動力學蠕變模擬方法。模擬了PBX9502在溫度和應力作用下的蠕變行為,獲得與實驗一致的結果。

本研究建立的近場動力學蠕變行為模擬方法可應用于同時計及溫度和應力作用的高聚物粘結炸藥的蠕變行為分析。

參考文獻:

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

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

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

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

[4] 楊挺青. 粘彈性力學[M]. 武漢: 華中理工大學出版社, 1992: 8-22.

YANG Ting-qing. Viscoelastic mechanics[M]. Wuhan: Press of Huazhong University of Science and Technology, 1992: 8-22.

[5] 周光泉, 劉孝敏. 粘彈性理論[M]. 合肥 : 中國科學技術大學出版社, 1996: 8-15.

ZHOU Guang-quan, LIU Xiao-min. Viscoelastic theory[M]. Hefei: Press of University of Science and Technology of China, 1996: 8-15.

[6] 唐明峰, 藍林鋼, 李明, 等. 以RDX為基的澆注PBX力學性能與本構模型[J]. 含能材料, 2014, 22(2): 215-220.

TANG Ming-feng, LAN Lin-gang, LI Ming, et al.Mechanical properties and constitutive models of RDX based cast PBX[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2014, 22(2): 215-220.

[7] 孟紅磊, 鞠玉濤. 含損傷非線性粘彈性本構模型及數(shù)值仿真應用[J]. 固體火箭技術, 2012, 35(6): 764-768.

MENG Hong-lei, JU Yu-tao.Nonlinear viscoelastic equation with cumulative damage and its application on numerical simulation[J].JournalofSolidRocketTechnology, 2012, 35(6): 764-768.

[8] 馮震宇, 王新軍, 王富生, 等. 朱-王-唐非線性粘彈性本構模型在有限元分析中的實現(xiàn)及其應用[J]. 材料科學與工程學報, 2007, 25(2): 269-272.

FENG Zhen-yu, WANG Xin-jun, WANG Fu-sheng, et al. Implementation and its application in finite element analysis of constitutive model for ZWT nonlinear viscoelastic material[J].JournalofMaterialsScience&Engineering, 2007, 25(2): 269-272.

[9] Mitchell J A. A non-local, ordinary-state-based viscoelasticity model for peridynamics. SAND2011-8064, 2011.

[10] Azizi M A, Ariffin A K, Nik M I, et al. The peridynamic model of viscoelastic creep and recovery[J].MultidisciplineModelinginMaterialsandStructures, 2015, 11(4): 579-597.

[11] Silling S A, Askari E. A meshfree method on the peridynamic model of solid mechanics[J].ComputersandStructures, 2005, 83: 1526-1535.

[12] Madenci E, Oterkus E. Peridynamic theory and its applications[M]. New York: Springer Science and Business Media, 2014.

[13] Underwood P. Dynamic relaxation[J].ComputMethTransAnal, 1983, 1: 245-265.

[14] Jazouli S, Wenbo Luo, Bremand F, Vu-Khanh T. Application of time-stress equivalence to nonlinear creep of polycarbonate[J].PolymerTesting, 2005, 24: 463-467.

[15] 王志方, 張國忠. 膠凝原油蠕變的時間-溫度-應力等效性[J]. 力學與實踐, 2008, 30(1): 62-65.

WANG Zhi-fang, ZHANG Guo-zhong.The time-temperature-stress equivalence of creep behavior of gelled crude oil[J].MechanicsinEngineering, 2008, 30(1): 62-65.

[16] 蔡峨. 粘彈性力學基礎[M]. 北京: 北京航空航天出版社, 1989: 54-81.

CAI E. Basis of viscoelastic mechanics[M]. Beijing: Press of Beijing University of Aeronautics and Astronautics, 1989: 54-81.

[17] Gagliardi F J, Cunningham B J. Creep testing plastic-bonded explosives in uni-axial compression[D]. United States: SEM XI international congress Orlando, FL, 2008.