張 曉，鄭 堅，彭 威，張前圖

(軍械工程學(xué)院，石家莊 050003)

?

HTPB復(fù)合固體推進劑粘彈性應(yīng)變能及非線性本構(gòu)模型

張 曉，鄭 堅，彭 威，張前圖

(軍械工程學(xué)院，石家莊 050003)

為了準確表征HTPB復(fù)合固體推進劑在有限變形條件下的力學(xué)性能，針對推進劑粘彈性應(yīng)變能及本構(gòu)模型進行研究。提出了推進劑粘彈性應(yīng)變能函數(shù)和非線性本構(gòu)方程的一般形式，并通過一元非線性回歸方法擬合不同應(yīng)變率下的拉伸試驗數(shù)據(jù)，得到了材料參數(shù)關(guān)于應(yīng)變率的函數(shù)，并由此建立了推進劑單軸拉伸變形下的應(yīng)變能函數(shù)和本構(gòu)方程，預(yù)測了不同應(yīng)變率下的應(yīng)力曲線，與試驗結(jié)果和已有模型的預(yù)測結(jié)果進行了對比。結(jié)果表明，材料參數(shù)與應(yīng)變率之間呈現(xiàn)冪函數(shù)關(guān)系；推進劑應(yīng)變能密度隨變形量的增大呈非線性單調(diào)增長，同一變形條件下，應(yīng)變率越高，推進劑的應(yīng)變能密度越大；本構(gòu)方程可準確描述推進劑拉伸變形的應(yīng)力應(yīng)變關(guān)系，且尤其適用于表征低應(yīng)變率下，材料在有限變形內(nèi)的粘彈特性。

HTPB推進劑；應(yīng)變能；非線性回歸；參數(shù)函數(shù)；本構(gòu)方程；變形

0 引言

端羥基聚丁二烯(HTPB)推進劑是一種以高聚物為基體的多組分、高延伸率復(fù)合材料，具有強烈的非線性粘彈特性[1]，且推進劑的力學(xué)性能主要取決于基體。針對高聚物非線性粘彈性行為的研究，Leaderman[2]最早指出非線性粘彈行為可以在小應(yīng)變時發(fā)生，并在修正Boltzmann線性疊加原理的基礎(chǔ)上給出了半經(jīng)驗公式；Ward和Onat進一步對疊加原理進行了修正，得到了非線性粘彈性的多重積分型本構(gòu)關(guān)系[3]；Shcapery[4-6]、Findley[7]、Bernstein等[8]做了大量工作，分別發(fā)展了熱力學(xué)本構(gòu)關(guān)系、冪律關(guān)系和BKZ理論。對于復(fù)合固體推進劑的線性粘彈性，國內(nèi)外已做了大量研究，并已相當(dāng)成熟；同時，國內(nèi)外在推進劑非線性力學(xué)行為方面，也做了很多卓有成效的工作，但無論理論方面、數(shù)值模擬方面，還是精確解方面，都遠不如前者充分和深入，尚有很多問題需要進一步探討。由于不滿足Boltzmann線性疊加原理，對推進劑材料的非線性粘彈性行為進行描述，要運用更加普遍的理論。Gyoo等[9]基于Vratsanos-Farris模型提出了一個含Mullins效應(yīng)和損傷脫濕效應(yīng)的三維非線性粘彈性模型；龔建良等[10]根據(jù)熱力學(xué)能量守恒定律，研究了HTPB推進劑增強粒子脫濕引起的非線性行為，確立了臨界脫濕應(yīng)變方程；鄧凱等[11]將Shcapery非線性粘彈本構(gòu)模型與Perzyna非線性粘塑本構(gòu)模型結(jié)合，建立了HTPB推進劑非線性粘彈塑性本構(gòu)方程；張永敬等[1]結(jié)合Monte-Carlo方法，用Shcapery蠕變型非線性粘彈本構(gòu)模型，分析了推進劑的蠕變行為；彭威等[12]曾經(jīng)從微觀尺度對推進劑的非線性本構(gòu)關(guān)系進行表征。

