王 鑫,高 鳴,伍 鵬，陳思彤，王玉峰

(1.海軍航空大學，山東 煙臺 264001；2.海軍裝備部，北京 100048)

引 言

固體火箭發動機從出廠、貯存至最后發射階段，期間經歷固化降溫、工作內壓以及軸向過載等多重載荷共同作用，在不同工況下甚至聯合載荷作用下發動機藥柱的結構完整性一直以來備受關注，結構完整性在一定程度上影響發動機的結構可靠性及壽命[1-2]。國內外對此進行了較多研究，朱衛兵[3]進行了藥柱在軸向高過載的可靠性分析、在公路運輸隨機振動條件下的累積損傷及可靠度計算；唐國金等[4]基于黏彈性隨機有限元和推進劑加速老化試驗，提出了藥柱概率貯存壽命預估模型，分析了不同貯存期藥柱的結構可靠度，并對其進行了壽命預估；劉中兵等[5]、田四朋等[6]和Long等[7]分別采用三維彈性-線黏彈性模型、三維黏彈性響應面隨機有限元和積分型隨機有限元方法，對不同載荷下固體發動機藥柱的結構完整性和可靠性進行了分析；高鳴等[8]以固體發動機力學性能變化引起藥柱點火瞬間結構可靠性降低為指標，對發動機壽命進行了預估；李毅等[9]將本構模型表示為彈性分項與黏性分項之和，并較好地實現了松弛試驗和拉伸試驗的數值模擬計算與試驗結果相互驗證；Deng等[10-11]先后利用三維熱黏彈性本構模型研究了復合固體推進劑顆粒的黏彈性老化行為及非線性黏彈性行為，并設計了對應的發動機藥柱結構分析模塊，可用于發動機結構分析。

立式貯存是固體發動機較為特殊的一種狀態。由于推進劑為黏彈性材料，除受溫度、內壓、軸向過載等，藥柱長期立式貯存期間由于自重可能會產生蠕變效應，使藥柱產生形變，藥柱下沉，嚴重時將使藥型發生變化，更甚之影響發動機內彈道，對于大型固體發動機更不可忽視。Zhang等[12-13]和李東等[14]開展了在20℃和50℃下5個應力水平的雙基推進劑蠕變實驗，并采用多種模型進行了驗證，結果表明，四參數固體模型可以更好地反映推進劑的蠕變行為；Bihari等[15]采用開爾文-沃伊特模型，研究了HTPB推進劑的黏彈特性；利用動態力學分析儀測試蠕變曲線，建立蠕變黏彈特性的數學關系式；王永帥等[16]對艦載導彈固體發動機蠕變損傷進行研究，得出蠕變占藥柱總變形的60%以上；袁軍等[17]開展了大型固體發動機立式貯存狀態下固化降溫、重力和內壓聯合作用的模擬計算分析以及立式貯存試驗，計算結果與試驗結果趨勢一致。但未對蠕變本構方程、蠕變過程等作具體分析，因此目前在大型固體發動機立式貯存的蠕變效應方面未見相關研究報道。

為了獲取發動機裝藥蠕變參數，本研究開展不同應力水平下的推進劑試件蠕變試驗，擬合Norton本構方程。基于本構參數，進行啞鈴型試件模擬計算驗證；并開展帶前后人工脫粘結構的翼柱型發動機固化降溫和立式自重聯合作用下的模擬計算分析，得到發動機藥柱應力、應變以及位移場，并給出關鍵點及關鍵路徑的應力、位移分布及損傷變化，并將模擬計算結果與某值班返廠發動機拆解測試結果進行對比。

1 蠕變試驗

1.1 試驗方案

本試驗研究對象為某型立式貯存固體火箭發動機藥柱，選用材料為其所用HTPB/AP/Al型推進劑，填充質量分數為91%，試件依照QJ924-85《復合固體推進劑單向拉伸試驗方法》規定執行，試驗在常溫下進行。

