艇載振動載荷下固體發(fā)動機裝藥有限元分析①

曲 凱，李金飛，陸巍巍，張旭東

(海軍航空大學 岸防兵學院，煙臺 264001)

0 引言

固體火箭發(fā)動機由于其結構簡單、易于維護、可靠性高、推重比大等優(yōu)點，已廣泛應用于各種型號的導彈中。為評估長期服役環(huán)境下固體火箭發(fā)動機的貯存壽命，必須研究溫度、濕度、艦艇振動等因素對發(fā)動機貯存壽命的影響。其中，艦艇振動載荷是決定服役環(huán)境下固體火箭發(fā)動機壽命的一個重要因素。潛射導彈在潛艇中長期經(jīng)歷復雜環(huán)境載荷的作用，準確獲取固體發(fā)動機藥柱在全壽命周期中載荷歷程，是對其進行技術狀態(tài)評估的基礎。目前，獲取發(fā)動機藥柱載荷歷程的方法主要有兩種：一是通過安裝嵌入式傳感器直接測量得到；二是通過測量導彈所承受環(huán)境載荷，然后利用有限元仿真軟件計算得到。

BUSWELL等[1]和GERHARD等[2-3]采用內(nèi)嵌式應力和化學傳感器來監(jiān)測固體火箭發(fā)動機粘接界面在溫度循環(huán)、降溫時的應力變化和化學老化。研究表明，該應力傳感器可以用來研究交變溫度應力作用下，復合固體推進劑粘接界面的損傷規(guī)律。張守誠等[4]和張波等[5]分別采用殼體外表面和內(nèi)表面貼片傳感器，可有效地測量相應的載荷和脫粘情況，實現(xiàn)發(fā)動機局部健康監(jiān)測。雖然嵌入式傳感器已經(jīng)取得了很大進步，但是到目前為止，該方法還沒有大規(guī)模運用到實際發(fā)動機監(jiān)測中，所取得的成果還需要進一步通過解剖發(fā)動機、取出裝藥進行試驗的方法對其破壞方式進行驗證，而且傳感器目前大部分只是針對溫度載荷產(chǎn)生的應力，對于振動載荷所產(chǎn)生的頻率較快的應力測量還沒有投入使用。

國內(nèi)外學者[6-10]針對固體火箭發(fā)動機長期貯存、公路運輸、載艦值班等情況進行了有限元分析研究，并根據(jù)藥柱破壞判據(jù)對其的可靠性進行了判斷，為發(fā)動機壽命評估提供了一定依據(jù)。但是，仿真計算過程中振動載荷輸入條件都為低頻振動載荷，特別是艦船運動周期為10 s，沒有將實測振動數(shù)據(jù)作為輸入條件進行仿真計算。王永帥等[11]針對艦載立式貯存固體發(fā)動機藥柱和粘接界面的蠕變損傷分別展開了研究。王鑫等[12]開展了發(fā)動機固化降溫、長期重力和低頻振動聯(lián)合作用下的有限元計算，得到了藥柱應力應變分布規(guī)律。文獻[6,8]研究表明，雖然艦船搖擺運動所產(chǎn)生的加速度載荷比較小，達不到固體火箭發(fā)動機裝藥破壞的極限載荷，但是這種小應力載荷的長期作用也可能產(chǎn)生較大損傷。因為推進劑藥柱在構成上是多組分混聚物，主要由高分子聚合物基體和摻入其中的大量固體氧化劑顆粒及金屬燃料顆粒組成。由于固體顆粒在基體中是一種機械混合，在外界各種小應力載荷作用下，推進劑內(nèi)部將發(fā)生不可逆的損傷。從細觀層面上講，該損傷包括基體結構的斷裂以及顆粒與基體界面的“脫濕”，而且這些細觀損傷可能會聚合和匯集成宏觀裂紋、脫粘等缺陷[13-14]。

綜上所述，本文利用粘彈性有限元分析方法計算艇載固體發(fā)動機在實測振動載荷作用下的應力應變場，通過分析確定發(fā)動機裝藥比較容易出現(xiàn)問題的部位，并計算研究危險部位Mises應力和剪切應力隨時間變化曲線，最后通過雨流計數(shù)法獲取發(fā)動機裝藥危險部位的疲勞載荷特性，可為后續(xù)開展裝藥實驗室加速老化試驗和艇載固體發(fā)動機壽命評估奠定基礎。

