鄧 斌，賀元吉，趙宏偉，江增榮，余慶波

(1. 中國(guó)人民解放軍96901部隊(duì)，北京 100094；2. 北京理工大學(xué)，北京 100081)

0 引言

長(zhǎng)期以來(lái)，因限于認(rèn)知、測(cè)量手段等因素，通常為簡(jiǎn)化問(wèn)題人們往往將火炸藥、推進(jìn)劑等粘彈性材料的泊松比視作常數(shù)處理[1-5]。此外，針對(duì)火炸藥、推進(jìn)劑等粘彈性材料的泊松比測(cè)量，目前國(guó)內(nèi)現(xiàn)有行業(yè)標(biāo)準(zhǔn)如QJ 3228—2005[6]和GJB 770B—2005[7]，均將其泊松比視作常數(shù)進(jìn)行處理，并采用了彈性泊松比試驗(yàn)測(cè)量手段。

大量研究表明[8-15]，實(shí)際的粘彈性材料泊松比為時(shí)間或頻率的函數(shù)，且與溫度緊密相關(guān)。采用泊松比作為計(jì)算參數(shù)進(jìn)行結(jié)構(gòu)分析時(shí)，泊松比的微小變化可導(dǎo)致計(jì)算結(jié)果的顯著差異，尤其是對(duì)于火炸藥、固體推進(jìn)劑等近似不可壓粘彈性材料，以往將其視作常數(shù)的處理方式，在某些情況下會(huì)導(dǎo)致計(jì)算結(jié)果的較大誤差，尤其是對(duì)于某些特定條件如短時(shí)沖擊載荷以及變溫過(guò)程，或當(dāng)粘彈性材料處于狀態(tài)轉(zhuǎn)變區(qū)時(shí)，粘彈性材料泊松比的變化往往較大[11-12]。此時(shí)用定泊松比模型已無(wú)法準(zhǔn)確表征這種變化往往導(dǎo)致分析結(jié)果的較大誤差[12]，難以準(zhǔn)確反映粘彈性材料真實(shí)力學(xué)行為。隨著人們對(duì)粘彈性結(jié)構(gòu)分析精度要求的不斷提高，與時(shí)間相關(guān)粘彈性泊松比的研究因而受到重視。目前，關(guān)于粘彈性泊松比理論和試驗(yàn)研究的成果也較多[9-16]，粘彈性問(wèn)題計(jì)算方法研究也主要集中于定泊松比模型[17-18]，而關(guān)于變化泊松比粘彈性結(jié)構(gòu)數(shù)值分析方法及其應(yīng)用研究則鮮見(jiàn)報(bào)道。

為此，考慮到粘彈性材料泊松比的時(shí)間相關(guān)性以及泊松比參數(shù)變化對(duì)計(jì)算結(jié)果的顯著影響，本文采用一種時(shí)間相關(guān)粘彈性材料泊松比形式建立三維線性粘彈性本構(gòu)模型，研究相應(yīng)工程應(yīng)用的數(shù)值分析方法，旨在為火炸藥、推進(jìn)劑等含能材料裝藥精細(xì)結(jié)構(gòu)完整性分析提供支持。

1 粘彈性本構(gòu)模型

1.1 粘彈性時(shí)間相關(guān)泊松比形式

(1)

一般地，縱向應(yīng)變并非階躍形式，而是隨時(shí)間的連續(xù)變化量。若在任意τ時(shí)刻施加一縱向階躍微應(yīng)變dεx(τ)，則持續(xù)作用至t時(shí)刻的橫向應(yīng)變響應(yīng)為dεy(t)=-ν(t-τ)dεx(τ)，當(dāng)τ在整個(gè)[0,t]上連續(xù)變化時(shí)，根據(jù)線性系統(tǒng)的Boltzmann疊加原理，可得t時(shí)刻的總橫向應(yīng)變響應(yīng)：

(2)