蠕變試驗分為3種方案，即分別采用位移傳感器、計時器和試驗機。由于采樣頻率限制，傳感器法獲取數據點連續性不好；由于工作時間限制，試驗機法無法實現低應力水平下試件斷裂；而計時器方法無法獲得試件拉伸量變化規律。因此，采用計時器法獲得不同應力下推進劑試件的蠕變破壞時間，采用位移傳感器法獲取不同應力水平的位移—時間變化規律，輔以拉伸機實現高應力水平下的試驗驗證。

通過對5組試件進行單向拉伸試驗，獲取推進劑的最大抗拉強度為1.05MPa，因此將蠕變應力水平控制在1MPa以下。在試驗過程中，試件中間段橫截面積會隨著蠕變過程逐漸減小，若砝碼質量保持不變，應力會逐漸增大。因此，由于試驗設備限制，試驗過程實際上是定載荷試驗，而非定應力試驗。在此，忽略橫截面積的改變，將試驗看作定應力試驗，后續所涉及的應力均為初始應力。

1.2 試驗結果分析

1.2.1 計時器法

為了準確地尋找蠕變時間與應力水平的關系，探究定應力載荷下推進劑的蠕變損傷，采取計時器的方法記錄不同應力水平下蠕變破壞時間。為避免因樣本數量過少而產生誤差，選取30組試件進行試驗，同一應力水平范圍內取4～5個試件，加載定應力(σ)與蠕變破壞時間(t)對應的散點圖見圖1。

圖1 蠕變試驗數據散點圖Fig.1 ScatterFigure of creep experimental data

通過分析發現，試驗滿足Bills驗證的藥柱蠕變斷裂時間與加載應力之間的對數線性關系[18]，即：

σ=B×lgtf+σt0

(1)

式中：σ為藥柱所受應力；B為損傷指數；tf為蠕變時間；σt0為常數。

試驗數據點與擬合關系如圖2所示。

圖2 蠕變試驗擬合結果Fig.2 Fitting results of creep experiment

對試驗數據進行擬合，得σ=-0.1124×lgt+1.1817，相關系數為0.9809，線性關系良好。

1.2.2 位移傳感器法

采用SDMSA磁致伸縮線性位移傳感器對數據進行采集并實時存儲，一端用細繩通過滾輪將傳感器磁環與砝碼相連，試驗開始后，磁環隨著細繩移動，磁環運動的位移即可轉化為試件的伸長位移，并通過計算機實時輸出與存儲。根據計時器法應力水平范圍，從0.55～0.80MPa之間等間隔(0.05MPa)取6個應力水平，并配置相應砝碼質量。為減小試驗誤差，每個應力水平做3組試驗，選取最接近擬合方程的一組數據作為最終試驗結果，各應力水平下的應變—時間變化曲線見圖3。

圖3 不同應力下的應變—時間曲線Fig.3 Curves of strain vs. time at different stresses

從圖3可以看出，試驗結果符合蠕變一般規律。各應力水平規律一致，即分為應變率隨時間減少的瞬時期、應變率恒定的穩定期和應變率急劇上升的破壞階段，大約分別占據總蠕變時間的15%、80%和5%；最后一個階段在應力水平較低時比較明顯，如圖3(d)、(e)和(f)所示。分析其可能原因如下：外載荷作用較大時，推進劑基體和局部顆粒的粘接界面破壞加劇，引起微裂紋萌生擴展直至宏觀斷裂速度更快，進而蠕變的第三階段不夠明顯。

一般情況下，當應力水平低于某一值時，無論蠕變多長時間，都不會產生破壞。對于大型固體發動機而言，其立式貯存重力載荷的作用不會導致藥柱破壞，因而只關注前兩階段即可。

兩種蠕變實驗方法對應的誤差值見表1。

表1 兩種蠕變試驗結果對比Table 1 Comparison of the two creep experiments results

由表1可知，兩種試驗結果一致性較好，可以相互驗證。

1.2.3 試驗機法