從應(yīng)變能函數(shù)出發(fā)研究彈性理論，在橡膠、生物膜等材料中已得到廣泛應(yīng)用[13-15]。粘彈性材料由于具有率相關(guān)性、歷史相關(guān)性和應(yīng)力松弛等特性，其應(yīng)變能不僅是應(yīng)變不變量(或主伸長比)的函數(shù)，還應(yīng)強烈依賴于時間。Hrapko等[16]基于運動學(xué)框架導(dǎo)出了非彈性變形的發(fā)展方程，并在Mooney-Rivlin模型的基礎(chǔ)上進行修正，結(jié)合純剪切試驗，建立了腦細胞薄膜的粘彈應(yīng)變能函數(shù)。真實材料或多或少存在蠕變、松弛、遲滯等現(xiàn)象，即表現(xiàn)出粘彈性，Yang等[17]就曾指出，橡膠在高應(yīng)變率變形條件下，同時表現(xiàn)出超彈性和粘彈性。因此，研究粘彈性應(yīng)變能表達式并用于建立本構(gòu)方程，對描述材料的真實特性，發(fā)展非線性粘彈性理論，具有重要的現(xiàn)實意義。HTPB推進劑是典型的粘彈性材料，在許多實際應(yīng)用中，其力學(xué)行為表現(xiàn)出明顯的非線性。因此，研究其變形行為的應(yīng)變能，準確表征應(yīng)力和應(yīng)變之間的關(guān)系，對于裝藥結(jié)構(gòu)完整性的評估，具有很高的應(yīng)用價值，而國內(nèi)外關(guān)于這方面的報道很少。

本文給出了HTPB推進劑粘彈性應(yīng)變能的一般形式，并運用應(yīng)變能函數(shù)表征復(fù)合固體推進劑應(yīng)力和應(yīng)變之間的關(guān)系，避開了復(fù)雜的物理本質(zhì)和發(fā)展過程，建立的本構(gòu)模型能準確描述推進劑材料的力學(xué)行為，且比以往更有利于工程運用。

1 變形的幾何描述

對于任意連續(xù)介質(zhì)，在外力的作用下發(fā)生連續(xù)變形，物質(zhì)內(nèi)任意質(zhì)點從參考構(gòu)形變化到瞬時構(gòu)形，在t時刻的瞬時位置矢量可表示為

x=x(X,t)

(1)

變形梯度張量F定義為

(2)

由式(2)有：

dx=F·dX

(3)

式(3)表明，變形梯度張量F把參考構(gòu)形中的每一個微線元矢量dX映射成瞬時構(gòu)形中的微線元矢量dx。因此，張量F能夠反映一點處變形的全部信息。對于任何一個可逆的張量F，可以有如下的極分解形式：

F=R·U=V·R

(4)

式中R為正交張量，表示純轉(zhuǎn)動；U與V為正定對稱張量，代表純變形。

但在實際運用中，通常用左、右柯西-格林變形張量B和C描述物體變形，B、C分別定義為

B=F·FT，C=FT·F

(5)

以I1、I2和I3分別代表左、右柯西-格林變形張量B、C的第1、第2和第3不變量，則有

I1=tr(C)=λ12+λ22+λ32

(6)

I3=det(C)=λ12λ22λ32

式中 tr和det分別表示張量的跡和行列式；變量λ1、λ2、λ3是3個主伸長比；λi=1+εi，i=1,2,3，εi表示名義應(yīng)變的主應(yīng)變。

為全面描述一點處的變形情況，定義一個量E為

(7)

E稱為Lagrange或Green應(yīng)變張量，求其隨體導(dǎo)數(shù)，得：

(8)

2 模型的建立

2.1 應(yīng)變能函數(shù)

材料變形的應(yīng)變能(或變形能)是在外載荷的作用下發(fā)生變形，而儲存在材料中的能量。研究表明，應(yīng)變能函數(shù)能較好表征物體的變形特征，在超彈性材料中，常用的應(yīng)變能模型有描述橡膠變形行為的Mooney-Rivlin模型、Ogden模型等及描述生物薄膜變形行為的Fung模型、GPR模型等。運用應(yīng)變能函數(shù)描述材料的變形，關(guān)鍵在于應(yīng)變能函數(shù)模型的選擇，不同的應(yīng)變能模型的適用條件和使用范圍是不同的。因此，所選應(yīng)變能函數(shù)應(yīng)盡量準確描述變形行為的真實情況，所建立的應(yīng)變能函數(shù)模型也應(yīng)詳細說明使用條件及范圍；另一方面，函數(shù)中應(yīng)包含盡量少的待求參數(shù)，以減少擬合所需試驗數(shù)據(jù)。

