HTPB推進(jìn)劑蠕變本構(gòu)方程研究*

羅一智，沙寶林，侯 曉

(中國航天科技集團(tuán)有限公司四院四十一所,西安 710025)

符號(hào)說明

σ——應(yīng)力，MPa

εcr——蠕變應(yīng)變

t——時(shí)間，s

S-H——應(yīng)變硬化方程縮寫

T-H——時(shí)間硬化方程縮寫

M-T-H——改進(jìn)時(shí)間硬化方程縮寫

M-S-H——改進(jìn)應(yīng)變硬化方程縮寫

g1、g2…g4——表示函數(shù)關(guān)系，無具體形式

G1、G2、G3——表示函數(shù)關(guān)系，無具體形式

a、b、c、d——為簡(jiǎn)化方程形式引入的中間量

e、f、r、p——為簡(jiǎn)化方程形式引入的中間量

A、B、C——為簡(jiǎn)化方程形式引入的中間量

0 引言

固體火箭發(fā)動(dòng)機(jī)經(jīng)歷長(zhǎng)期的立式貯存工況，承受溫度載荷、自重載荷，藥柱產(chǎn)生蠕變變形發(fā)生下沉，甚至阻塞排氣通道，這些因素都對(duì)發(fā)動(dòng)機(jī)的結(jié)構(gòu)可靠性和壽命產(chǎn)生影響[1]。對(duì)于藥柱長(zhǎng)期立式貯存的特殊工況及蠕變效應(yīng)，近些年愈來愈得到重視。國防科技大學(xué)的沈懷榮[2]，研究了復(fù)合固體推進(jìn)劑的溫度相關(guān)的蠕變損傷模型；針對(duì)大力神-4固體助推器的立式貯存問題[3]，國外研究人員利用Rocstar程序包進(jìn)行模擬，研究大變形情況；南京理工大學(xué)的李東等[4-5]，針對(duì)雙基固體推進(jìn)劑，基于冪律蠕變模型，結(jié)合試驗(yàn)數(shù)據(jù)，研究藥柱本構(gòu)關(guān)系和蠕變特性；袁軍等[6]針對(duì)大型固體發(fā)動(dòng)機(jī)，考慮溫度載荷和內(nèi)壓載荷的作用，進(jìn)行燃燒室有限元分析；王永帥等[7]基于時(shí)間硬化蠕變方程，研究艦載導(dǎo)彈發(fā)動(dòng)機(jī)的蠕變效應(yīng)；Bihari等[8]采用經(jīng)典的Kelvin-Voigt模型，對(duì)推進(jìn)劑的粘彈特性進(jìn)行研究，建立了蠕變粘彈的數(shù)學(xué)模型；海軍航空大學(xué)的王鑫等[9-11]，對(duì)于HTPB推進(jìn)劑，基于時(shí)間硬化模型，研究了考慮蠕變效應(yīng)的立式貯存發(fā)動(dòng)機(jī)在多種載荷作用下的應(yīng)力應(yīng)變狀態(tài)。然而，至今有關(guān)固體發(fā)動(dòng)機(jī)的蠕變效應(yīng)研究，多采用時(shí)間硬化蠕變方程。前人關(guān)于蠕變方程的研究則多集中于金屬和巖石領(lǐng)域，對(duì)于固體推進(jìn)劑的蠕變本構(gòu)方程及其適用性的全面研究尚未見報(bào)道。

本文針對(duì)某配方HTPB推進(jìn)劑發(fā)動(dòng)機(jī)的立式貯存問題，開展了多應(yīng)力水平的蠕變?cè)囼?yàn)。根據(jù)試驗(yàn)數(shù)據(jù)，擬合了多種蠕變本構(gòu)方程，并對(duì)其適用性進(jìn)行了研究。

1 試驗(yàn)

1.1 試驗(yàn)儀器與方案