1 振動載荷作用下裝藥結構有限元分析

1.1 三維實體建模

對某型艇載導彈固體火箭發(fā)動機進行三維有限元分析，首先必須建立其相應三維結構模型。本文所研究的固體火箭發(fā)動機由金屬殼體、絕熱層、襯層、翼柱形推進劑、人工脫粘層等部分組成，如圖1所示。

圖1 某型導彈發(fā)動機結構示意圖

論文研究重點是分析發(fā)動機裝藥危險部位的應力應變分布，為提高運算效率，對發(fā)動機剛性殼體結構進行一些必要簡化。根據(jù)某型導彈固體發(fā)動機的結構特點，忽略對應力分析不產(chǎn)生明顯影響的復雜結構，比如發(fā)動機殼體的前裙、后裙和點火器等結構，利用建立有限元仿真軟件的模型構建模塊建立發(fā)動機模型，如圖2所示。

圖2 某型導彈發(fā)動機三維模型示意圖

為更清楚地觀察所建固體火箭發(fā)動機模型的內(nèi)部結構，現(xiàn)在取該發(fā)動機三維模型的1/2進行分析。由圖3可見，發(fā)動機內(nèi)部結構相對比較復雜，它由頭部和尾部開槽的翼柱形推進劑藥柱、人工脫粘、襯層和發(fā)動機殼體組成。該發(fā)動機為了減少應力集中，在藥柱頭部和尾部均設有人工脫粘絕熱層，用來消除頭部和尾部的應力集中。

圖3 某型導彈發(fā)動機1/2建模結構示意圖

1.2 固體火箭發(fā)動機材料參數(shù)的選取

進行有限元分析，還需要提供所建模型材料的基本性能參數(shù)。其中，該型發(fā)動機藥柱的彈性模量通過開展發(fā)動機所使用的推進劑松弛模量測量實驗得到。該固體火箭發(fā)動機的裝藥、殼體、襯層材料參數(shù)見表1。

表1 材料性能參數(shù)

由于這里推進劑裝藥采用線性粘彈材料，有限元分析軟件為了建立本構關系，在計算過程中需要用到體積模量函數(shù)K(t)和剪切模量函數(shù)G(t)，其計算公式按照式(1)和式(2)計算。

(1)

(2)

通過推進劑應力松弛試驗[15]，可獲得裝藥的應力松弛模量函數(shù)Prony級數(shù)形式為

(3)

式中Ee為平衡模量；τi為松弛時間。

瞬態(tài)模量E0可由下式計算：

(4)

根據(jù)應力松弛試驗原理，n、t和τi可以任意選取。這里n取4，t取1、4、40、400 s和600 s，τi取0.05、0.5、5和50。因此，測量推進劑裝藥在各時間點上的松弛模量數(shù)據(jù)，如表2所示。

表2 推進劑松弛模量數(shù)據(jù)

根據(jù)上述條件，式(4)里面有五個未知數(shù)，分別為Ee、E1、E2、E3和E4。將表2中的試驗數(shù)據(jù)代入到式(4)中，可得到5個方程，因此可將5個未知數(shù)求出。最終得到松弛模量函數(shù)曲線，如圖4所示，方程如式(5)所示。

圖4 推進劑裝藥松弛模量曲線

E(t)=1.561+1.685e-t/0.05+2.167e-t/0.5+

3.034e-t/5+3.731e-t/50

(5)

1.3 載荷施加與邊界條件

由于進行有限元計算所需的邊界條件和載荷條件與選取的固體火箭發(fā)動機參考坐標系相關，因此首先確定有限元分析的參考坐標系。為建模研究方便，本文采用的坐標系與潛艇航行坐標系相同，如圖5所示。

圖5 固體發(fā)動機坐標系示意圖

不論選取何種坐標系，在進行有限元分析時，任何微元的受力可分為體積力和表面力。有限元分析中每個微元所受到的應力大小僅取決于微元形變的大小，與所選取的坐標系無關，因為無論在慣性坐標性下還是非慣性坐標系下其物體發(fā)生形變量是相同的。微元所承受體積力與坐標系的選取密切相關，因為在不同的坐標系中微元運動規(guī)律不同，產(chǎn)生的加速度體積力也不相同，而且在有限元分析中，體積力不是未知要求解的載荷，而是作為已知載荷條件輸入到模型當中。