對于各向同性材料，應(yīng)變能函數(shù)可用左(或右)柯西-格林變形張量B(或C)的主不變量I1、I2和I3或變形梯度張量F的特征值，即主伸長比λ1、λ2、λ3寫出：

W=W(I1,I2,I3)或W=W(λ1,λ2,λ3)

(9)

由于復(fù)合固體推進劑為接近不可壓縮各向同性材料，因此假設(shè)推進劑材料不可壓縮，且在變形過程中各向同性仍然有效。對于不可壓縮材料(I3=1)，可得到：

W=W(I1,I2)

(10)

鑒于推進劑單向拉伸的應(yīng)力應(yīng)變曲線與橡膠類材料拉伸響應(yīng)的S型曲線的前半段形狀相同，如圖1所示。本文考慮以橡膠類材料的應(yīng)變能表達式作為函數(shù)模型，描述推進劑的拉伸變形，但由于推進劑具有明顯的粘彈特性，其應(yīng)變能函數(shù)中的材料參數(shù)將不再為恒定值，而是與應(yīng)變率(或時間)等相關(guān)的參數(shù)函數(shù)。

圖1 推進劑與橡膠類材料單向拉伸應(yīng)力應(yīng)變曲線對比

Fig.1 Comparison between uniaxial tensile stress-strain curves of propellant and rubber-like material

Hart-Smith曾在文獻[18]中介紹的用于描述橡膠類材料力學(xué)行為的“Exp-Ln”方程形式為

(11)

式中W為應(yīng)變能；A、a、b為待求材料參數(shù)。

“Exp-Ln”方程能夠較好描述橡膠類材料的變形行為，且具有較少的材料參數(shù)。因此，以式(11)作為推進劑應(yīng)變能的函數(shù)表達式。

函數(shù)中的待求參數(shù)A、a、b僅為試驗數(shù)據(jù)的擬合系數(shù)，不具有任何物理意義。因此，本文所建立應(yīng)變能函數(shù)是對推進劑的力學(xué)行為進行唯象描述，不涉及任何有關(guān)結(jié)構(gòu)的解釋，這對于工程應(yīng)用和宏觀力學(xué)研究是足夠的。

所建立模型應(yīng)具有如下特性：(1)能夠準確描述推進劑的變形響應(yīng)，適用于確定的變形模式；(2)待求參數(shù)函數(shù)的確定需要盡量少的試驗。以上將在試驗部分進行驗證和討論。

2.2 本構(gòu)模型

純力學(xué)變形條件下，考慮加載開始時刻為時間零點，t<0時的變形歷史的影響可忽略不計。所以，推進劑粘彈應(yīng)變能函數(shù)在同一應(yīng)變率下的材料參數(shù)(A、a、b)將為定值，具有與橡膠等的彈性應(yīng)變能函數(shù)相同的特征。因此，認為此時應(yīng)變能函數(shù)與應(yīng)力的關(guān)系可寫成[17]：

(12)

式中σ為Cauchy應(yīng)力張量；p為靜水壓力；I為單位張量。

式(12)為不可壓縮材料的真實應(yīng)力(或Cauchy應(yīng)力)與主不變量(或主伸長比)之間的關(guān)系表達式。其中，真實應(yīng)力σ與變形狀態(tài)有關(guān)，靜水壓力p可取任意值，其對應(yīng)的應(yīng)力分量是不確定的。

3 試驗與分析

采用不同加載速率下的單軸拉伸試驗數(shù)據(jù)擬合材料參數(shù)，試驗在微機控制五頭電子式萬能試驗機上進行，單軸等速拉伸速率分別為0.5、5、20、100、500 mm/min；同時，進行了2、2 000 mm/min下的單軸等速拉伸試驗，用于驗證所建立的數(shù)學(xué)模型。

試驗均在環(huán)境溫度20 ℃、濕度40%條件下進行，試驗材料為標準啞鈴型試件，標距70 mm，截面10 mm×10 mm。以上試驗各進行5組，結(jié)果分別為5組試驗的平均值。

3.1 材料參數(shù)擬合

