張硯珅，強洪夫，王學仁，耿挺京

(火箭軍工程大學 導彈工程學院，西安 710025)

0 引言

固體火箭發(fā)動機在貯存階段發(fā)動機主要受自重載荷以及固化降溫載荷聯(lián)合作用，長期受應力作用可導致推進劑力學性能下降，影響裝藥結構完整性，使發(fā)動機失效[1]。

針對臥式貯存條件下的固體發(fā)動機裝藥結構完整性問題，唐國金等[2]分析了藥柱在長期自重載荷作用下的位移水平，認為以半年為翻轉周期能夠較好消除藥柱積累的蠕變；楊月誠等[3]通過數(shù)值計算得到真空、高壓和實際貯存情況下藥柱界面危險部位應力。但以上研究均未考慮藥柱在貯存條件下?lián)p傷緩慢積累導致的破壞問題，目前針對臥式貯存定期翻轉條件下固體發(fā)動機藥柱損傷分析的相關研究相對較少。

本文基于熱粘彈性本構和累積損傷理論建立了考慮損傷的發(fā)動機藥柱有限元計算模型，對臥式貯存條件下某型固體火箭發(fā)動機定期翻轉進行數(shù)值仿真，獲得藥柱的應力和損傷分布情況，分析了是否翻轉和不同翻轉周期對藥柱損傷的影響，為準確評估定期翻轉對臥式貯存固體火箭發(fā)動機藥柱的影響提供參考。

1 考慮損傷的發(fā)動機藥柱有限元計算模型

本文研究的發(fā)動機藥柱采用丁羥復合固體推進劑，該推進劑是一種含無機氧化劑和金屬添加劑等高填充固體顆粒的高分子聚合物，基體為粘合劑母體。其力學性質主要取決于高分子粘合劑，因而該發(fā)動機藥柱的力學性質主要表現(xiàn)為粘彈性。

1.1 熱粘彈性本構關系

假設固體推進劑為各向同性材料，考慮變化溫度情形，并基于熱流變簡單材料假設，三維線性粘彈性材料本構方程的一般形式為

(1)

假設其泊松比為常數(shù)ν，則

(2)

(3)

ξ和ξ′形式如下：

(4a)

(4b)

其中，aT為溫度移位因子，且滿足如下的WLF方程：

(5)

式中C1和C2為材料常數(shù)；T為當前溫度；Tr為參考溫度。

材料的松弛模量E(t)表達式如下：

(6)

1.2 藥柱累積損傷模型

貯存過程中，發(fā)動機藥柱長期在外載荷作用下，內部結構會緩慢發(fā)生變化，推進劑力學性能下降，這是一個損傷積累的過程。當損傷發(fā)展到一定程度，會產(chǎn)生明顯裂紋，導致裝藥失效。MINER[4]最早提出累積損傷的概念，他在研究金屬材料循環(huán)載荷下的損傷時，提出了線性累積損傷模型：

(7)

其中，t*(σi)為在恒定應力水平σi作用下材料破壞時間。當D=1時，材料破壞。

BILLS等[2]開展了相關試驗，認為累積損傷理論能有效預測發(fā)動機藥柱的貯存壽命。一些學者也在固體發(fā)動機裝藥結構完整性分析時，考慮了藥柱累積損傷的影響[5-7]。針對固體推進劑材料，試驗研究表明，其應力水平與失效時間具有如下關系[8]：

(8)

則針對連續(xù)變化的應力，其損傷可表示為

(9)

其中，N為關于應力的勒貝格范數(shù)：

N=σ0t01/β

(10)

式中σ0為參考蠕變應力；t0為對應的失效破壞時間。

通過開展固體推進劑不同應力加載水平σi下的定應力拉伸破壞試驗，獲得其破壞時間ti，代入式(9)即可擬合得到損傷參數(shù)N和β。

1.3 增量有限元形式