1.3.1 邊界條件

由于某艇載導彈固體火箭發(fā)動機裝載于潛艇的發(fā)射筒中，如圖5所示。為防止發(fā)動機在發(fā)射筒中左右、前后移動，發(fā)動機殼體與發(fā)射筒內(nèi)壁之間裝有適配器。除此之外，為防止發(fā)動上下移動，在發(fā)動機圓筒段與發(fā)射筒內(nèi)壁之間還裝限位剪切裝置。因此，發(fā)動機殼體圓筒段可以近似認為位移為0，即

ui=0(i=r,z,θ)

(6)

由于發(fā)動機前后封頭與發(fā)射筒之間沒有物理連接，因此可視為自由表面。

發(fā)動機裝藥內(nèi)孔，在未點燃燃燒時，其與空氣的接觸邊界為自由表面。

綜上所述，選取的發(fā)動機坐標系建立裝藥有限元分析模型時，施加的邊界條件為

(1)發(fā)動機殼體圓筒段為固定邊界條件；

(2)發(fā)動機前后封頭、裝藥內(nèi)孔表面為自由邊界條件。

1.3.2 載荷條件

發(fā)動機裝藥內(nèi)孔，在發(fā)動機未點火時，徑向所受外力為零，即

σr=0

(7)

由于在潛艇運動過程中，固體火箭發(fā)動機裝藥變形量相對于發(fā)動機的位移運動非常小，可以忽略其影響。因此，可假定發(fā)動機裝藥和發(fā)動機殼體的位移運動規(guī)律相同。

根據(jù)達朗貝爾原理，可求解發(fā)動機裝藥微元在慣性坐標下的體積力，計算公式如下：

(8)

其具體數(shù)值由實測潛艇發(fā)射筒中固體火箭發(fā)動機的振動數(shù)據(jù)確定，根據(jù)文獻[16]研究表明，發(fā)動機所遭受的振動載荷主要集中在振動幅值比較小的區(qū)域內(nèi)，而且通過對載荷頻譜分析表明，10 s載荷的頻譜與1000 s振動載荷的頻譜基本一致。考慮到振動載荷頻率快，載荷步比較小，計算量較大，所以綜合考慮上述因素，以10 s為時間單位作為載荷輸入，如圖6所示。

(a) X axis acceleration-time curve

1.4 網(wǎng)格劃分

根據(jù)有限元分析步驟，還需要對所建立三維實體模型進行網(wǎng)格劃分，采用自由四面體劃分方式，劃分單元總數(shù)為967 888個，為了更好說明發(fā)動機各個組成部分的網(wǎng)格劃分情況，選取其網(wǎng)格的1/2進行展示，如圖7所示。

圖7 某型導彈發(fā)動機三維模型網(wǎng)格劃分圖

完成上述步驟后，就可以利用有限元計算軟件對固體發(fā)動機在振動載荷下進行應力應變分析。

1.5 計算結果

利用有限元計算軟件進行顯式動力學分析，仿真計算機采用8個CPU處理器(每個處理器主頻2.33 GHz)，內(nèi)存24 GB，計算的載荷步為10 000步，步長時間為0.01 s，得到了艇載導彈發(fā)動機藥柱在執(zhí)行戰(zhàn)備值班任務時遭受典型振動載荷0～10 s的應力應變分布。圖8～圖11分別為4.8 s時刻發(fā)動機裝藥Mises應力場、最大主應變場、τyz剪應力場和γvz剪應變場。

圖8 4.8 s藥柱Mises應力云圖

圖9 4.8 s藥柱最大主應變云圖

圖10 4.8 s藥柱yoz平面剪應力云圖

圖11 4.8 s藥柱yoz平面剪應變云圖

根據(jù)計算結果可知，在4.8 s時，裝藥的Mises應力最大值為15.88 kPa，其主應變最大值為4.599×10-3，yoz平面剪應力最大值為3.394 kPa，yoz平面剪應變最大值為3.382×10-3。其中，Mises應力和最大剪應變位置分別位于發(fā)動機藥柱頭部尾部外表面，而yoz平面剪應力和剪應變位置位于發(fā)動機藥柱內(nèi)表面上，而且它們都遠遠低于材料相應的極限破壞值。因此，推進劑裝藥不可能由于瞬時受力超過極限臨界值而發(fā)生破壞。但這種振動載荷對艇載固體發(fā)動機裝藥結構可能發(fā)生累積損傷，而且其破壞形式主要包括以下兩方面：