為了驗證傳感器試驗中應變變化規律的準確性，采用MTS臺式微機控制電子萬能試驗機進行指定應力水平下的蠕變試驗，并與前兩種方法進行對比。由于拉伸機工作時間的限制，僅選取0.75MPa和0.80MPa進行試驗。為保持方案一致性，加載方式為力保載，直至試件斷裂。兩種應力水平下試驗機與傳感器方法的對比結果如圖4所示。

圖4 兩種應力水平下位移傳感器法與試驗機法試驗結果對比Fig.4 Comparison of the results obtained by sensor with testing machine at different stresses

由圖4可知，傳感器法與試驗機法獲取的應變—時間曲線一致。進而通過3種方法相互驗證，證明了傳感器方法所獲取應變—時間曲線的準確性。

2 蠕變本構方程擬合及驗證

2.1 方程選取與參數確定方法

一般情況下蠕應變εc可表示為時間t、溫度T和應力σ的函數，即：

εc=f(σ,t,T)=f1(σ)f2(t)f3(T)

(2)

溫度與負載固定的蠕變行為模式應變簡化為：

εc=f(σ,t)=f1(σ)f2(t)

(3)

本研究采用Norton冪次方法進行蠕變分析，蠕變可表示為：

(4)

(5)

根據應變—時間關系發現，當采用Origin對某一應力水平下的三參數進行擬合時，m與擬合相關系數R固定不變，而A和n有多種組合，并且一一對應，但是A和n的變化會影響蠕變模擬計算結果，因此，采用以下方法對三參數進行確定：

(1)將各應力擬合mi取平均值作為本構方程中的m值；

(2)對于推進劑而言，給定特定蠕變時間t，對于應力σi和σj，可以得到以下關系：

(6)

(7)

二式相除，得：

f(σi,t)/f(σj,t)=(σi/σj)n

(8)

(9)

(10)

在通過兩應力求解n時，由于不同應力水平斷裂時間不同，因此以高應力對應的蠕變穩定段時間進行取點。對n各行(或列)求和取均值，即作為該應力水平下的ni。

(3)根據mi和ni值擬合Ai。

2.2 擬合結果

對各試驗數據進行擬合，將得到mi及相關系數Ri，見表2。

表2 方程參數匯總表Table 2 Summary of equation parameters

根據上述方法分別求得n和A，同時列于表2。

2.3 模擬計算驗證

對啞鈴型試件進行三維建模，并分別于模具加載界面上施加不同大小應力，圖5為啞鈴型試件施加的邊界條件及載荷，上端兩弧面施加固定約束，下端兩弧面施加面力，與試驗條件一致。

由于在試驗過程中，細繩與砝碼相連，同時試件下端伸長位移幾近相等，因此，選取圖示O點作為比較點。

圖5 啞鈴型試件邊界與載荷Fig.5 Boundaries and loads of dumbbell specimen

圖6為試件受0.55MPa定應力時的Von-Mises應力云圖和位移云圖，豎直向下為正，試件應力集中部位為兩固定曲面，最大軸向位移為39.18mm。

圖6 試件模擬計算結果Fig.6 Simulation results of specimen

同時，繪制各應力水平下試件底端中心點O的應變—時間曲線，并與試驗結果進行對比分析，如圖7所示。

由于ABAQUS模擬計算過程中，并未設置試件斷裂，因此模擬計算僅能反應蠕變前兩個階段。從模擬計算結果來看，各個應力水平試驗數據與模擬計算數據雖有一定的偏差，但總體趨勢一致性較好，蠕變參數選取較為恰當，其構成的蠕變本構方程可以模擬推進劑及發動機裝藥在恒溫定載荷下的應力應變。從結果看，0.75MPa時m值與平均值相差較大，因此擬合結果相差較為明顯，可能由于試驗過程中試件切割及砝碼裝配問題造成的試驗數據不準確造成的。

圖7 模擬計算與蠕變試驗曲線對比Fig.7 Comparison of the simulation curves and the experimental ones

對于0.75和0.80MPa兩個應力水平，將傳感器試驗數據、試驗機試驗數據和ABAQUS模擬計算數據進行對比，如圖8所示。