對于粘彈性材料本構方程的結構有限元分析，由于其力學響應不僅與當前的載荷狀態(tài)有關，而且還與整個載荷作用歷史緊密相關，因此這類問題一般需采用增量法進行求解。首先，對1.1節(jié)的本構模型進行離散，在時間區(qū)間[0，t]共劃分n個分析時間步，即

(11)

而折算時間ξ(t)和損傷變量D(t)也劃分為相應的時間步：

(12)

(13)

同理，可以得到離散化的本構方程：

(14)

其中

ΔSij(tm)=Sij(tm)-Sij(tm-1)

(15)

Δσkk(tm)=σkk(tm)-σkk(tm-1)

(16)

在每個時間步結束時，粘彈性應變εve可以通過對應變率求積分的方法求得[9-10]：

(17)

類似地，一個時間步結束時，粘彈性應力的計算通常通過對應力變化率的數(shù)值積分的方法求得：

(18)

式中σve為粘彈性應力。

由于推進劑不具有明顯的方向性，在工程中也常用各向同性損傷模型來近似簡化，在三向應力狀態(tài)下可選用Mises應力σvon作為等效應力來計算推進劑損傷情況。假設等效應力在每個時間步內近似呈線性變化，則D(t)由式(9)可得：

(19)

2 發(fā)動機計算模型

2.1 發(fā)動機幾何構型與有限元模型

本文所研究的固體火箭發(fā)動機長約4 m，外徑1.8 m，內徑0.5 m。因其具體結構較復雜，建模時根據(jù)實際進行了一定簡化：省略發(fā)動機前后裙與噴管結構，假設殼體與藥柱緊密粘接，中間無氣體存在，計算過程中不考慮界面脫粘情況。藥柱前后端均為八翼柱，中孔為圓管形，采用C3D8R單元對發(fā)動機各結構進行網(wǎng)格離散，所建發(fā)動機有限元模型單元總數(shù)227 600個，節(jié)點248 862個，如圖1所示。

(a)Integral structure of the SRM (b)Section structure of the SRM

2.2 材料屬性

發(fā)動機殼體與絕熱層均按線彈性材料進行考慮，藥柱為粘彈性材料，發(fā)動機殼體、絕熱層、藥柱材料屬性如表1和表2所示。

表1 發(fā)動機材料參數(shù)

表2 藥柱松弛模量系數(shù)

推進劑零應力溫度331.15 K，取參考溫度Tr=298.15 K時WLF方程系數(shù)C1=20.804，C2=458.987。推進劑損傷參數(shù)取β=5.455，N=3.514。

2.3 載荷條件

發(fā)動機臥式貯存需要考慮溫度載荷與重力載荷的聯(lián)合作用情況，對于溫度載荷，發(fā)動機由初始零應力溫度331.15 K緩慢降至298.15 K，并保持一段時間，整個過程72 h，降溫曲線如圖2所示。隨后，發(fā)動機水平放置進行臥式貯存。貯存期間重力方向與發(fā)動機翻轉方向如圖3所示，從彈體頭部向尾部看，在每個翻轉周期結束后沿發(fā)動機中軸線翻轉90°，依次將發(fā)動機I-IV象限轉至豎直向下的重力位置，并如此循環(huán)。

圖2 發(fā)動機溫度曲線

圖3 發(fā)動機翻轉方向示意圖

3 計算結果分析

按照固體發(fā)動機貯存時間和翻轉周期的不同進行取值，參數(shù)條件設置如表3所示，共分為6個計算工況。其中，工況1為不翻轉貯存0.5 a，用于分析臥式貯存條件下發(fā)動機藥柱的力學響應規(guī)律；工況2、3分別為不翻轉、定期翻轉條件下貯存2 a，計算結果用于分析是否翻轉對藥柱損傷的影響；工況4～6為不同翻轉周期條件下貯存8 a，計算結果用于分析翻轉周期對藥柱損傷規(guī)律的影響。

表3 計算工況參數(shù)條件設置