(1)拉伸、壓縮的交變應力循環(huán)會對藥柱造成機械損傷，降低藥柱力學性能；

(2)另外，振動產(chǎn)生的慣性力會在固體發(fā)動機藥柱粘接界面產(chǎn)生交變的剪切應力，降低粘接界面剪切強度，從而造成粘接界面的破壞。

2 艇載固體發(fā)動機裝藥危險部位的確定

2.1 發(fā)動機裝藥應力分布

危險部位可按照裝藥的破壞形式來確定，關注Mises應力和yoz平面界面最大剪切應力較大的部位。首先定義兩條路徑和一個截面A-A。兩條路徑分別為用紅“▲”點和紅“■”點來表示，并命名為路徑1、路徑2，如圖12(a)所示。路徑1和路徑2分別位于發(fā)動機外表面和內(nèi)表面，并貫通整個裝藥的軸向。發(fā)動機截面A-A靠近發(fā)動機尾部，該位置為發(fā)動機4.8 s時內(nèi)孔應力最大的截面。選取的路徑3為發(fā)動機截面A-A中紅點沿順時針沿裝藥內(nèi)孔行進，并繞內(nèi)孔一周最終回到該點的路徑，如圖12(b)所示。

(a) Paths of 1/2 grain (b) Section A-A

2.2 根據(jù)Mises應力確定危險點

為了研究4.8 s時發(fā)動機的Mises等效應力分布規(guī)律，首先研究發(fā)動機裝藥沿各個路徑上各點的應力分布情況，因為最大應力部位會出現(xiàn)在這些路徑上。每條路徑上Mises應力分布，見圖13。

(a)Mises stress of path 1

由圖13(a)、(b)可見，在發(fā)動機頭部A點和尾部B點、C點Mises等效應力較大，其都沒有位于裝藥的兩端，而位于翼槽和人工脫粘層后面，曲線中間呈現(xiàn)明顯雙峰特點。這是由于頭部和尾部人工脫粘設計釋放了這兩個部位的應力集中。由圖13(c)可見，在裝藥內(nèi)孔表面上，應力最大部位均在翼柱形藥柱的8個星角處，而且位于發(fā)動機裝藥下部翼柱的星角應力較大，考慮對稱性，取C點作為裝藥星角易發(fā)生破壞部位。

根據(jù)上述分析，取Mises等效應力較大的點A、B、C作為發(fā)動機裝藥易于破壞位置來重點研究。

2.3 根據(jù)粘界面最大剪切應力確定危險點

研究4.8 s時發(fā)動機粘接界面的最大剪切等效應力分布場規(guī)律，確定沿發(fā)動機裝藥沿路徑1和路徑2的剪應力分布情況,如圖14所示。可見，在發(fā)動機頭部的A點、發(fā)動機尾部的B點和C點等局部點的界面剪切等效應力較大。

(a)Shear stress of path 1

綜合考慮裝藥最大等效應力和粘接界面最大剪切應力兩個方面，確定發(fā)動機裝藥容易發(fā)生破壞的危險點為A點、B點和C點。計算過程中，記錄這些危險點每一個載荷步的應力應變數(shù)值。

2.4 危險部位的應力曲線

根據(jù)所得仿真計算結果和所確定的發(fā)動機裝藥危險部位，研究A點、B點和C點的Mises等效應力隨時間變化的規(guī)律，如圖15所示。同時，研究粘接界面A點、B點和C點的最大剪應力隨時間變化的規(guī)律，如圖16所示。

(a) A point (b) B point (c) C point

(a) A point (b) B point (c) C point

通過觀察圖15可知，危險點A的Mises應力幅值最大，而C點的Mises應力幅值最小。另外，每一個危險點的Mises應力加載周期差別不大。通過觀察圖16可知,危險點C的剪應力幅值最大，而A點的剪應力幅值最小，但其最大值小于對應的Mises等效應力。不論是Mises等效應力還是剪應力，都遠低于推進劑及相應粘接界面的載荷破壞極限。因此，推進劑裝藥不可能由于瞬時受力超過極限臨界值而發(fā)生破壞。雖然，長期振動載荷引起的載荷較小，但是其作用時間長，其疲勞損傷累積效應予以關注。