試驗(yàn)儀器為CMT4203-3A蠕變應(yīng)力-松弛儀。該儀器的最大試驗(yàn)力為2 kN，準(zhǔn)確度等級(jí)為0.5級(jí)，試驗(yàn)力測(cè)量范圍為10%～100%，試驗(yàn)力示值分辨力為最大試驗(yàn)力的1/300 000，試驗(yàn)力示值相對(duì)誤差為示值的±0.5%以內(nèi)，位移示值相對(duì)誤差±0.5%，位移分辨為0.03 μm。試驗(yàn)箱使用溫度為-30～100 ℃，溫度波動(dòng)度為±2 ℃。

試驗(yàn)對(duì)象是某配方HTPB推進(jìn)劑試件，試件依照QJ 924—1985《復(fù)合固體推進(jìn)劑單向拉伸試驗(yàn)方法》規(guī)定執(zhí)行，推進(jìn)劑配方如表1所示。

表1 典型HTPB推進(jìn)劑配方

蠕變?cè)囼?yàn)在常溫下進(jìn)行，考慮到HTPB試件最大抗拉強(qiáng)度在1 MPa左右、藥柱貯存過程中的界面最大應(yīng)力以及試驗(yàn)過程所需大量時(shí)間的限制，分別進(jìn)行0.4、0.5、0.6、0.7 MPa應(yīng)力水平的蠕變?cè)囼?yàn)，每個(gè)應(yīng)力水平開展3組試驗(yàn)，通過計(jì)算機(jī)程序?qū)崟r(shí)記錄應(yīng)力和試件的伸長(zhǎng)量。根據(jù)標(biāo)距長(zhǎng)度計(jì)算蠕變應(yīng)變，對(duì)各應(yīng)力水平的3組試驗(yàn)取平均值，得到試件的蠕變曲線。蠕變?cè)囼?yàn)裝置如圖1所示。

圖1 蠕變?cè)囼?yàn)裝置Fig.1 Set-up for creep test

1.2 試驗(yàn)結(jié)果

圖2為不同應(yīng)力水平下的應(yīng)變-時(shí)間曲線。試件受到載荷作用時(shí)，首先產(chǎn)生瞬時(shí)彈性應(yīng)變；定載荷的持續(xù)作用下，試件應(yīng)變不斷增大。初始階段，應(yīng)變率隨時(shí)間的增長(zhǎng)而減小；第二階段也稱為穩(wěn)態(tài)階段，應(yīng)變率趨于恒定，應(yīng)變穩(wěn)定增大；第三階段，應(yīng)變率變大，應(yīng)變快速增大，試件破壞。應(yīng)力越大，試件破壞時(shí)間越短。

(a)σ=0.4 MPa (b)σ=0.5 MPa (c)σ=0.6 MPa (d)σ=0.7 MPa

2 蠕變本構(gòu)方程

固體推進(jìn)劑屬于高聚物材料，具有粘彈性特性，在長(zhǎng)期自重載荷的作用下，產(chǎn)生蠕變效應(yīng)。一般情況，蠕變應(yīng)變率是時(shí)間、應(yīng)力、應(yīng)變、溫度的函數(shù)，如式(1)：

(1)

對(duì)方程進(jìn)行積分、移項(xiàng)處理，將方程化為

εcr=G1(σ)G2(t)G3(T)

(2)

本文研究的蠕變本構(gòu)方程[12-13]匯總見表2。

表2 蠕變本構(gòu)方程

蠕變本構(gòu)方程包括應(yīng)力相關(guān)項(xiàng)、時(shí)間相關(guān)項(xiàng)和溫度相關(guān)項(xiàng)。本文研究對(duì)象基于保溫筒中發(fā)動(dòng)機(jī)燃燒室，溫度保持恒定。常溫下，不同組蠕變?cè)囼?yàn)的溫度保持一致，可忽略方程中的溫度相關(guān)項(xiàng)，對(duì)于溫度相關(guān)項(xiàng)的系數(shù)，進(jìn)行歸零處理；對(duì)于某個(gè)應(yīng)力水平下的蠕變?cè)囼?yàn)，應(yīng)力為定值，G1(σ)為常數(shù)，通過G2(t)項(xiàng)描述某定應(yīng)力水平下的蠕變過程；對(duì)于不同應(yīng)力水平下的蠕變過程，通過G1(σ)項(xiàng)描述規(guī)律性。

