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

羅一智，沙寶林，侯 曉

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

符號說明

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

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

t——時間，s

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

T-H——時間硬化方程縮寫

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

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

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

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

a、b、c、d——為簡化方程形式引入的中間量

e、f、r、p——為簡化方程形式引入的中間量

A、B、C——為簡化方程形式引入的中間量

0 引言

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

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

1 試驗

1.1 試驗儀器與方案

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

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

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

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

圖1 蠕變試驗裝置Fig.1 Set-up for creep test

1.2 試驗結(jié)果

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

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

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

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

(1)

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

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

(2)

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)

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

(4)

(2)時間硬化方程

(5)

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

(6)

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

(7)

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

(8)

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

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

(9)

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

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

為統(tǒng)一描述各應(yīng)力水平的蠕變過程，結(jié)合各組試驗數(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所示。最終擬合曲線與試驗數(shù)據(jù)對比如圖4所示，擬合值與試驗值的殘差平方和RSS如表5所示。

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

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

表5 冪律類方程擬合值與試驗值殘差平方和

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

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

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

εcr=-C1σC2e-rt+B

r=C5σC3e-C4/T

(10)

對于某組應(yīng)力水平的蠕變試驗，應(yīng)力相關(guān)項為定值，上式等價于

εcr=-Ae-rt+B

r=C5σC3e-C4/T

(11)

式中A=C1σC2。

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

(a)Values of A

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

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

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

(12)

于是，可簡寫為

εcr=atb+ctd+etf

(13)

其中

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

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

根據(jù)a、c、e與應(yīng)力的冪函數(shù)關(guān)系，對表中數(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

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

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

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

(1)廣義伽羅伐洛方程

(14)

(2)指數(shù)方程

(15)

(3)諾頓方程

(16)

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

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

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

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

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

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

(17)

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

εcr=AtB+Ct

(18)

針對各應(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ù)合時間硬化蠕變本構(gòu)方程適用于此推進(jìn)劑蠕變過程。最終結(jié)果如表7及圖8所示。

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

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

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

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

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

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

(19)

根據(jù)各應(yīng)力水平的蠕變試驗數(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é)合圖形與計算結(jié)果表明，關(guān)于C的冪函數(shù)擬合計算的決定系數(shù)R2<<1，不具備此有理多項式蠕變本構(gòu)方程要求的冪函數(shù)形式，此方程不適用于推進(jìn)劑的蠕變過程。

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

對于某應(yīng)力水平的蠕變試驗以及常溫條件，原方程可簡化為

εcr=ftr

(20)

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

針對各應(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)方程要求，廣義時間硬化蠕變本構(gòu)方程適用于此推進(jìn)劑蠕變過程。最終確定所有參數(shù)的值，結(jié)果如表9所示。

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

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

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

表10 廣義時間硬化蠕變方程擬合值與試驗值殘差平方和

4 結(jié)論

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

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

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

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