含能材料變泊松比粘彈性本構及其數(shù)值實現(xiàn)方法

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

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

0 引言

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

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

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

1 粘彈性本構模型

1.1 粘彈性時間相關泊松比形式

(1)

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

(2)

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

1.2 采用時間相關泊松比的粘彈性本構方程

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

(3)

其中

(4)

根據彈性-粘彈性對應原理，得到式(3)、式(4)在相域內的對應線粘彈性關系式，再對其作Laplace逆變換，并基于熱流變簡單材料假設，可得時域內的線熱粘彈性本構方程：

(5)

其中

(6)

(7)

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

(8)

(9)

(10)

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

(11)

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

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

2.1 本構方程的增量形式

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

(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)

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

(16)

其中

(17)

(18)

(19)

(20)

(21)

(22)

(23)

其中

(24)

(25)

(26)

應用中值定理，對式(26)進行積分計算，可得

(27)

其中

(28)

(29)

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

(30)

(31)

(32)

(33)

(34)

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

(35)

其中

(36)

(37)

(38)

(39)

(40)

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

(41)

(42)

(43)

(44)

(45)

其中

(46)

(47)

(48)

2.2 本構方程的切線剛度

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

(49)

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

(50)

(51)

(52)

3 算例分析

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

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

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

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

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

3.1 驗證算例：粘彈性薄板的應力松弛分析

采用變泊松比線性粘彈性本構方程進行分析。為模擬單軸應力松弛試驗過程，采用尺寸為10 mm×50 mm×5 mm的方形薄板模型，其有限元網格模型見圖1。

圖1 單向拉伸有限元模型

(53)

(54)

(55)

將式(54)兩邊對t求導并結合式(55)，可得其橫向應εy(t)變表達式：

(56)

圖2 縱向應力沿隨時間變化對比曲線

圖3 橫向應變沿隨時間變化對比曲線

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

3.2 合理性算例：內壓載荷下的裝藥結構分析

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

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

以受均勻內壓載荷下的細長固體火箭發(fā)動機為例，模型的幾何尺寸為內徑200 mm，外徑397 mm，取其一段所建計算模型見圖4。

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

圖4 藥柱有限元模型

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

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

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

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

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

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

4 結論

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

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

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