3 發(fā)動機裝藥危險部位疲勞載荷特性分析

為后續(xù)研究振動載荷引起的疲勞累積損傷，需要對發(fā)動機裝藥進行疲勞損傷試驗研究。因此，首先需要對危險點A和危險點C分別研究其Mises等效應力和剪應力的疲勞載荷特性，才能設計合理的裝藥結構件損傷試驗，為后續(xù)發(fā)動機疲勞損傷評估奠定基礎。

在各種疲勞特性分析方法中雨流計數(shù)法[17-18]由于其原理與材料疲勞損傷機理相一致，而被廣泛使用。利用雨流計數(shù)法，對危險點A的Mises等效應力載荷特性的分析結果，如圖17所示。把循環(huán)加載Mises應力幅值按照1000 Pa為一個等級，對其進行計數(shù)統(tǒng)計，結果如圖18所示。

圖17 發(fā)動機A點雨流計數(shù)結果曲線

圖18 發(fā)動機A點Mises循環(huán)應力幅值分布曲線

由圖17可見，在0～10 s范圍內(nèi)，發(fā)動機危險點A累計作用94個Mises應力循環(huán)，其中最大循環(huán)應力幅值為15.573 kPa，最大循環(huán)平均應力為8.891 kPa。由圖18中可見，危險點A點Mises應力循環(huán)大多數(shù)集中在小應力幅值對應區(qū)域，比如小于循環(huán)應力幅值低于7 kPa的個數(shù)為73個，占對應應力循環(huán)總數(shù)的77.7%，隨著加載循環(huán)應力幅值的持續(xù)上升，循環(huán)的數(shù)量不斷減少。為了方便說明把Mises循環(huán)應力幅值低于7 kPa定義為小應力加載循環(huán)。

利用雨流計數(shù)法，對危險點C的剪應力載荷特性的分析結果，如圖19所示。把循環(huán)加載剪應力幅值按照1000 Pa為一個等級，對其進行計數(shù)統(tǒng)計，結果如圖20所示。

圖19 發(fā)動機C點雨流計數(shù)結果曲線

圖20 發(fā)動機C點循環(huán)剪應力幅值分布曲線

由圖19可見，在0～10 s范圍，內(nèi)發(fā)動機危險點C累計作用114個剪應力循環(huán)，其中最大循環(huán)應力幅值為10.11 kPa，最大循環(huán)平均剪應力為1.8 kPa，并在0點附近振蕩，遠低于危險點A的最大循環(huán)平均剪應力。從圖20中可以看出，危險點A點平均剪應力循環(huán)大多數(shù)也集中在小應力幅值對應區(qū)域，比如小于循環(huán)應力幅值低于6 kPa的個數(shù)為102個，占對應應力循環(huán)總數(shù)的89.5%，隨著加載循環(huán)應力幅值的持續(xù)上升，循環(huán)的數(shù)量不斷減少。為了方便說明，把剪切循環(huán)應力幅值低于6 kPa定義為小應力加載循環(huán)。

4 結論

對某艇載固體火箭發(fā)動機翼柱型裝藥戰(zhàn)備值班期間振動環(huán)境載荷作用下的應力應變情況進行了有限元模擬，并對裝藥危險部位的疲勞載荷特性進行了分析，研究結果表明：

(1)該發(fā)動機翼柱型裝藥在振動載荷作用下產(chǎn)生的應力應變遠遠小于其破壞極限，不會發(fā)生突然斷裂破壞。

(2)長期振動載荷作用下應力呈現(xiàn)隨機性變化，確定固體發(fā)動機頭部外表面A點遭受的Mises應力最大，藥柱后翼槽危險點C的剪應力幅值最大。

(3)根據(jù)雨流計數(shù)法，危險點Mises等效應力循環(huán)加載幅值大于剪應力循環(huán)加載幅值，而且Mises等效應力循環(huán)加載均值明顯大于剪應力循環(huán)加載均值。

(4)不論是Mises應力還是剪應力，循環(huán)大多數(shù)集中在小應力幅值(Mises應力幅值低于7 kPa，剪切應力幅值低于6 kPa)對應區(qū)域，A點Mises應力和C點剪應力的小應力循環(huán)數(shù)分別占對應應力循環(huán)總數(shù)的77.7%和89.5%。