針對(duì)蠕變本構(gòu)方程研究的總體思路(見圖3)：

圖3 總體研究思路Fig.3 General research process

(1)選擇蠕變本構(gòu)方程，采用最小二乘法針對(duì)不同應(yīng)力的蠕變?cè)囼?yàn)進(jìn)行數(shù)據(jù)擬合，確定不同應(yīng)力水平下的G1(σ)和G2(t)。對(duì)于模型參數(shù)是非線性的函數(shù)，采用Levenberg-Marquardt方法[14]迭代處理。Levenberg-Marquardt方法是非線性最小二乘估計(jì)的一種估計(jì)方法，在其中引入了阻尼因子，結(jié)合了高斯-牛頓法和梯度下降法的優(yōu)勢(shì)。

(2)研究分析不同應(yīng)力水平下G1(σ)與G2(t)的規(guī)律性和一致性，判斷其與本構(gòu)方程的匹配度，確定適用于所有應(yīng)力水平的方程參數(shù)。

(3)比較分析各方程與試驗(yàn)數(shù)據(jù)的擬合程度，最終確定適用于推進(jìn)劑的蠕變本構(gòu)方程。

3 蠕變本構(gòu)方程的適用性

蠕變是時(shí)間相關(guān)的非線性過程，蠕變本構(gòu)方程形式多樣，數(shù)學(xué)特征具有異同性。針對(duì)不同方程，采取不同的策略，研究方程的推進(jìn)劑適用性。

3.1 冪律類型蠕變本構(gòu)方程

此類蠕變本構(gòu)方程的特點(diǎn)是，通過對(duì)蠕變方程積分、移項(xiàng)、合并等處理方式，可將本構(gòu)方程化為εcr=atb形式，且b為常數(shù)，a為σ的冪函數(shù)形式。

(1)應(yīng)變硬化方程

(3)

對(duì)方程進(jìn)行積分、移項(xiàng)、合并等處理，方程可化為

(4)

(2)時(shí)間硬化方程

(5)

對(duì)方程進(jìn)行積分、移項(xiàng)、合并等處理，方程可化為

(6)

(3)改進(jìn)時(shí)間硬化方程

(7)

(4)改進(jìn)應(yīng)變硬化方程

(8)

對(duì)方程進(jìn)行積分、移項(xiàng)、合并等處理，方程可化為

εcr=C1σC2(C3+1)C3tC3+1

(9)

于是，a=C1(C3+1)C3σC2,b=C3+1。

確定各方程a、b的具體形式，再針對(duì)各應(yīng)力水平的蠕變?cè)囼?yàn)進(jìn)行數(shù)據(jù)擬合，確定適用于各自應(yīng)力水平的a、b值，結(jié)果如表3所示。

為統(tǒng)一描述各應(yīng)力水平的蠕變過程，結(jié)合各組試驗(yàn)數(shù)據(jù)擬合a與σ的函數(shù)關(guān)系，結(jié)果為

a=0.02042×σ2.63112，R2=0.94306

其中，R2為決定系數(shù)，表征回歸分析的擬合程度，其值越接近1,則模型擬合度越高。R2接近1，表明a與σ的函數(shù)關(guān)系滿足蠕變本構(gòu)方程要求，冪律類型蠕變本構(gòu)方程適用于此推進(jìn)劑蠕變過程。

上述四種方程最終參數(shù)結(jié)果如表4所示。最終擬合曲線與試驗(yàn)數(shù)據(jù)對(duì)比如圖4所示，擬合值與試驗(yàn)值的殘差平方和RSS如表5所示。

