羅一智，沙寶林，侯 曉

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

符號說明

σ——應力，MPa

εcr——蠕變應變

t——時間，s

S-H——應變硬化方程縮寫

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

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

M-S-H——改進應變硬化方程縮寫

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

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

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

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

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

0 引言

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

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

1 試驗

1.1 試驗儀器與方案

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

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

表1 典型HTPB推進劑配方

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

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

1.2 試驗結果

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

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

2 蠕變本構方程

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

(1)

對方程進行積分、移項處理，將方程化為

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

(2)

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

表2 蠕變本構方程

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

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

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

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

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

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

3 蠕變本構方程的適用性

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

3.1 冪律類型蠕變本構方程

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

(1)應變硬化方程

(3)

對方程進行積分、移項、合并等處理，方程可化為

(4)

(2)時間硬化方程

(5)

對方程進行積分、移項、合并等處理，方程可化為

(6)

(3)改進時間硬化方程

(7)

(4)改進應變硬化方程

(8)

對方程進行積分、移項、合并等處理，方程可化為

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

(9)

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

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

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

a=0.02042×σ2.63112，R2=0.94306

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

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

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

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

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

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

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

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

εcr=-C1σC2e-rt+B

r=C5σC3e-C4/T

(10)

對于某組應力水平的蠕變試驗，應力相關項為定值，上式等價于

εcr=-Ae-rt+B

r=C5σC3e-C4/T

(11)

式中A=C1σC2。

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

(a)Values of A

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

3.3 廣義格雷厄姆蠕變本構方程

對原方程進行積分，并取C8=0，可化為

(12)

于是，可簡寫為

εcr=atb+ctd+etf

(13)

其中

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

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

根據(jù)a、c、e與應力的冪函數(shù)關系，對表中數(shù)據(jù)進行冪函數(shù)擬合，結果如下：

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

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

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

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

(1)廣義伽羅伐洛方程

(14)

(2)指數(shù)方程

(15)

(3)諾頓方程

(16)

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

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

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

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

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

3.5 復合時間硬化蠕變本構方程

(17)

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

εcr=AtB+Ct

(18)

針對各應力水平擬合確定A、C的值，再確定A、C與應力的冪函數(shù)關系，結果如下：

A=0.01351×σ1.64377,R2=0.99383

C=0.0172×σ18.82979,R2=0.99645

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

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

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

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

表8 復合時間硬化蠕變本構方程擬合值與試驗值殘差平方和

3.6 有理多項式蠕變本構方程

有理多項式蠕變本構方程具有極其復雜的形式，本文對其簡化形式進行研究。對方程積分簡化后：

(19)

根據(jù)各應力水平的蠕變試驗數(shù)據(jù)確定C、b、p值，再研究分析C、b、p與應力的關系。C、b、p結果如圖9所示。

(a)Values of C

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

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

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

3.7 廣義時間硬化蠕變本構方程

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

εcr=ftr

(20)

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

針對各應力水平，擬合確定f、r的值，進一步確定f、r與應力的關系，結果如下：

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

R2=0.92548

r=0.19786+0.2221×σ,R2=1

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

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

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

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

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

4 結論

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

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

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

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