3.1 臥式貯存下應力計算結果

圖4給出了臥式貯存0.5 a后的藥柱外表面Mises應力分布情況。從圖4(a)可以看出，應力水平相對較大的區(qū)域位于藥柱上部外表面前后端人工脫粘層根部位置附近，其中III象限藥柱前端人工脫粘層根部位置處應力水平最大(0.787 MPa)。這是由于藥柱膨脹系數(shù)遠大于殼體，溫度下降后在人工脫粘層根部附近受拉力作用，而III象限藥柱人工脫粘層根部位置在重力載荷作用下同樣產(chǎn)生拉力作用，載荷疊加導致該位置出現(xiàn)應力集中。圖4(b)給出了藥柱前端局部的Mises應力分布?？梢姡瑧^大區(qū)域主要分布在藥柱表面，并沿人工脫粘層根部位置呈周向分布。

(a)Mises stress distribution of the whole structure

藥柱沿XOY對稱面的剖面Mises應力分布如圖5所示?？梢?，除了前后人工脫粘層根部位置有應力集中外，藥柱中孔內表面整體應力水平也相對較大，其中內表面中部位置最大Mises應力為0.040 8 MPa，這主要是降溫過程藥柱收縮產(chǎn)生的熱應力。

圖5 臥式貯存0.5 a后藥柱剖面Mises應力分布

根據(jù)計算結果，如圖6所示，選取III象限藥柱剖面上的內輪廓線path1、外輪廓線path2以及藥柱前端人工脫粘層根部附近環(huán)形區(qū)域path3等3條路徑，對其應力分布情況進行進一步分析。

(a)Distribution position of path1 and path2 (b)Distribution position of path3

藥柱沿特征線的Mises應力分布情況如圖7所示。由圖7(a)可知，在臥式貯存條件下，path1路徑上500～3700 mm段(即中孔區(qū)域)應力水平相對較大，最大值0.040 8 MPa；path2路徑上大部分區(qū)域應力水平小于0.02 MPa，但在790 mm和4136 mm處(即前后人工脫粘層根部位置)存在峰值，分別為0.078 7 MPa和0.064 8 MPa。圖7(b)給出了藥柱沿path3路徑的Mises應力水平分布曲線?？梢姡琁II象限中間位置人工脫粘層根部附近藥柱應力水平最大，為0.078 7 MPa，最小應力水平位于I象限中間位置，為0.026 3 MPa。

(a)Distance vs Mises stress curves on path1 and path2 (b)Distance vs Mises stress curve on path3

綜合分析可知，臥式貯存條件下，推進劑藥柱的Mises應力遠低于其載荷破壞極限。因此，推進劑裝藥不可能由于瞬時受力超過臨界值而發(fā)生破壞。然而，由于載荷作用時間長，藥柱的損傷累積效應不可忽視。其中，藥柱前端人工脫粘層根部位置整體應力水平相對較大，在貯存過程中更容易發(fā)生損傷積累，需重點關注。

3.2 臥式存儲累積損傷計算結果

按照工況2、3的參數(shù)設定，分別對不翻轉和定期翻轉的固體發(fā)動機進行仿真計算，獲得對應的藥柱累積損傷結果。圖8為臥式貯存2 a后藥柱損傷場分布。圖8(a)為不翻轉條件下的藥柱損傷場，其累積損傷較大區(qū)域位于III象限藥柱前端人工脫粘層根部位置，呈周向分布，最大值為6.23%；圖8(b)為0.5 a定期翻轉條件下的藥柱損傷場，累積損傷較大區(qū)域沿藥柱前端人工脫粘層根部呈環(huán)形分布，最大值1.93%，為不翻轉條件下的31.0%。