可知，橫向應(yīng)變?chǔ)舮(t)要滯后于縱向應(yīng)變?chǔ)舩(t)響應(yīng)歷程，且粘彈性泊松比是橫向變形的記憶函數(shù)，而非簡(jiǎn)單的代數(shù)關(guān)系。為區(qū)分傳統(tǒng)的定泊松比 ，后續(xù)也稱時(shí)間相關(guān)泊松比ν(t)為變泊松比。

1.2 采用時(shí)間相關(guān)泊松比的粘彈性本構(gòu)方程

對(duì)于各向同性線彈性材料，采用楊氏模量E和泊松比ν表示的本構(gòu)形式如下：

(3)

其中

(4)

根據(jù)彈性-粘彈性對(duì)應(yīng)原理，得到式(3)、式(4)在相域內(nèi)的對(duì)應(yīng)線粘彈性關(guān)系式，再對(duì)其作Laplace逆變換，并基于熱流變簡(jiǎn)單材料假設(shè)，可得時(shí)域內(nèi)的線熱粘彈性本構(gòu)方程：

(5)

其中

(6)

(7)

其中，泊松比ν(t)和松弛模量E(t)Prony級(jí)數(shù)形式分別為

(8)

(9)

(10)

其中，aT為溫度平移因子，滿足如下WLF方程：

(11)

式中C1、C2為材料常數(shù)；T為當(dāng)前溫度；T0為參考溫度。

2 本構(gòu)方程的數(shù)值模型

2.1 本構(gòu)方程的增量形式

針對(duì)式(5)所示的本構(gòu)方程，采用增量有限元法對(duì)其進(jìn)行數(shù)值離散，給出其增量形式。將分析時(shí)間[0,t]劃分為[0,t1],[t1,t2], …,[tm,tm+1],…,[tM-1,tM]共M個(gè)子時(shí)間增量步，則可將其在任意增量步[tm,tm+1]內(nèi)的增量形式寫成：

(12)

其中

Δσij(tm+1)=σij(tm+1)-σij(tm)

(13)

ΔSij(tm+1)=Sij(tm+1)-Sij(tm)

(14)

Δσkk(tm+1)=σkk(tm+1)-σkk(tm)

(15)

針對(duì)式(6)所示本構(gòu)方程偏張量部分，結(jié)合Stieltjes卷積定理，可得其在[tm,tm+1]內(nèi)的增量形式：

(16)

其中

(17)

(18)

(19)

(20)

(21)

(22)

(23)

其中

(24)

(25)

(26)

應(yīng)用中值定理，對(duì)式(26)進(jìn)行積分計(jì)算，可得

(27)

其中

(28)

(29)

根據(jù)式(8)和式(9)，式(18)和式(20)可寫成如下形式：

(30)

(31)

(32)

(33)

(34)

對(duì)于式(7)所示本構(gòu)部分，類似式(6)處理方法，可得其在[tm,tm+1]內(nèi)增量形式：

(35)

其中

(36)

(37)

(38)

(39)

(40)

針對(duì)式(36)～式(40)，可采用類似式(17)～式(21)的方法進(jìn)行數(shù)值離散，得到式(36) ～式(40)的表達(dá)式：

(41)

(42)

(43)

(44)

(45)

其中

(46)

(47)

(48)

2.2 本構(gòu)方程的切線剛度

利用一致切線剛度可保證有限元分析所采用的Newton法具有二階收斂速率。根據(jù)切線剛度的定義，對(duì)于小變形或者小體變的大變形問(wèn)題，其一致切線剛度的形式為

(49)

聯(lián)立式(12)和式(49)，可得在第tm+1時(shí)刻的切線剛度張量各個(gè)分量：

(50)

(51)

(52)

3 算例分析

針對(duì)上述變泊松比粘彈性本構(gòu)增量模型，基于Abaqus軟件平臺(tái)二次開發(fā)技術(shù)，采用Fortran語(yǔ)言編寫相應(yīng)UMAT材料子程序，并依托軟件平臺(tái)前后處理和求解功能，實(shí)現(xiàn)該本構(gòu)的工程應(yīng)用。

