高聚物粘結(jié)炸藥蠕變行為的近場(chǎng)動(dòng)力學(xué)模擬方法

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

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

1 引 言

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

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

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

本研究利用Burgers粘彈性模型表征蠕變?nèi)岫戎髑€(xiàn),結(jié)合非線(xiàn)性粘彈體的時(shí)間-應(yīng)力等效原理,得到不同應(yīng)力水平下蠕變?nèi)岫鹊谋磉_(dá)式。將其代入PD本構(gòu)力函數(shù),推導(dǎo)出非線(xiàn)性粘彈體的PD蠕變本構(gòu)力函數(shù),從而建立起一種可應(yīng)用于高聚物粘結(jié)炸藥的近場(chǎng)動(dòng)力學(xué)蠕變模擬方法。利用該方法對(duì)PBX9502在不同溫度和不同應(yīng)力作用下的蠕變行為進(jìn)行模擬,獲得與實(shí)驗(yàn)一致的結(jié)果。

2 PD方法的基本理論

2.1 運(yùn)動(dòng)方程

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

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

(1)

根據(jù)牛頓第二定律可得物質(zhì)點(diǎn)xk的運(yùn)動(dòng)方程為

b(xk,t)

(2)

式中,ρ為物質(zhì)點(diǎn)密度,u為物質(zhì)點(diǎn)的位移,b為單位體積物質(zhì)所受的外載荷。

圖1物質(zhì)點(diǎn)與其近場(chǎng)范圍內(nèi)其它物質(zhì)點(diǎn)的相互作用

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

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

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

(3)

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

(4)

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

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

2.3 方程求解方法

PD方法的運(yùn)動(dòng)方程式(2)可離散為

(5)

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

可見(jiàn),近場(chǎng)動(dòng)力學(xué)方法給出的運(yùn)動(dòng)方程是動(dòng)力學(xué)形式的,用它來(lái)計(jì)算靜力學(xué)問(wèn)題時(shí)還需做一定的處理,如可采用引入人工阻尼的動(dòng)態(tài)松弛方法來(lái)求解。Underwood[13]提出了一種自適應(yīng)動(dòng)態(tài)松弛法求解近場(chǎng)動(dòng)力學(xué)運(yùn)動(dòng)方程,其中阻尼隨物質(zhì)點(diǎn)的位移變化,能較快地使結(jié)果收斂。

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

(6)

由中心差分法得到

(7)

初始條件為

(8)

式中,F為物質(zhì)點(diǎn)所受合力,D為密度矩陣,其對(duì)角元素滿(mǎn)足

(9)

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

(10)

(11)

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

3 非線(xiàn)性粘彈體的PD蠕變本構(gòu)力函數(shù)

3.1 非線(xiàn)性粘彈體的時(shí)間-應(yīng)力等效原理

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

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

(12)

式中,aσ和bσ分別為應(yīng)力水平和豎直移位因子,具有與溫度移位因子類(lèi)似的形式,J(σ0,t)為參考應(yīng)力σ0下的蠕變?nèi)岫戎髑€(xiàn)。

3.2 蠕變?nèi)岫戎髑€(xiàn)的表征

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

參考應(yīng)力σ0下的蠕變?nèi)岫戎髑€(xiàn)可以通過(guò)Burgers粘彈性模型來(lái)描述,如圖2所示。

圖2Burgers粘彈性模型

Fig.2Burgers viscoelastic model

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

(13)

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

令蠕變?nèi)岫菾(σ0,t)=ε(σ0,t)/σ0,則:

(14)

由此,畫(huà)出Burgers模型的蠕變曲線(xiàn),如圖3所示。

圖3Burgers模型的蠕變曲線(xiàn)

Fig.3Creep curve of Burgers model

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

結(jié)合式(12)和式(14)可以得到不同應(yīng)力水平下,非線(xiàn)性粘彈體的蠕變?nèi)岫缺磉_(dá)式

(15)

當(dāng)時(shí)間t為定值時(shí),材料的蠕變?nèi)岫缺３植蛔?。因?可以將式(15)代入式(4),再代入式(3),得到蠕變各個(gè)時(shí)刻對(duì)應(yīng)的各參數(shù)均具有明確物理意義的PD方法本構(gòu)力函數(shù)

(16)

利用Fortran自編程序?qū)崿F(xiàn)算法,首先將結(jié)構(gòu)離散為均勻分布的物質(zhì)點(diǎn),各個(gè)時(shí)刻物質(zhì)點(diǎn)間的相互作用力用式(16)表示,采用動(dòng)態(tài)松弛法求得物質(zhì)點(diǎn)的位移。最后將該方法應(yīng)用于PBX9502材料的蠕變行為模擬。

4 應(yīng)用實(shí)例

4.1 應(yīng)力移位因子的確定

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

圖4PBX9502的蠕變實(shí)驗(yàn)結(jié)果

Fig.4Creep test results of PBX9502

a. 24 ℃

b. 50 ℃

c. 70 ℃

圖5不同應(yīng)力下PBX9502的對(duì)數(shù)蠕變?nèi)岫惹€(xiàn)

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

