童 心,許進(jìn)升,李 博

(1.中國船舶科學(xué)研究中心,無錫 214082;2.南京理工大學(xué) 機(jī)械工程學(xué)院,南京 210094)

0 引言

以固體火箭發(fā)動機(jī)為動力的空空導(dǎo)彈在完整的壽命期內(nèi)一般要經(jīng)歷運(yùn)輸、貯存、掛飛、機(jī)動飛行和自主飛行等階段,在長途運(yùn)輸以及掛機(jī)飛行過程中的振動頻率高達(dá)幾千赫茲,作為導(dǎo)彈動力裝置的固體火箭發(fā)動機(jī)也隨之產(chǎn)生高頻振動,發(fā)動機(jī)內(nèi)的復(fù)合推進(jìn)劑藥柱隨著振動載荷的作用會發(fā)生疲勞熱耗散,導(dǎo)致的溫度驟升和熱軟化效應(yīng)對發(fā)動機(jī)的安全性有著顯著影響。

由于粘合劑大分子間的纏結(jié)、相互吸引以及化學(xué)作用,再加上粘合劑分子與填充顆粒間的物理和化學(xué)吸附,復(fù)合推進(jìn)劑內(nèi)部形成了多尺度、多相互作用的多網(wǎng)絡(luò)結(jié)構(gòu),力學(xué)行為可用非線性粘彈性理論來進(jìn)行描述。PARK等[1]建立了含損傷熱粘彈性本構(gòu)模型,利用折算時間將不同溫度下的材料函數(shù)統(tǒng)一至參考溫度下,隨后HA等[2]將模型拓展為三維形式,HINTERHOELZL等[3]完整地提出了相關(guān)的數(shù)值算法,將其應(yīng)用于有限元軟件Abaqus,成功地模擬了端羥基聚丁二烯(HTPB)推進(jìn)劑在單軸和多軸復(fù)雜應(yīng)力條件下的力學(xué)響應(yīng)。張永敬等[4]通過HTPB推進(jìn)劑蠕變回復(fù)實驗來擬合了Schapery非線性本構(gòu)模型的參數(shù),通過與線性粘彈性本構(gòu)計算的結(jié)果和實驗數(shù)據(jù)對比,證明了Schapery非線性本構(gòu)模型能夠很好反映復(fù)合推進(jìn)劑的力學(xué)響應(yīng)。

復(fù)合推進(jìn)劑作為一種溫度敏感的顆粒填充類聚合物,其力學(xué)表現(xiàn)受到變形中能量耗散的影響,能量耗散的宏觀表現(xiàn)是變形過程中自身溫度上升[5]。然而,目前大多數(shù)非線性粘彈性本構(gòu)模型的研究仍集中于單純力學(xué)狀態(tài)的解釋與預(yù)測,對變形過程中的能量耗散及其影響則關(guān)注較少。在高應(yīng)變率、交變加載等動態(tài)條件下,熱交換的特征時間比載荷作用時間大幾個數(shù)量級,由機(jī)械功轉(zhuǎn)化而來的熱能帶來較大的溫升,因而材料因變形而生熱的性質(zhì)(自熱效應(yīng))不容忽視,在這些情況下必須考慮力學(xué)耗散引起的熱源。由于復(fù)合推進(jìn)劑特殊的組成與結(jié)構(gòu),這種自熱現(xiàn)象在高應(yīng)變率沖擊加載和高頻率循環(huán)加載下尤為顯著。在熱力學(xué)框架下,材料的時效變形行為本質(zhì)上是一種非平衡演化過程,在熱力學(xué)理論框架下建立的本構(gòu)模型與守恒定律相一致,具有嚴(yán)格的熱力學(xué)背景和理論基礎(chǔ),不僅滿足熱力學(xué)基本定律,還能表征材料內(nèi)部的不可逆能量耗散過程和材料的熱力學(xué)狀態(tài)、性質(zhì)[6-8],對熱力耦合狀態(tài)有直觀的描述[9]。RODAS等[10]提出了適用于丁苯橡膠有限變形的非線性粘彈性本構(gòu)模型,以不可逆熱力學(xué)內(nèi)變量理論為基礎(chǔ),引入了Helmholtz自由能,并構(gòu)建了內(nèi)變量的演化方程,很好地模擬了丁苯橡膠低周疲勞的生熱現(xiàn)象。KRAIRI等[11]將總應(yīng)變分解為粘彈性應(yīng)變、熱應(yīng)變和粘塑性應(yīng)變,結(jié)合時溫等效原理,在不可逆熱力學(xué)框架下將MILED等[12]的粘彈-粘塑性等溫本構(gòu)模型拓展為非等溫形式,發(fā)展了熱粘彈-粘塑性本構(gòu)模型,并在熱塑性聚合物的有限元仿真中進(jìn)行了完全耦合形式的熱力分析。SHEN等[13]進(jìn)行了類似工作,模型為熱彈-粘塑性本構(gòu)模型。目前,關(guān)于復(fù)合推進(jìn)劑循環(huán)加載下的自熱效應(yīng)未見相關(guān)報道。