單軸等速拉伸對試件施加一維應(yīng)力，假設(shè)拉伸方向上主伸長比為λ，則由不可壓條件導(dǎo)出：

(13)

Cauchy-Green變形張量的主不變量：

I1=λ2+2λ-1，I2=2λ+λ-2，I3=1

(14)

試件在拉伸方向上σ11=σ，垂直于拉伸方向上σ22=σ33=0。

將式(5)、式(11)、式(13)代入(12)，得：

σ11=-p+2Aλ2{exp[a(I1-3)]-bln(I1-2)}

(15)

σ22=σ33=-p+2Aλ-1{exp[a(I1-3)]-bln(I1-2)}

(16)

又因為σ22=σ33=0

(17)

σ11為真實應(yīng)力(Cauchy應(yīng)力)，不可壓條件下，名義(工程)應(yīng)力f11可表示為

(18)

用非線性回歸分析方法，擬合拉伸應(yīng)變率為1.19×10-4s-1(0.5 mm/min)、11.9×10-4s-1(5 mm/min)、47.6×10-4s-1(20 mm/min)、0.023 8s-1(100 mm/min)、0.119 s-1(500 mm/min)的單軸拉伸試驗數(shù)據(jù)，擬合曲線與試驗曲線對比見圖2。

(a)100、500 mm/min

(b)0.5、5、20 mm/min

由圖2可看出，擬合結(jié)果與試驗結(jié)果基本吻合，說明本文所提模型能夠較好地表征推進劑的拉伸特性。另一方面，圖2(a)、圖2(b)對比表明，模型對于數(shù)據(jù)的擬合效果與拉伸速率有關(guān)，拉伸速率愈低，擬合效果愈好。因此，所提模型更加適用于描述低應(yīng)變率下推進劑的粘彈性力學(xué)行為。不同應(yīng)變率下模型參數(shù)的取值不同，但呈規(guī)律性分布，參數(shù)擬合結(jié)果見表1。不同應(yīng)變率下模型參數(shù)的散點分布如圖3所示。由圖3可看出，參數(shù)取值與應(yīng)變率之間呈現(xiàn)非線性相關(guān)關(guān)系，可考慮作一元非線性回歸。根據(jù)散點分布趨勢線的走勢，選取以下2種函數(shù)作為理論回歸方程：

表1 不同應(yīng)變率下的模型參數(shù)Table 1 Model parameters at different strain rate

(a)參數(shù)A

(b)參數(shù)a

(c)參數(shù)b 圖3 模型參數(shù)回歸曲線Fig.3 Regression curves of model parameters

同時，對多種不同配方HTPB推進劑的參數(shù)函數(shù)進行了擬合。結(jié)果表明，上述方法具有良好的普遍適用性，多種配方的材料參數(shù)均可表示為關(guān)于應(yīng)變率的冪函數(shù)形式。

(19)

表2 參數(shù)函數(shù)擬合結(jié)果Table 2 Fitting results of parametric functions

代入式(11)，則應(yīng)變能函數(shù)表示為

基于BIM模型，進行空調(diào)系統(tǒng)的自動設(shè)計。主要由三步組成：（1）空調(diào)分區(qū)與負荷計算。由BIM模型通過計算導(dǎo)出輕量化模型，進行自動分區(qū)與負荷計算，輸出房間負荷計算結(jié)果、空調(diào)系統(tǒng)分區(qū)結(jié)果、建筑基本幾何信息；（2）理想系統(tǒng)的生成。輸出理想系統(tǒng)描述文件（包含設(shè)備明細，連接關(guān)系）；（3）實際系統(tǒng)的生成。修正理想系統(tǒng)描述文件，生成實際空調(diào)系統(tǒng)圖。在BIM大力發(fā)展的時代，通過空調(diào)系統(tǒng)的自動設(shè)計，來減輕工程師在設(shè)計中的繁重勞動，能將更多精力投入在方案的選擇上。在實際工程中，減少與建筑、結(jié)構(gòu)設(shè)計的時間滯后性，促使各部門之間及時地互聯(lián)互通。

[1-ln(λ2+2λ-1-2)]-1}

(20)

3.2 模型驗證