以σ0為參考應(yīng)力,將其對(duì)應(yīng)的對(duì)數(shù)蠕變?nèi)岫戎髑€(xiàn)向其他應(yīng)力水平下的蠕變曲線(xiàn)做相應(yīng)的豎直移位,使兩曲線(xiàn)重合,得到的蠕變曲線(xiàn)如圖6所示。圖6中顯示,經(jīng)過(guò)移位的對(duì)數(shù)蠕變主曲線(xiàn)與其他應(yīng)力水平下的蠕變曲線(xiàn)幾乎完全重合,這是由于不同應(yīng)力下的蠕變?nèi)岫惹€(xiàn)具有相似性。因此,應(yīng)力水平移位因子aσ=1。

a. 24 ℃

b. 50 ℃

c. 70 ℃

圖6豎直移位后PBX9502的對(duì)數(shù)蠕變?nèi)岫惹€(xiàn)

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

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

根據(jù)如圖4所示的PBX9502的蠕變?cè)囼?yàn)數(shù)據(jù),擬合得到Burgers模型參數(shù),結(jié)果見(jiàn)圖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.

通過(guò)蠕變?nèi)岫戎髑€(xiàn)與其他應(yīng)力水平下蠕變?nèi)岫惹€(xiàn)初始柔度的比值,得到的應(yīng)力豎直移位因子bσ的取值列于表2。

表2應(yīng)力豎直移位因子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所示,蠕變實(shí)驗(yàn)采用圓柱體試樣[17],直徑12.7 mm,高25.4 mm,上下受均勻分布的壓應(yīng)力。模型被離散為均勻的物質(zhì)點(diǎn),計(jì)算參數(shù)如下: dt=1.0 s,密度ρ=1900 kg/m3,泊松比ν=0.25,物質(zhì)點(diǎn)間距dx=0.6 mm,近場(chǎng)半徑δ=3dx。

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

a. geometric modelb. discrete model

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

Fig.8Numerical simulation model

a. change of initial strain with time step

b. change of stress with time step

圖9應(yīng)力和應(yīng)變隨載荷步的變化

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

將參考應(yīng)力σ0對(duì)應(yīng)蠕變?nèi)岫戎髑€(xiàn)的Burgers模型參數(shù)以及應(yīng)力豎直移位因子bσ代入式(16)的PD蠕變本構(gòu)力函數(shù),得到的模擬蠕變曲線(xiàn)與實(shí)驗(yàn)曲線(xiàn)的對(duì)比情況見(jiàn)圖10,可以看出,模擬曲線(xiàn)與實(shí)驗(yàn)曲線(xiàn)吻合很好。

a. 24 ℃

b. 50 ℃

c. 70 ℃

圖10PBX9502模擬蠕變曲線(xiàn)與試驗(yàn)曲線(xiàn)的對(duì)比

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

5 結(jié) 論

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

本研究通過(guò)非線(xiàn)性粘彈體的時(shí)間-應(yīng)力等效原理,與Burgers粘彈性模型得到了不同應(yīng)力作用下材料蠕變?nèi)岫鹊谋磉_(dá)式,推導(dǎo)出非線(xiàn)性粘彈體的PD蠕變本構(gòu)力函數(shù),從而建立起可應(yīng)用于高聚物粘結(jié)炸藥的近場(chǎng)動(dòng)力學(xué)蠕變模擬方法。模擬了PBX9502在溫度和應(yīng)力作用下的蠕變行為,獲得與實(shí)驗(yàn)一致的結(jié)果。

本研究建立的近場(chǎng)動(dòng)力學(xué)蠕變行為模擬方法可應(yīng)用于同時(shí)計(jì)及溫度和應(yīng)力作用的高聚物粘結(jié)炸藥的蠕變行為分析。

參考文獻(xiàn):

[1] 黃丹, 章青, 喬丕忠, 等. 近場(chǎng)動(dòng)力學(xué)方法以及應(yīng)用[J]. 力學(xué)進(jìn)展, 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] 楊挺青. 粘彈性力學(xué)[M]. 武漢: 華中理工大學(xué)出版社, 1992: 8-22.

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

[5] 周光泉, 劉孝敏. 粘彈性理論[M]. 合肥 : 中國(guó)科學(xué)技術(shù)大學(xué)出版社, 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] 唐明峰, 藍(lán)林鋼, 李明, 等. 以RDX為基的澆注PBX力學(xué)性能與本構(gòu)模型[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] 孟紅磊, 鞠玉濤. 含損傷非線(xiàn)性粘彈性本構(gòu)模型及數(shù)值仿真應(yīng)用[J]. 固體火箭技術(shù), 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)性粘彈性本構(gòu)模型在有限元分析中的實(shí)現(xiàn)及其應(yīng)用[J]. 材料科學(xué)與工程學(xué)報(bào), 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] 王志方, 張國(guó)忠. 膠凝原油蠕變的時(shí)間-溫度-應(yīng)力等效性[J]. 力學(xué)與實(shí)踐, 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] 蔡峨. 粘彈性力學(xué)基礎(chǔ)[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.