圖8 兩種應力水平下3種方法的應變—時間曲線對比Fig.8 Comparison of the strain—time curves obtained by three methods at different stress

3 發動機立式貯存模擬計算分析

某型固體火箭發動機三維建模如圖9所示。由于載荷及邊界條件的對稱性，為節省計算時間，選取1/16模型進行計算，并對模型整體施加重力載荷，方向頭部指向尾部為正，即沿X軸方向為正。邊界條件為對殼體施加固定約束、對稱面施加對稱約束，推進劑蠕變參數選取Norton參數擬合結果，即表2最后一行中的數據。

圖9 發動機結構與邊界載荷Fig.9 Structure, boundaries and loads of motor

3.1 固化降溫載荷模擬計算分析

在進行固化降溫下的發動機有限元分析時，采用基于小變形理論下的積分型本構關系[19]，即

(11)

(12)

式中：ξ和ξ′為等效時間；G(t)和K(t)為剪切模量和體積模量，可用松弛模量E(t)和泊松比μ表示；E(t)采用Prony級數形式表示。

圖10為固化降溫載荷下發動機藥柱的Von-Mises應力和應變云圖，發動機經固化降溫后產生的預應變和預應力不可忽視，并將持續影響后續值班環境的發動機立式貯存。

圖10 固化降溫模擬計算結果Fig.10 Simulation results of curing cooling

3.2 固化降溫和重力聯合載荷

在固化降溫和立式貯存自重載荷耦合作用下，發動機裝藥的應力—應變云圖如圖11所示。

圖11 固化降溫和重力載荷模擬計算結果Fig.11 Simulation results under curing cooling and gravity loading

由圖11可知，藥柱中孔及前人工脫粘前緣應力較大，應力最大部位發生在圖示位置，大小為0.07710MPa，最大應變為5.072%。

圖12為發動機軸向位移和徑向位移云圖，分別記為U1和U2。結果表明，在固化降溫和立式自重載荷耦合作用下，發動機前、后人工脫粘層均有不同程度的張開，前人工脫粘層自分離面至開口端張開位移范圍在9～23mm左右，最大位移處為前人脫前開口端部，位移為22.96mm。后人工脫粘層張開位移值在-10～-5mm左右，最大張開位置為后人工脫粘層前開口端部，最大位移為-9.904mm。藥柱徑向位移最大位置發生在藥柱中孔，為-7.715mm，前人工脫粘層位置處最大徑向位移為5.62mm。

圖12 固化降溫和重力載荷位移云圖Fig.12 Displacement clouds of curing cooling and gravity loading

3.3 發動機返廠測試與模擬計算對比

對值班返廠的某枚導彈進行拆解，并對固體發動機結構尺寸進行返廠測試。結果表明，值班環境下發動機前、后人工脫粘層均有張開，且前人工脫粘層張開縫隙大于后人工脫粘層。與此同時，藥柱總長度與藥柱內孔直徑也比出廠前有一定的變化。此外，對發動機裝藥前后開口部位及前后封頭部位進行高能加速直線探傷，結果表明推進劑/絕熱層界面無脫粘，推進劑內部無氣孔，U型脫粘結構完整，各粘接界面無貫通至襯層的缺陷，并成功經過地面熱試車考核。

模擬計算結果與實際結果對比表明，前后人工脫粘層最大張開位移相對值在15%以內，藥柱長度變化在1%以內，后錐段相交處藥柱內孔直徑誤差值在1%以內。前后人工脫粘層張開最大位移相對誤差值較大，誤差產生的原因可能與沒有考慮藥柱實際貯存環境的溫度變化和值班振動環境等載荷有關，將在后續作進一步分析。

4 蠕變效應分析及損傷計算

4.1 蠕變效應分析