為了精確描述復(fù)合推進(jìn)劑循環(huán)加載下的熱力耦合特性,本文基于不可逆熱力學(xué)理論建立了適用于復(fù)合推進(jìn)劑的非線性熱粘彈性本構(gòu)模型,并通過實驗驗證了模型和數(shù)值算法的正確性,本構(gòu)模型成功表征了復(fù)合推進(jìn)劑在不同加載模式下的非線性力學(xué)行為和自熱效應(yīng),為高能固體聚合物的熱力學(xué)表征提供參考。

1 本構(gòu)模型的構(gòu)建

1.1 非線性熱粘彈性本構(gòu)模型

高分子聚合物材料的力學(xué)行為通常具有高度的時間和環(huán)境依賴性。Schapery基于不可逆熱力學(xué)原理,通過引入應(yīng)力或應(yīng)變相關(guān)標(biāo)量,使粘彈性模型非線性化,并滿足熱力學(xué)基本定律。因此,該模型也被稱為熱力學(xué)本構(gòu)模型。在Schapery模型的基礎(chǔ)上,提出以下非線性熱粘彈性本構(gòu)模型：

(1)

(2a)

(2b)

式中aT為時間-溫度等效因子;am為時間-濕度等效因子;aa為時間-老化等效因子。

若濕度等環(huán)境因素保持恒定,可認(rèn)為am=aa=1。本文認(rèn)為復(fù)合推進(jìn)劑屬于各向同性的均質(zhì)材料,滿足連續(xù)介質(zhì)力學(xué)的基本假設(shè),因此有

Cve(t,σ)=2G(t,σ)Idev+3K(t,σ)Ivol

(3)

其中,Idev為四階偏斜張量;Ivol為四階體積張量;G(t,σ)和K(t,σ)分別為剪切模量和體積模量,具有Prony級數(shù)的形式：

(4)

本構(gòu)模型的關(guān)鍵是選擇合適的非線性粘彈性應(yīng)力函數(shù)g1(σ)和g2(σ),它們是構(gòu)成模型的非線性元素。HAJ-ALI等[14]曾采用多項式形式的應(yīng)力函數(shù),該函數(shù)雖考慮了應(yīng)力的影響,但函數(shù)的參數(shù)繁多,難以準(zhǔn)確獲取。此外,MULIANA等[15]也提出了具有指數(shù)形式的非線性函數(shù),但他們所提出的本構(gòu)模型未應(yīng)用于循環(huán)加載實驗。此外,這些模型中函數(shù)參數(shù)對最終結(jié)果的影響不明確。為了克服上述不足,提出具有指數(shù)函數(shù)形式的非線性應(yīng)力函數(shù)：