為了驗證所建應(yīng)變能模型和本構(gòu)模型的準確性，將應(yīng)變率為4.76×10-4(2 mm/min)、0.476 s-1(2 000 mm/min)的單軸拉伸試驗結(jié)果與本文模型預(yù)測結(jié)果、數(shù)值計算結(jié)果、Mooney-Rivlin模型預(yù)測結(jié)果以及BKZ單積分本構(gòu)模型預(yù)測結(jié)果進行對比，如圖4所示。

圖4表明，所建模型可較好地預(yù)測拉伸速率為2、2 000 mm/min下的名義應(yīng)力，預(yù)測結(jié)果優(yōu)于Mooney-Rivlin模型和數(shù)值仿真，與能有效反映應(yīng)力松弛等粘彈特性的BKZ模型預(yù)測結(jié)果基本相同，且在2 000 mm/min時，一定程度上優(yōu)于BKZ模型，驗證了模型的準確性。由于所提模型避免考慮復(fù)雜的物理本質(zhì)和發(fā)展過程，因此比以往模型更有利于工程應(yīng)用。

(a)2 mm/min

(b)2 000 mm/min

在較大變形時，模型誤差較大的原因是從粘彈性應(yīng)變能函數(shù)出發(fā)，研究推進劑的本構(gòu)關(guān)系，并未考慮在較大變形下的損傷和破壞效應(yīng)，單純運用此模型，無法準確描述推進劑較大變形時的力學(xué)特性。描述推進劑大變形下的力學(xué)行為，需要計及損傷演化和破壞性能的更加普遍的理論模型。

以上同時驗證了所建立應(yīng)變能函數(shù)的有效性，由應(yīng)變能函數(shù)表達式可得應(yīng)變能密度與變形之間的關(guān)系，如圖5所示。從圖5可見，隨變形量增大，推進劑的應(yīng)變能密度呈非線性單調(diào)增長；同一變形條件下，應(yīng)變率越高，推進劑的應(yīng)變能密度越大。

圖5 推進劑應(yīng)變能密度與變形的關(guān)系曲線Fig.5 Relationship between strain energy density and deformation of propellant

4 結(jié)論

(1)所提出的推進劑粘彈性應(yīng)變能函數(shù)和非線性本構(gòu)方程的一般形式，可較好反映材料在有限變形條件下的力學(xué)特性，選取冪函數(shù)作為回歸方程，對不同應(yīng)變率下的模型參數(shù)作一元非線性回歸，得到了材料參數(shù)關(guān)于應(yīng)變率的函數(shù)，由此建立的本構(gòu)方程能夠準確描述拉伸變形條件下的應(yīng)力應(yīng)變關(guān)系，且求解方法和方程形式較為簡單，比以往更有利于工程應(yīng)用。

(2)推進劑在不同拉伸速率下的材料參數(shù)不同，但呈規(guī)律性分布，擬合的參數(shù)函數(shù)反映了材料特性與應(yīng)變率的關(guān)系呈現(xiàn)冪函數(shù)形式；推進劑應(yīng)變能密度隨變形的增大呈非線性單調(diào)增長，且在同一變形條件下，應(yīng)變率越高，推進劑的應(yīng)變能密度越大。

(3)模型對于復(fù)合推進劑應(yīng)力應(yīng)變行為的預(yù)測結(jié)果與試驗結(jié)果吻合較好，但仍然存在不足之處，即較大變形時，模型的誤差較大，且低應(yīng)變率下的預(yù)測效果優(yōu)于較高應(yīng)變率下的預(yù)測效果。因此，所建模型尤其適用于表征低應(yīng)變率拉伸條件下，材料在有限變形內(nèi)的粘彈特性。

[1] 張永敬, 趙光輝, 陽建紅, 等. HTPB復(fù)合固體推進劑非線性本構(gòu)研究[J]. 推進技術(shù), 2000, 21(4): 69-72.

[2] Leaderman H. Elastic and creep properties of filamentous materials and other high polymers[M]. Washington D C: Textile Foundation, 1943.

[3] Ward I M. Mechanical properties of solid polymers[M]. New York: Wiley, 2012.

[4] Shcapery R A. On the characterization of nonlinear viscoelastic materials[J]. Polymer Engineering and Science, 1969, 9(4): 295-310.