根據圖11和圖12的模擬計算分析，選取應力集中區域和位移變化具有代表性的5個關鍵點和3條關鍵路徑，如圖13所示。其中，A、C點分別代表藥柱前、后人工脫粘層前開口端位置，B點為藥柱中孔某位置，D、E分別為藥柱前后翼槽錐段某點。Path 1、Path 2和Path 3都為藥柱上節點所組成的路徑，分別代表前人工脫粘層-推進劑/襯層界面-后人工脫粘層、前翼槽-中孔-后翼槽，前翼-中孔-后翼。

圖13 發動機關鍵點及關鍵路徑Fig.13 The key points and key paths of motor

圖14(a)為A、B、C3點再經受固化降溫和自重聯合位移作用下的歷程，以A點為例作具體分析，A點固化降溫后軸向位移為15.2701mm，在經歷靜力自重載荷作用下藥柱下沉至18.9996mm，考慮重力的長期作用，即蠕變作用時，藥柱下沉至22.5704mm，考慮蠕變作用的總位移量是不考慮蠕變時的1.2倍。若單獨分析重力載荷的影響，考慮蠕變時，位移量比固化降溫增加了7.3003mm；忽略蠕變時，位移量比固化降溫增加了3.7295mm，前者是后者的1.95倍。隨著時間的不斷增加，蠕變位移也會隨之增加，所以蠕變對立式貯存固體發動機藥柱所帶來的損傷不可忽視。

圖14 關鍵點位移和應力變化規律Fig.14 Displacement and stress variation law of the key points

圖14(b)為D、E兩點的應力歷程，固化降溫后分別產生0.1515MPa和0.0336MPa的預應力，經自重載荷，應力呈現逐漸減小的趨勢。

圖15(a)為藥柱沿不同路徑的應力變化情況，從圖中可以看出，藥柱中孔的應力要高于前后翼及前后翼槽位置處，藥柱界面的應力集中部位為前人工脫粘層附近界面處。圖15(b)為藥柱沿Path 2和Path 3的軸向位移，用U1表示，位移值在前后翼槽處有所區別，藥柱中孔軸向段幾乎重合。

圖15 應力和位移沿不同路徑的分布Fig.15 Stress and displacement distribution along different paths

4.2 蠕變損傷計算

根據蠕變載荷作用原理，蠕變損傷Dc值可表示為試件在應力σi下的加載時間ti與該應力下蠕變破壞時間tf的比值[20]。

圖16(a)表示Path 2和Path 3兩條路徑的損傷分布規律，與應力規律一致。圖16(b)表示D、E兩點蠕變損傷隨蠕變時間的變化規律，蠕變損傷隨蠕變時間呈增長趨勢，D點蠕變損傷為0.49%。

圖16 關鍵點和路徑蠕變損傷變化規律Fig.16 Creep damage variation law of the key points and paths

5 結 論

(1)3種蠕變試驗結果表明，利用位移傳感器開展蠕變試驗是可行的；啞鈴型試件模擬計算結果證明了選取的參數擬合方法的有效性，Norton本構方程可以用來分析立式貯存固體發動機藥柱的蠕變效應。但是由于蠕變試驗中應力并非恒定，后續考慮通過使用補償彈簧的方式以減小試驗誤差。

(2)帶人工脫粘層的翼柱型立式貯存發動機在固化降溫和長期自重載荷的聯合作用下，前后人工脫粘層均張開，且前者大于后者；藥柱前翼錐段有明顯的應力集中，且中孔應力較高，中孔徑向位移最大。

(3)計算結果與返廠拆解發動機測試對比發現，模擬計算結果與實測結果一致性較好。誤差可能與模擬計算載荷與實際值班環境的偏差以及忽略了公路運輸、臥式貯存等環節的影響有關。

(4)發動機長期立式貯存藥柱會產生蠕變效應，帶來蠕變損傷，損傷值為0.49%，但多種載荷引起的累積損傷將大于該值。結果表明，藥柱內部發生破壞的可能性較小，因此應將重點關注在粘接界面上，藥柱的耦合損傷以及粘接界面的分析將在后續作進一步研究。同時長期貯存引起的裝藥變形對發動機內彈道的影響也將是后續研究的主要內容。