以某HTPB推進(jìn)劑為例，取式(11)中的C1=4.97，C2=156.1，T0=293.15 K；式(9)所示松弛模量取4項(xiàng)時(shí)的各參數(shù)見(jiàn)表1；根據(jù)某型推進(jìn)劑泊松比試驗(yàn)數(shù)據(jù)，擬合得到式(8)所示 取5項(xiàng)時(shí)的相關(guān)系數(shù)見(jiàn)表2；其他部件材料性能參數(shù)見(jiàn)表3。

表1 松弛模量Prony級(jí)數(shù)系數(shù)

表2 泊松比Prony級(jí)數(shù)系數(shù)

表3 其他材料參數(shù)

3.1 驗(yàn)證算例：粘彈性薄板的應(yīng)力松弛分析

采用變泊松比線性粘彈性本構(gòu)方程進(jìn)行分析。為模擬單軸應(yīng)力松弛試驗(yàn)過(guò)程，采用尺寸為10 mm×50 mm×5 mm的方形薄板模型，其有限元網(wǎng)格模型見(jiàn)圖1。

圖1 單向拉伸有限元模型

(53)

(54)

(55)

將式(54)兩邊對(duì)t求導(dǎo)并結(jié)合式(55)，可得其橫向應(yīng)εy(t)變表達(dá)式：

(56)

圖2 縱向應(yīng)力沿隨時(shí)間變化對(duì)比曲線

圖3 橫向應(yīng)變沿隨時(shí)間變化對(duì)比曲線

從圖2和圖3中可看出，采用變泊松比粘彈性本構(gòu)模型進(jìn)行分析時(shí)，本文有限元解與對(duì)應(yīng)解析解吻合良好，且具有很高的計(jì)算精度，這進(jìn)一步驗(yàn)證了本文算法和程序的有效性。

3.2 合理性算例：內(nèi)壓載荷下的裝藥結(jié)構(gòu)分析

相比于定泊松比本構(gòu)模型，同等條件下的變泊松比粘彈性本構(gòu)模型相關(guān)問(wèn)題的求解要復(fù)雜得多，一般難以給出其解析解，且因限于條件暫無(wú)法開展相關(guān)試驗(yàn)驗(yàn)證。下面就算法程序做一些合理性驗(yàn)證分析，其基本思路：將變泊松比粘彈性本構(gòu)模型所得結(jié)果與相應(yīng)的定泊松比本構(gòu)模型所得結(jié)果進(jìn)行比較，并根據(jù)力學(xué)規(guī)律來(lái)判斷結(jié)果的合理性來(lái)驗(yàn)證算法程序合理性。

根據(jù)式(8)及表2中的系數(shù)可知，該推進(jìn)劑泊松比為時(shí)間遞增函數(shù)，其初始值約為0.487，在0.3 s內(nèi)其值位于0.487～0.491之間，若本文算法程序有效，則對(duì)于變泊松比模型所得結(jié)果，應(yīng)位于定泊松比為0.487和0.491時(shí)所對(duì)應(yīng)結(jié)果之間，且開始時(shí)刻結(jié)果接近于0.487對(duì)應(yīng)值，而在0.3 s時(shí)結(jié)果接近于0.491對(duì)應(yīng)值。

以受均勻內(nèi)壓載荷下的細(xì)長(zhǎng)固體火箭發(fā)動(dòng)機(jī)為例，模型的幾何尺寸為內(nèi)徑200 mm，外徑397 mm，取其一段所建計(jì)算模型見(jiàn)圖4。