對藥柱前端人工脫粘層根部位置的損傷情況進行進一步分析，藥柱損傷沿path3的分布曲線如圖9所示。從分布曲線可以看出，在不翻轉條件下，藥柱損傷主要在III象限區(qū)域積累，最大值達到6.23%；而II、IV象限區(qū)域藥柱損傷相對較小，在0.07%～3.58%之間，I象限區(qū)域藥柱損傷最小，在0.07%以內；當采取定期翻轉的貯存條件時，藥柱損傷在各象限近似均勻分布，損傷值在1.73%～1.93%之間。

根據(jù)上述分析，取藥柱最大損傷值對應的P點作為發(fā)動機藥柱易于破壞位置來重點研究，該點位置如圖8(a)中標注所示。在臥式貯存2 a中，P點損傷隨時間增長曲線如圖10所示?？梢钥闯觯诓环D條件下，其累積損傷值隨貯存時間延長近似呈線性增長，這是由于貯存過程中，該點位置藥柱應力水平基本恒定的原因；而在定期翻轉條件下，P點的累積損傷值在第1個翻轉周期內迅速增長至1.56%，而在隨后3個翻轉周期中轉為緩慢增長，2 a后達到1.93%。

(a)Unturned motor

圖9 臥式貯存2 a藥柱累積損傷沿path3分布曲線

圖10 臥式貯存2 a藥柱最大損傷點的損傷增長曲線

3.3 翻轉周期對藥柱損傷影響

為研究翻轉周期對藥柱損傷規(guī)律的影響，按照工況4～6的參數(shù)設置，對固體發(fā)動機進行仿真計算。圖11給出了臥式貯存8 a后藥柱累積損傷場分布?？梢?，不同翻轉周期條件下，藥柱損傷分布規(guī)律基本一致，累積損傷較大區(qū)域沿藥柱前端人工脫粘層根部呈環(huán)形分布，最大損傷點位于III象限藥柱前端人工脫粘層根部位置，為7.71%。圖12中顯示了藥柱沿XOY對稱面的剖面累積損傷分布。其中，較大累積損傷區(qū)域主要分布在藥柱外表面人工脫粘層根部位置。

(a)Tt=0.5 a

臥式貯存8 a中，藥柱危險點P的累積損傷值隨時間增長曲線如圖13所示。其中，在貯存2 a后，2 a翻轉條件下的損傷值最大，達到6.24%，而采取1 a翻轉和0.5 a翻轉條件下的損傷值分別為3.49%和 1.96%；貯存4 a后，采取2、1、0.5 a翻轉條件的損傷值分別為6.97%、3.90%和3.90%；貯存8 a后，這3個計算工況中P點損傷值十分接近，在7.71%左右。

圖13 臥式貯存8 a藥柱最大損傷點累積損傷曲線

這是由于在貯存初期，當翻轉周期較大時，發(fā)動機III象限藥柱前端人工脫粘層根部位置長期處于較大應力狀態(tài)，損傷積累相對較快；而當翻轉周期較小時，藥柱重力載荷及內部應力分布頻繁變化，藥柱最大損傷積累也相對較慢。但在整個貯存過程結束后，采取0.5、1、2 a翻轉條件下，藥柱損傷分布規(guī)律基本一致。

4 結論

(1)臥式貯存條件下，藥柱前后端人工脫粘層根部附近和中孔內表面應力水平較大，藥柱應力危險點位于III象限藥柱前端人工脫粘層根部附近。

(2)臥式貯存不翻轉時，藥柱損傷主要發(fā)生在藥柱前端人工脫粘層根部III象限區(qū)域；而定期翻轉時，藥柱損傷發(fā)生在藥柱前端人工脫粘層根部環(huán)形區(qū)域。

(3)臥式貯存條件下，采取定期翻轉的策略時，藥柱最大累積損傷為不翻轉情況的31.0%。

(4)以0.5、1、2 a為翻轉周期，對固體火箭發(fā)動機進行水平臥式貯存8 a，藥柱損傷均發(fā)生在藥柱前端人工脫粘層根部環(huán)形區(qū)域，損傷分布規(guī)律基本一致。