(5)

其中,k1、m1、k2、m2和C2為待求參數(shù);σ0為線性粘彈性極限應(yīng)力,是區(qū)分線性粘彈性與非線性粘彈性應(yīng)力的臨界值;σvM為von Mises等效應(yīng)力,其表達(dá)式為

(6)

(7)

綜上可知,g1(σ)和g2(σ)是關(guān)于應(yīng)力歷史的函數(shù),也是關(guān)于時間t的函數(shù)。

1.2 產(chǎn)熱方程

材料的循環(huán)屬于不可逆的熱力學(xué)過程,伴隨著能量的耗散,并導(dǎo)致材料自身溫度場發(fā)生變化。為了描述材料的熱力耦合特性,需要建立能量平衡方程。一般情形下,Clausius-Duhem不等式為

(8)

式中q為熱流矢量;ψ為Helmholtz自由能;S為熵;T為溫度;φ為能量耗散率;φmec為由于外部機(jī)械載荷引起的能量耗散率;φT為熱耗散率。

結(jié)合熱力學(xué)第一定律可得到：

(9)

因此,產(chǎn)熱方程為

(10)

式中 ?包括熱彈性源與內(nèi)耦合源;φmec為固有耗散轉(zhuǎn)化為熱能,即自熱效應(yīng)。

自熱效應(yīng)在高速沖擊和高頻循環(huán)載荷下具有重要的意義。若?和φmec均為0,式(10)變?yōu)榉菬崃︸詈系臏囟绕胶夥匠?即溫度場不受材料本構(gòu)關(guān)系的影響。根據(jù)Fourier定律,熱流矢量表示為

q=-κ·T

(11)

式中κ為導(dǎo)熱張量,若導(dǎo)熱是各向同性的,則κ=κl,κ為熱導(dǎo)率。

在無外部熱源的情況下,rext=0,式(3)變?yōu)?/p>

(12)

式(12)表明,在無外部熱源作用時,材料在循環(huán)過程中的溫度變化主要受熱彈性響應(yīng)、非彈性響應(yīng)和熱傳遞(熱傳導(dǎo)、與環(huán)境的熱對流及熱輻射)三方面的影響。與非彈性耗散源相比,熱彈性效應(yīng)對試件表面溫度變化的影響很有限,僅僅引起溫度的微小波動;對處于粘彈性變形的復(fù)合推進(jìn)劑而言,內(nèi)部結(jié)構(gòu)改變不涉及相變,因而內(nèi)耦合源也可以忽略[16]。

綜上所述,在關(guān)于復(fù)合推進(jìn)劑能量耗散的分析中,可認(rèn)為?≈0。

1.3 能量耗散

線性熱粘彈性本構(gòu)模型能量耗散φtve的表達(dá)式為

φmec=φtve

(13)

其中,φ和m均為松弛函數(shù),反映了溫度因素的影響。若溫度的變化幅度不大,材料力學(xué)特性受溫度變化的影響較小,可假設(shè)φ和m均為0,則有

(14)

該式與等溫條件下線性粘彈性本構(gòu)模型推導(dǎo)的能量耗散表達(dá)式相似,因而：

(15)

式中si和σHj為應(yīng)力的粘性部分。

LéVESQUE等[17]已經(jīng)證明,建立本構(gòu)模型時在Schapery模型中引入不同的非線性函數(shù)依然滿足熱力學(xué)基本定律,且不對Clausius-Duhem不等式產(chǎn)生影響。因此,非線性粘彈性本構(gòu)模型中的能量耗散采用了和線性模型一致的計算公式。

2 本構(gòu)模型的二次開發(fā)

為模擬聚合物循環(huán)過程中的熱力耦合,有限元仿真主要采用了下列兩種方法：