在同等條件下進(jìn)行分析，當(dāng)本構(gòu)(5)退化至定泊松比線粘彈性本構(gòu)時(shí)，所得結(jié)果與Marc線粘彈性本構(gòu)解吻合良好，限于篇幅，此處不詳述。在其內(nèi)表面施加均布?jí)毫d荷P(t)=6(1-e-20t)MPa，計(jì)算時(shí)間t=0.3 s，分析時(shí)取定泊松比ν=0.487和0.491分別進(jìn)行計(jì)算，采用有限元網(wǎng)格模型、載荷及位移邊界條件，先后考慮本文變泊松比和定泊松比(取為0.487和0.491)粘彈性本構(gòu)模型，分別對(duì)上述問(wèn)題進(jìn)行分析，其中該問(wèn)題在采用定泊松比時(shí)所得環(huán)向、徑向的應(yīng)力應(yīng)變關(guān)系式見(jiàn)文獻(xiàn)[19]。得到圓管內(nèi)表面的應(yīng)力、應(yīng)變隨時(shí)間變化曲線分別見(jiàn)圖5和圖6所示；加載至0.3 s時(shí)的應(yīng)力、應(yīng)變沿徑向變化規(guī)律分別如圖7和圖8所示。

圖4 藥柱有限元模型

(a)徑向應(yīng)力 (b)環(huán)向應(yīng)力

(a)徑向應(yīng)變 (b)環(huán)向應(yīng)變

由圖5和圖6可知，在加載過(guò)程的任意時(shí)刻，變泊松比所對(duì)應(yīng)的應(yīng)力應(yīng)變結(jié)果，均位于定泊松比取0.487和0.491時(shí)所對(duì)應(yīng)結(jié)果之間，且在開始階段，變泊松比所對(duì)應(yīng)結(jié)果接近于定泊松比為0.487時(shí)的對(duì)應(yīng)結(jié)果，而隨著時(shí)間的增加，逐漸向定泊松比取0.491時(shí)的對(duì)應(yīng)結(jié)果靠近，在加載至0.3 s時(shí)，其各部位的應(yīng)力應(yīng)變結(jié)果均已非常接近于定泊松比取0.491時(shí)的對(duì)應(yīng)結(jié)果(圖7和圖8)，這與變泊松比值從0.487變化到0.491趨勢(shì)是一致的，符合力學(xué)變化規(guī)律。

(a)徑向應(yīng)力 (b)環(huán)向應(yīng)力

(a)徑向應(yīng)變 (b)環(huán)向應(yīng)變

綜上，隨著分析時(shí)間的加長(zhǎng)，變泊松比值逐漸趨于接近平衡值，可以預(yù)測(cè)，對(duì)于長(zhǎng)時(shí)間持續(xù)加載分析過(guò)程，采用取平衡泊松比值的定泊松比參數(shù)與變泊松比參數(shù)所得對(duì)應(yīng)結(jié)果間幾乎無(wú)差異，因而此時(shí)，在缺乏變泊松比參數(shù)時(shí)，可采用取平衡泊松比值的定泊松比粘彈性本構(gòu)模型進(jìn)行分析。然而，對(duì)于泊松比值變化較大的過(guò)程，如短時(shí)加載、劇烈變溫過(guò)程，或當(dāng)粘彈性材料處于狀態(tài)轉(zhuǎn)變區(qū)時(shí)，采用泊松比取某一定值所得結(jié)果與實(shí)際差異會(huì)較大，有必要考慮變泊松比粘彈性本構(gòu)模型進(jìn)行分析。

4 結(jié)論

(1)采用變泊松比粘彈性本構(gòu)的本文數(shù)值解與對(duì)照解析解吻合良好，合理性算例結(jié)果符合力學(xué)規(guī)律，驗(yàn)證了本文算法程序的正確性。

(2)此粘彈材料泊松比為時(shí)間遞增函數(shù)，在任意給定時(shí)間內(nèi)存在相應(yīng)上下限值，且該變泊松比模型所得結(jié)果，均位于定泊松比模型分別取該上下限值時(shí)所對(duì)應(yīng)結(jié)果之間，且其變化趨勢(shì)符合規(guī)律，驗(yàn)證了本文算法程序的合理性。

(3)利用本文算法及材料子程序方法，有效實(shí)現(xiàn)了時(shí)間相關(guān)泊松比粘彈性本構(gòu)方程的數(shù)值分析與工程應(yīng)用，可為后續(xù)火含能材料裝藥精細(xì)結(jié)構(gòu)分析提供支持。