表4 冪律類方程各參數(shù)值

(a)σ=0.4 MPa (b)σ=0.5 MPa

表5 冪律類方程擬合值與試驗(yàn)值殘差平方和

需要指出的是，冪律類各蠕變方程在物理意義與方程形式上存在明顯差異，但由于各方程均能化為時(shí)間的冪函數(shù)形式，最終擬合曲線具有相同結(jié)果。結(jié)果表明，冪律類蠕變方程擬合曲線與試驗(yàn)曲線雖有一定偏差，但趨勢(shì)一致性較好，能夠反映蠕變過程。

3.2 廣義指數(shù)蠕變本構(gòu)方程

對(duì)原方程進(jìn)行積分、移項(xiàng)、合并等處理方式，可將方程化為

εcr=-C1σC2e-rt+B

r=C5σC3e-C4/T

(10)

對(duì)于某組應(yīng)力水平的蠕變?cè)囼?yàn)，應(yīng)力相關(guān)項(xiàng)為定值，上式等價(jià)于

εcr=-Ae-rt+B

r=C5σC3e-C4/T

(11)

式中A=C1σC2。

對(duì)各應(yīng)力水平的試驗(yàn)數(shù)據(jù)進(jìn)行擬合，確定A、B的值，如圖5所示。根據(jù)圖5，擬合計(jì)算A=C1σC2中C1、C2的值，擬合結(jié)果為A=0.11396×σ0.62424，R2=0.2821。

(a)Values of A

關(guān)于A=C1σC2的擬合計(jì)算中，R2遠(yuǎn)小于1，故A無法滿足A=C1σC2函數(shù)條件，廣義指數(shù)蠕變本構(gòu)方程不適用于此推進(jìn)劑的蠕變過程。

3.3 廣義格雷厄姆蠕變本構(gòu)方程

對(duì)原方程進(jìn)行積分，并取C8=0，可化為

(12)

于是，可簡(jiǎn)寫為

εcr=atb+ctd+etf

(13)

其中

根據(jù)已有經(jīng)驗(yàn)和相關(guān)文獻(xiàn)[15]，取b=0.2，d=1，f=2，針對(duì)各應(yīng)力水平的蠕變?cè)囼?yàn)數(shù)據(jù)進(jìn)行擬合，確定a、c、e的值，如表6所示。

表6 廣義格雷厄姆方程的a、c、e的值

根據(jù)a、c、e與應(yīng)力的冪函數(shù)關(guān)系，對(duì)表中數(shù)據(jù)進(jìn)行冪函數(shù)擬合，結(jié)果如下：

a=1.1149×10-2×σ0.45845,R2=0.56025

c=2.2045×10-3×σ0.45845,R2=0.13102

e=-3.5289×10-9×σ0.45845,R2=0.11934

三個(gè)中間量的關(guān)于冪函數(shù)擬合計(jì)算的決定系數(shù)R2均遠(yuǎn)小于1，故a、c、e不具備應(yīng)力的冪函數(shù)關(guān)系，說明此廣義格雷厄姆方程不適用于推進(jìn)劑的蠕變過程。

3.4 穩(wěn)定階段蠕變本構(gòu)方程

蠕變過程根據(jù)蠕變應(yīng)變率的變化，劃分為三個(gè)階段。蠕變第二階段蠕變率相對(duì)恒定，稱為穩(wěn)定階段。穩(wěn)定階段蠕變本構(gòu)方程只針對(duì)于蠕變第二階段，描述蠕變第二階段蠕變率與應(yīng)力水平的關(guān)系，即蠕變-時(shí)間曲線中穩(wěn)定階段的斜率與應(yīng)力的關(guān)系。

(1)廣義伽羅伐洛方程

(14)

(2)指數(shù)方程

(15)

(3)諾頓方程

(16)

式中 溫度相關(guān)項(xiàng)的參數(shù)取為0，則式(14)式～(16)均為關(guān)于應(yīng)力的函數(shù)。