(1)將材料受載后的自熱視為內(nèi)部熱源,通過有限元軟件子程序的拓展,先進(jìn)行力學(xué)分析,然后將計算而來的內(nèi)部熱源應(yīng)用到溫度場的求解中[18-19]。RIAHI等[20]、SONG等[21]和LUO等[22]針對不同聚合物材料進(jìn)行了相關(guān)研究。

(2)建立準(zhǔn)確可靠的本構(gòu)模型,通過熱力學(xué)基本定律,可以從本構(gòu)模型中獲取能量耗散的表達(dá)式,從而在有限元分析中同時求解熱學(xué)和力學(xué)問題[23]。依此方法,針對具體研究對象選取適當(dāng)?shù)淖杂赡芎瘮?shù)形式,引入適當(dāng)?shù)膬?nèi)變量來描述不可逆變形過程中的材料內(nèi)部微結(jié)構(gòu)的變化,即可獲得材料的本構(gòu)關(guān)系,然后就可以從能量平衡方程中推導(dǎo)出熱力學(xué)變形過程中的溫度控制方程[24]。

為了實現(xiàn)完全耦合形式的熱力分析,本文選擇第二種方法。

2.1 數(shù)值算法

應(yīng)力表達(dá)式為σ(t)=s(t)+σH(t)l,其中：

(16)

其中,s∞(t)和σH∞(t)為彈性部分;si(t)和σHj(t)為粘性部分,它們的表達(dá)式如下：

(17)

根據(jù)gi(σ)的表達(dá)式,gi(σ)也是關(guān)于t的函數(shù),為了簡化,記gi(t)≡gi(σ(t))。給定時間增量[tn,tn+1],則tn+1時的應(yīng)力偏量s∞(t)為

(18)

si(tn+1)為

(19)

為了獲得應(yīng)力更新方程,令

(20)

由于t=tn+1時應(yīng)力是未知的,因而gi(tn+1)也是未知的。若時間增量步Δt很小,可以忽略非線性函數(shù)在這一時間區(qū)間內(nèi)的變化,作出以下假設(shè)：g1(tn+1)≈g1(tn)、g2(tn+1)≈g2(tn)。應(yīng)力偏量部分可簡化為

(21)

其中

(22)

同樣地,對于應(yīng)力的球量部分,應(yīng)力更新方程為

(23)

(24)

其中

(25)

綜合可得

(26)

其中

(27)

Jacobian矩陣定義為

(28)

2.2 非線性模型子程序

Abaqus擁有用戶材料子程序,可方便地對材料本構(gòu)模型進(jìn)行二次開發(fā),在高聚物力學(xué)特性的數(shù)值模擬中得到了廣泛的應(yīng)用。本文選擇UMATHT與UMAT聯(lián)合使用來實現(xiàn)完全耦合的熱力分析。在時間增量步中,UMAT通過狀態(tài)變量“rpl”(rpl的值即為式(12)中的φmec)存儲積分點(diǎn)處由于機(jī)械加載引起的能量耗散。UMATHT提供了熱流矢量,定義了材料的熱本構(gòu)關(guān)系,根據(jù)產(chǎn)熱方程求解熱傳遞問題。圖1是執(zhí)行子程序的流程圖。

圖1 子程序UMAT和UMATHT流程圖Fig.1 Flowchart of UMAT and UMATHT

由于復(fù)合推進(jìn)劑是對溫度敏感的材料,應(yīng)力場和溫度場相互影響,因而本文采取了完全耦合的熱力分析方法。在完全的熱力耦合分析中,耦合方程的矩陣表達(dá)式為

(29)

式中 Δu和Δθ分別為位移增量和溫度增量;Kij為分矩陣;Ru和Rθ分別為力學(xué)和熱學(xué)殘余向量。

3 仿真結(jié)果與分析