[5] Schapery R A. Correspondence principles and a generalized J intergral for large deformation and fracture analysis of viscoelastic media[J]. International Journal of Fracture, 1984, 25(3): 195-223.

[6] Schapery R A. Nonlinear viscoelastic and viscoplastic constitutive equations with growing damage[J]. International Journal of Solids and Structures,1999, 97(1-4): 33-66.

[7] Findley W N, Lai J S, Onaran K. Creep and relaxation of nonlinear viscoelastic materials[M]. Amsterdam: Noth-holland Publishing Company, 1976.

[8] Bernstein B, Kearsley E A, Zapas L J. A study of stress relaxation with finite strain[J]. Trans. Soc. Rheol.,1963, 7(1): 391-410.

[9] Gyoo D J, Sung K Y, Bong K K. A three-dimentional nonlinear viscoelastic constitutive model of solid propellant[J]. International Journal of Solids and Structures, 2000, 37(34): 4715-4732.

[10] 龔建良, 劉佩進, 李強. 基于能量守恒的HTPB推進劑非線性本構(gòu)關(guān)系[J]. 含能材料, 2013, 21(3): 325-329.

[11] 鄧凱, 陽建紅, 陳飛, 等. HTPB復(fù)合固體推進劑本構(gòu)方程[J]. 宇航學(xué)報, 2010, 31(7): 1815-1818.

[12] 彭威, 周建平, 任均國. 復(fù)合固體推進劑非線性粘彈本構(gòu)方程的微觀力學(xué)分析[J]. 固體火箭技術(shù), 1999, 22(4): 23-25.

[13] Khajehsaeid H, Arghavani J, Naghdabadi R. A hyperelastic constitutive model for rubber-like materials[J]. European Journal of Mechanics A/Solids, 2013, 38: 144-151.

[14] Singh F, Katiyar V K, Singh B P. A new strain energy function to characterize apple and potato tissues[J]. Journal of Food Engineering, 2013, 118(2): 178-187.

[15] Garikipati K, Ktepe S G, Miehe C. Elastica-based strain energy functions for soft biological tissue[J]. Journal of the Mechanics and Physics of Solids, 2008, 56(4): 1693-1713.

[16] Hrapko M, Van Dommelen J A W, Peters G W M, et al. The mechanical behaviour of brain tissue: Large strain response and constitutive modelling[J]. Biorheology, 2006, 43(5): 623-636.

[17] Yang L M, Shim V P W, Lim C T. A visco-hyperelastic approah to modelling the constitutive behaviour of rubber[J]. International Journal of Impact Engineering, 2000, 24(6-7): 545-560.

[18] Hart-smith L J. Elasticity parameters for finite deformations of rubber-like materials[J]. The Journal of Applied Mathematics and Physics, 1966, 17(5): 608-626.

(編輯：劉紅利)

Viscoelastic strain energy and nonlinear constitutive model for HTPB composite solid propellant

ZHANG Xiao，ZHENG Jian，PENG Wei，ZHANG Qian-tu

(Ordnance Engineering College，Shijiazhuang 050003，China)

To characterize the mechanical property of HTPB composite solid propellant under finite deformation condition, viscoelastic strain energy and nonlinear constitutive model were studied in this paper. The general forms of strain energy function(SEF) and nonlinear constitutive model were proposed, and the tensile test data were fitted by unitary nonlinear regression analysis. As a consequence, the material parametric function of strain rate was obtained, and on this basis the uniaxial tensile SEF and constitutive equation were established. The nominal stress curves at different strain rate were predicted, furthermore, they were compared with experimental results and prediction results of existing models. The results show the power function relationship between material parameters and strain rate, and the strain energy density proves to rise nonlinearly with strain, the higher the strain rate is, the larger the strain energy density will be at the same deformation. Moreover, the constitutive equation was demonstrated to be accurate on describing the tensile stress-strain relation of propellant, especially on representing the viscoelastic behavior at finite deformation and low strain rate.

HTPB propellant；strain energy；nonlinear regression；parametric function；constitutive equation；deformation

2014-09-01；

：2014-10-22。

張曉(1990—)，男，碩士，研究方向為固體推進劑力學(xué)性能。E-mail:erebuss@outlook.com

V512

A

1006-2793(2015)06-0827-06

10.7673/j.issn.1006-2793.2015.06.014