針對(duì)各應(yīng)力水平的蠕變?cè)囼?yàn)，確定其穩(wěn)定階段的斜率，再根據(jù)斜率對(duì)以上三方程進(jìn)行數(shù)據(jù)擬合，擬合結(jié)果如圖6所示。結(jié)果表明，三種方程均能描述蠕變穩(wěn)定階段的斜率變化趨勢(shì)。

圖6 第二階段斜率的擬合曲線Fig.6 Fitting curves of secondary period slopes

圖7為三種方程擬合結(jié)果與蠕變?cè)囼?yàn)數(shù)據(jù)的對(duì)比結(jié)果。三種方程均能描述蠕變穩(wěn)定階段，但不具備描述蠕變初始階段的能力。本文研究的是描述蠕變前兩階段的蠕變方程，故此類方程不適用于推進(jìn)劑的蠕變過程。但在具體應(yīng)用過程中，若確定推進(jìn)劑的蠕變過程處在第二階段，可考慮采用以上三種方程。

(a)σ=0.4 MPa (b)σ=0.5 MPa

3.5 復(fù)合時(shí)間硬化蠕變本構(gòu)方程

(17)

復(fù)合時(shí)間硬化蠕變本構(gòu)方程無需進(jìn)行積分處理，溫度相關(guān)項(xiàng)參數(shù)取為0，即C4=0、C7=0。對(duì)于某應(yīng)力水平的蠕變?cè)囼?yàn)而言，式(17)可寫為

εcr=AtB+Ct

(18)

針對(duì)各應(yīng)力水平擬合確定A、C的值，再確定A、C與應(yīng)力的冪函數(shù)關(guān)系，結(jié)果如下：

A=0.01351×σ1.64377,R2=0.99383

C=0.0172×σ18.82979,R2=0.99645

其中,R2均接近1，表明A、C與σ的函數(shù)關(guān)系滿足蠕變本構(gòu)方程要求，復(fù)合時(shí)間硬化蠕變本構(gòu)方程適用于此推進(jìn)劑蠕變過程。最終結(jié)果如表7及圖8所示。

表7 復(fù)合時(shí)間硬化蠕變本構(gòu)方程系數(shù)值

(a)σ=0.4 MPa (b)σ=0.5 MPa

最終擬合曲線與試驗(yàn)數(shù)據(jù)對(duì)比見圖8，擬合值與試驗(yàn)值的殘差平方和RSS如表8所示。結(jié)果表明，復(fù)合時(shí)間硬化蠕變本構(gòu)方程擬值線與試驗(yàn)數(shù)據(jù)偏差最小，擬合曲線與試驗(yàn)曲線一致性最好，復(fù)合時(shí)間硬化蠕變本構(gòu)方程能夠很好地描述推進(jìn)劑的蠕變過程。

表8 復(fù)合時(shí)間硬化蠕變本構(gòu)方程擬合值與試驗(yàn)值殘差平方和

3.6 有理多項(xiàng)式蠕變本構(gòu)方程

有理多項(xiàng)式蠕變本構(gòu)方程具有極其復(fù)雜的形式，本文對(duì)其簡(jiǎn)化形式進(jìn)行研究。對(duì)方程積分簡(jiǎn)化后：

(19)

根據(jù)各應(yīng)力水平的蠕變?cè)囼?yàn)數(shù)據(jù)確定C、b、p值，再研究分析C、b、p與應(yīng)力的關(guān)系。C、b、p結(jié)果如圖9所示。

(a)Values of C

基于上述數(shù)據(jù)，擬合確定C、b、p與應(yīng)力的冪函數(shù)關(guān)系：

C=5.887×10-2×σ0.3606,R2=0.6092

b=8.0202×10-4×σ9.2175,R2=0.9938

p=1.99×10-2×σ2.5639,R2=0.9788