模型參數(shù)的獲取方法見文獻(xiàn)[25]。利用HTPB推進(jìn)劑單軸拉伸實驗和應(yīng)變控制的單軸循環(huán)實驗對本構(gòu)模型的仿真結(jié)果進(jìn)行驗證,試驗材料、設(shè)備和試驗流程的具體說明可參考文獻(xiàn)[26],在此不再贅述。仿真時物理模型的幾何形狀與尺寸與真實啞鈴型試件一致,邊界條件與實驗一致,試件的一端固定,另一端施加位移載荷,見圖2,選擇試件標(biāo)距段中部作為分析區(qū)域。環(huán)境溫度根據(jù)實驗要求設(shè)定。試件與環(huán)境之間為自然對流換熱,試件表面的熱邊界條件為

圖2 Abaqus/CAE中試件邊界條件示意圖Fig.2 Schematic diagram of boundary conditions of the specimen in Abaqus/CAE

q=-h(Ts-T∞)n

(30)

式中h為試件與環(huán)境之間的對流換熱系數(shù);Ts為試件表面溫度;T∞為環(huán)境溫度;n為表面的法線方向。

3.1 單軸加載

聯(lián)合使用子程序UMAT和UMATHT,以模擬復(fù)合推進(jìn)劑由于機(jī)械載荷引起的自熱效應(yīng),探討其對復(fù)合推進(jìn)劑力學(xué)響應(yīng)的影響。在自熱效應(yīng)的有限元仿真中,UMATHT中的導(dǎo)熱系數(shù)和比熱容是常數(shù),不隨溫度變化而變化。HTPB推進(jìn)劑試件與空氣之間的對流換熱系數(shù)h設(shè)置為10 W/(m2·K)。

圖3 單軸拉伸實驗中的自熱效應(yīng)Fig.3 Self-heating effect in uniaxial tensile test

圖4顯示了自熱效應(yīng)對于單軸應(yīng)力-應(yīng)變關(guān)系的影響,未考慮自熱效應(yīng)的預(yù)測曲線是在等溫條件下獲取的。由于能量耗散導(dǎo)致了溫升,HTPB推進(jìn)劑受熱而軟化,因此預(yù)測的應(yīng)力整體偏小。

圖4 自熱效應(yīng)對單軸應(yīng)力-應(yīng)變關(guān)系的影響Fig.4 Effect of self-heating on uniaxial tensile stress-strain relations

為反映自熱效應(yīng)與應(yīng)變率之間的關(guān)系,還進(jìn)行了常溫下其他應(yīng)變率下的仿真,見圖5。隨著應(yīng)變率的提高,外載荷做功越劇烈,能量耗散越來越多,自熱效應(yīng)也越來越顯著。該發(fā)現(xiàn)揭示了高聚物材料變形中自熱效應(yīng)的應(yīng)變率相關(guān)性,為了更精確地描述材料的力學(xué)性能,建立本構(gòu)模型時需根據(jù)變形與加載速率等條件考慮材料自身溫度的改變。對復(fù)合推進(jìn)劑而言,在高應(yīng)變率加載的大應(yīng)變情形下,自熱效應(yīng)較顯著[27]。

圖5 不同應(yīng)變率下的自熱效應(yīng)Fig.5 Self-heating effect under various strain rates

3.2 循環(huán)加載

為了明確自熱效應(yīng)對循環(huán)行為的影響,隨后對受循環(huán)加載的HTPB推進(jìn)劑試件進(jìn)行了完全熱力耦合形式的有限元仿真。圖6是試件加載過程中典型的溫度云圖(加載條件為f=50 Hz,εa=0.03)?？梢?試件的標(biāo)距段的溫度由于材料力學(xué)耗散而升高,軸向(與試件加載平行的方向)上存在溫度梯度,這與實驗觀測一致。圖6中還給出了中心區(qū)域截面的溫度云圖,截面核心溫度與邊緣溫度相差極小,因而可認(rèn)為試件內(nèi)部實現(xiàn)了溫度平衡,各空間點(diǎn)的溫度趨向一致,溫度場在這一截面是均勻分布的,這也進(jìn)一步證明了此前實驗分析中利用表面溫度表征HTPB推進(jìn)劑力學(xué)性能的合理性。

圖6 循環(huán)加載中試件典型的溫度云圖Fig.6 Typical temperature contour of the specimen in cyclic loading

采用對稱建模可以減少工作量,既不影響計算結(jié)果又可提升計算效率。由于需要同時使用兩個子程序,考慮試件的對稱性,在后續(xù)的仿真中選取試件的一半進(jìn)行分析。從文獻(xiàn)[28]可知,由于HTPB推進(jìn)劑力學(xué)性質(zhì)的溫度相關(guān)性,自熱效應(yīng)引起的溫升影響著HTPB推進(jìn)劑動態(tài)應(yīng)力-應(yīng)變關(guān)系的預(yù)測。本文進(jìn)行了多種加載條件下的自熱效應(yīng)模擬,得到了試件標(biāo)距段表面溫度隨時間的變化,結(jié)果見圖7。

(a)Strain amplitude 0.01～0.03

由圖7可見,由于開始時試件和環(huán)境的溫差較小,對流過程中熱量損失也比較少,試件的產(chǎn)熱速率大于試件和環(huán)境的熱交換率,因而溫度在一段時間逐漸上升;隨著循環(huán)時間的增加,試件表面溫度與環(huán)境溫度通過熱對流而達(dá)到平衡,表面溫度保持相對穩(wěn)定的狀態(tài)。綜上,仿真預(yù)測結(jié)果與實驗中測得的溫升較吻合,驗證了本構(gòu)模型和數(shù)值算法的正確性。

圖8給出了不同仿真條件下循環(huán)應(yīng)力-應(yīng)變關(guān)系的對比圖?？梢园l(fā)現(xiàn),與不考慮自熱效應(yīng)相比,滯回圈的斜率降低,表明力學(xué)性能已受到自身溫度升高的影響。在熱力耦合分析中,考慮自熱效應(yīng)后,滯回圈的形狀與實驗結(jié)果更為吻合,非線性熱粘彈性本構(gòu)模型的預(yù)測精度得到提高,表明計及自熱效應(yīng)的熱力耦合型本構(gòu)模型能更準(zhǔn)確地描述HTPB推進(jìn)劑的循環(huán)力學(xué)行為。

為了定量評估模型的預(yù)測效果,計算了滯回圈斜率的誤差,見圖9。由圖9可知,模型的預(yù)測與實驗之間的誤差均在10%以內(nèi),表明與文獻(xiàn)[28]相比,本構(gòu)模型的預(yù)測能力得到了提高,有力佐證了自熱效應(yīng)在本構(gòu)模型應(yīng)用中的重要性。

(a)f=10 Hz,εa=0.01 (b)f=1 Hz,εa=0.03 (c)f=10 Hz,εa=0.03

圖9 滯回圈斜率的預(yù)測誤差Fig.9 Prediction errors of the slope of the hysteresis loop

4 結(jié)論

(1)對于循環(huán)加載下的復(fù)合推進(jìn)劑,自熱效應(yīng)較單軸加載狀態(tài)更顯著(最大溫升為15 K),需要在本構(gòu)模型中加以考慮。

(2)自熱效應(yīng)對復(fù)合推進(jìn)劑的力學(xué)行為有著直接影響,與不考慮自熱效應(yīng)相比,循環(huán)加載下應(yīng)力-應(yīng)變滯回圈的斜率降低,其力學(xué)性能的劣化受到自熱效應(yīng)的直接影響。

(3)本構(gòu)模型成功表征了復(fù)合推進(jìn)劑在不同加載模式下的非線性力學(xué)行為和自熱效應(yīng),模型預(yù)測的滯回圈斜率與實驗之間的誤差均在10%以內(nèi)。該研究為高能固體聚合物的熱力學(xué)表征提供了參考。