結(jié)合圖形與計(jì)算結(jié)果表明，關(guān)于C的冪函數(shù)擬合計(jì)算的決定系數(shù)R2<<1，不具備此有理多項(xiàng)式蠕變本構(gòu)方程要求的冪函數(shù)形式，此方程不適用于推進(jìn)劑的蠕變過程。

3.7 廣義時(shí)間硬化蠕變本構(gòu)方程

對(duì)于某應(yīng)力水平的蠕變?cè)囼?yàn)以及常溫條件，原方程可簡(jiǎn)化為

εcr=ftr

(20)

式中f=C1σ+C2σ2+C3σ3；r=C4+C5σ。

針對(duì)各應(yīng)力水平，擬合確定f、r的值，進(jìn)一步確定f、r與應(yīng)力的關(guān)系，結(jié)果如下：

f=0.02756×σ-0.07577×σ2+0.07208×σ3

R2=0.92548

r=0.19786+0.2221×σ,R2=1

其中，R2均接近1，表明f、r與σ的函數(shù)關(guān)系滿足蠕變本構(gòu)方程要求，廣義時(shí)間硬化蠕變本構(gòu)方程適用于此推進(jìn)劑蠕變過程。最終確定所有參數(shù)的值，結(jié)果如表9所示。

表9 廣義時(shí)間硬化蠕變本構(gòu)方程系數(shù)值

最終擬合曲線與試驗(yàn)數(shù)據(jù)對(duì)比如圖10所示，擬合值與試驗(yàn)值的殘差平方和RSS如表10所示。結(jié)果表明，廣義時(shí)間硬化蠕變本構(gòu)方程擬合值與試驗(yàn)數(shù)據(jù)偏差大于冪律類方程與復(fù)合時(shí)間硬化方程，但在0.4、0.6、0.7 MPa下的偏差小于冪律類方程，擬合曲線與試驗(yàn)曲線一致性較好，廣義時(shí)間硬化蠕變本構(gòu)方程能夠較好地描述推進(jìn)劑的蠕變過程。

(a)σ=0.4 MPa (b)σ=0.5 MPa

表10 廣義時(shí)間硬化蠕變方程擬合值與試驗(yàn)值殘差平方和

4 結(jié)論

本文針對(duì)某配方HTPB推進(jìn)劑發(fā)動(dòng)機(jī)的立式貯存問題，開展了多應(yīng)力水平的蠕變?cè)囼?yàn)。根據(jù)試驗(yàn)數(shù)據(jù)，研究擬合了多種蠕變本構(gòu)方程，并對(duì)各種蠕變的本構(gòu)方程的適用性進(jìn)行了研究。結(jié)論如下：

(1)開展了某配方HTPB推進(jìn)劑試件的蠕變?cè)囼?yàn)，獲取蠕變數(shù)據(jù)。試驗(yàn)結(jié)果表明，HTPB推進(jìn)劑的蠕變過程符合蠕變一般規(guī)律。

(2)對(duì)多種蠕變本構(gòu)方程進(jìn)行研究：冪律類型蠕變本構(gòu)方程、復(fù)合時(shí)間硬化蠕變本構(gòu)方程和廣義時(shí)間硬化蠕變本構(gòu)方程均適用于HTPB固體推進(jìn)劑的蠕變過程。其中，復(fù)合時(shí)間硬化蠕變本構(gòu)方程的擬合結(jié)果與試驗(yàn)數(shù)據(jù)誤差最小。

(3)冪律類蠕變本構(gòu)方程的擬合曲線與試驗(yàn)曲線在各應(yīng)力水平都具有較好的一致性，能夠反映HTPB固體推進(jìn)劑的蠕變過程，且形式簡(jiǎn)單、處理方便，可在初步分析或特定應(yīng)力水平下使用。本文研究成果可用于固體發(fā)動(dòng)機(jī)立式貯存問題的仿真分析，為預(yù)示發(fā)動(dòng)機(jī)立式貯存過程的蠕變效應(yīng)提供指導(dǎo)。