基于編制載荷譜的固體發動機藥柱疲勞損傷分析

劉磊，李高春，劉著卿，李金飛

（1.海軍航空大學，山東 煙臺 264001；2.91458部隊，海南 三亞 572000）

對于裝載于海上作戰平臺的固體發動機而言，其先后要經歷生產制造、無損檢測、道路運輸、洞庫貯存、海上值班和發射等過程，其中，洞庫貯存和海上值班經歷的時間最長。海上值班時，部分大型固體發動機能夠實現對溫度和濕度較好的控制，但裝載于平臺的固體發動機時刻受到海洋波浪的作用，產生的低頻振動是無法人為控制的。海洋波浪引起的振動是不具有明顯規律的復雜隨機載荷，難以通過數學公式對其進行表述，研究海洋值班條件下振動對固體發動機裝藥帶來的影響往往需要較多的實測數據。由于艦艇上設備存儲容量有限，進行長時間振動數據監測難度較大。

工程上常用載荷譜來進行結構的疲勞設計以及振動疲勞計算。編制固體發動機振動載荷譜主要通過對部分實測載荷數據進行數據處理和統計學分析，獲得長周期內所承受載荷的幅值和頻次信息[1]。自Ernst[2]提出標準程序載荷譜以來，程序疲勞試驗法在工程中得到了廣泛應用。隨著液壓伺服控制精度的提高，利用隨機疲勞載荷譜對不規則載荷進行模擬的技術得到重視[3]。Johannesson[4]在2001年提出了基于雨流計數的非參數密度估計方法，并通過汽車測量信號進行了驗證。高天宇等[5]應用非參數雨流外推方法編制了某型裝載機液壓缸的載荷譜，進而設計了全作業段程序加載譜。李妍[6]使用改進的非參數雨流外推方法建立了輪式裝載機的載荷譜，證明了全帶寬非參數模型外推性能優于單帶寬非參數模型。

本文利用某固體發動機實測振動載荷數據，通過非參數雨流外推方法獲得了長周期條件下發動機循環載荷的幅值-頻次信息，應用雨流重構法構建了固體發動機振動時的域載荷譜。

1 數據采集與預處理

1.1 數據采集

為獲取某固體發動機值班時的環境數據，首先對發射貯運箱內的振動載荷進行實測，采用1C301型電容式三向加速度傳感器采集振動數據。安裝時，將傳感器粘接于發射貯運箱內壁?？紤]到長期監測的需求與存儲容量的限制，對振動數據采用間斷性采集，設置采樣頻率為400 Hz。振動傳感器三個方向的坐標設定符合笛卡爾坐標系，x軸正向表示主航向，y軸正向表示水平面內主航向逆時針旋轉 90°方向，z軸正向表示垂直向上方向。

通過前方數據傳輸，獲得了某固體發動機在海上常規巡航時的6個采集樣本，樣本量見表1。

表1 監測數據Tab.1 Monitored data

1.2 數據預處理

1.2.1 消除趨勢項

趨勢項是指在長期環境溫度等因素影響下造成的傳感器監測數據偏離原始基線的現象，數據處理的第一步就是要消除長周期趨勢項，使信號在零線附近變化。采用較為常用的多項式最小二乘法來消除趨勢項，該方法通過解線性方程組來求出待定系數，具體算法見文獻[7]。

1.2.2 數字濾波

研究表明，艦船在海浪、潮涌等載荷作用下的低頻振動頻率范圍為0～5 Hz[8]。因此，使用帶通濾波器提取監測樣本在 0～5 Hz頻率范圍內的振動信息。以樣本1的監測數據為例，圖1給出了消除趨勢項及數字濾波處理前后的數據結果。

分析圖1可知，原始采集數據由于疊加了高頻噪聲信號而顯得雜亂無序，同時在x軸、y軸采集信號中存在多處極端幅值點。經過消除趨勢項和濾波處理后，x、y、z三軸振動數據呈現出較為明顯的振蕩波形，部分極端大幅值載荷也得到濾除。

從數據處理后的PSD曲線圖可知，在樣本1的采集時間段內，振動信號的主要頻率集中在0～0.5 Hz之間。其中，x、y軸功率譜密度極大值對應的頻率為0.2 Hz和0.4 Hz，z軸功率譜密度極大值對應的頻率為0.3 Hz。

1.2.3 消除奇異點

處理信號時出現幅值遠遠偏離左右兩點的數據點稱為奇異點，奇異點的存在會使振動曲線產生尖銳毛刺。Ncode軟件中的奇異值探測模塊提供了對奇異點的探測和消除功能（如圖2所示），探測到的奇異點可以直接刪除或者用周期信號代替[9]。載荷數據中奇異點的探測方法主要有梯度門限法、幅值門限法和標準方差檢驗法，分別對應軟件模塊中的Differential、Amplitude、CrestFactor三種算法。通常需要使用兩種及以上的方法來清除毛刺，對采集的樣本數據使用多種算法進行處理，共消除奇異點71個。

圖1 原始信號與數據處理后的信號Fig.1 Raw signal and signal with processed data: a) raw data; b) x axial data after processing;b) y axial data after processing; c) z axial data after processing

圖2 軟件濾除奇異點操作過程Fig.2 Filter-out of singular points with software

1.2.4 數據拼接

將經過完全數據處理的 6個樣本的振動數據進行拼接，組成非參數雨流外推的輸入載荷。拼接時，為避免每段載荷連接處出現奇異點，對拼接數據段需要再一次做剔除奇異點處理。剔除奇異點后的拼接載荷段如圖3所示。

圖3 振動加速度載荷Fig.3 Vibration acceleration load

2 載荷譜編制

以拼接載荷數據為基礎，進行振動載荷譜的編制。載荷譜的編制主要包含三個內容：載荷頻次統計、載荷外推和載荷重構。載荷外推是載荷譜編制流程中的重要一環，是指運用數理統計的方法，將測得的短時載荷擴展至全壽命周期載荷。雨流矩陣外推法能夠實現頻次和載荷的雙向外推，其中，雨流矩陣非參數外推法可以避免參數法難以用已知分布擬合雨流矩陣的問題，對于隨機載荷的外推效果更好。因此，采用雨流矩陣非參數外推法進行振動載荷的外推。

2.1 載荷頻次統計

載荷頻次統計一般采用的是循環計數法，該方法將復雜的時間-載荷歷程進行簡化，運用統計學方法，重點研究復雜波形中某些量值出現的頻次。常見的幅值計數法有峰值計數法、雨流計數法、穿級計數法等。其中，雨流計數法能同時統計載荷的幅值和均值，且能較好地表現載荷變化對疲勞損傷的效應，因而得到了廣泛應用。分別對x、y、z三軸“from-to”形式的雨流結果進行統計，得到的雨流直方圖如圖4所示。

圖4 三軸雨流矩陣直方圖Fig.4 Triaxial rainflow matrix histogram: a) x axix;b) y axis; c) z axis

2.2 載荷外推

非參數雨流外推的過程建立在核密度估計方法和蒙特卡洛法之上。核密度估計法本質上是以二維“from-to”形式雨流直方圖上所有的實測數據點為核心、以帶寬h為邊界的概率密度估計算法[10]。非參數雨流外推算法的具體過程分為三步。

1）計算初始估計。

通過固定帶寬的二維核估計方法獲得所有數據點(x,y)處的初始概率f(x,y)。

式中：n為數據點個數；h為初始帶寬，一般通過經驗公式獲得。

2）對雨流矩陣進行自適應核密度估計。

在非參數外推過程中，極端載荷附近的數據點比較稀疏，一般需要設置較大的帶寬，小幅值區域的載荷密集，一般設置較小的帶寬。依據以上思路，使用自適應核密度估計方法對樣本數據的概率密度分布進行重新估計。

式中：λi為自適應帶寬因子；h*為自適應帶寬。

λi和h*是自適應核估計算式中最主要的兩個參數。自適應帶寬因子λi通過第一步中的初始核密度估計來計算，相應的計算公式如式(3)。

與初始估計帶寬h不同，自適應帶寬h*的確定需要以幅值-頻次曲線中的載荷極值Rmax為限制條件，通過迭代運算求解得到。

3）運用蒙特卡洛法進行模擬，得到外推后的雨流直方圖。

進行載荷外推時，需要選取一個合適核函數來進行概率密度估計。Ncode雨流子程序中內置了四種在水平面內投影的核函數：圓形核函數、基于均值的橢圓形核函數、基于幅值的橢圓形核函數及 Epanechenikov核函數，分別如圖5所示。由于Epanechenikov核函數綜合考慮了載荷分布中均值和幅值對外推結果的影響因素，且其概率分布模型與雨流矩陣具有良好的適應性[11]，因此采用Epanechenikov核函數對雨流矩陣進行概率密度分布估計。

經驗理論認為，經歷106次循環載荷，對包括極少發生的高幅值載荷之內的全部載荷具有足夠代表性[12]。將載荷累計頻次外推至106次，外推因子k為：

圖5 四種形式的核函數Fig.5 Four forms of kernel functions: a) circular kernel function; b) elliptic kernel function based on mean value;c) elliptic kernel function based on amplitude;d) epanechenikov kernel function

式中：Ni為某工況原始累積頻次，為某工況外推后的累積頻次。

得到外推因子k后，運行雨流外推子程序，得到x、y、z三軸外推頻次曲線和外推后的雨流直方圖，分別如6和圖7所示。由圖6可知，非參數雨流外推實現了載荷頻次與幅值的雙重外推，外推后的x、y、z三軸的極值載荷分別提高了51%、83%、115%，實現了分段監測數據中未涵蓋的更加嚴酷的振動載荷的量化表示。

由圖7可知，外推前后“from-to”雨流分布圖載荷數據點均主要分布在小幅值的中心區域；外推前較大的載荷幅值主要呈沿幅值和均值的交叉十字分布，外推后十字形分布得到平滑，呈現類似矩形的區域分布，且中心區域的小幅值載荷更加集中。

2.3 載荷重構

對雨流外推后的載荷進行時域外推的方法稱為雨流矩陣模擬法。該方法基于以下基本原則：新序列必須呈現與先前序列相同的雨流計數結果，載荷重構過程相當于載荷循環提取的逆過程，將提取的載荷循環插入到時間序列中。如圖8所示，將幅值劃分為等間隔的幾行，從下到上行數依次遞增，載荷起始點所在位置稱為行數（row），載荷轉折點所在位置稱為列數（col.）。一般稱插入到時間序列的載荷循環為插入循環，ri、ci分別是循環載荷的行數與列數；接收插入循環的時間序列為接收循環，rr、cr分別是其行數與列數。雨流重構主要有如下4種規則[13]。

圖6 外推頻次曲線Fig.6 Extrapolation frequency curve: a) x axis; b) y axis; c) z axis

圖7 外推前后的雨流直方圖對比Fig.7 Comparison of rainflow histogram before and after extrapolation: a) before extraplolation of x axis; b) before extraplolation of y axis; c) before extraplolation of z axis; d) after extraplolation of x axis; e) after extraplolation of y axis; f) after extraplolation of z axis

1）如果ri＞ci、接收循環是有序的“谷-峰”值，則必須滿足：ri≤ci且cr≥ri，如圖8a所示。

2）如果ri＞ci、接收循環是有序的“峰-谷”值，則必須滿足：rr≥ri且cr≤ci，如圖8b所示。

3）如果ri＜ci、接收循環是有序的“峰-谷”值，則必須滿足：ri≥ci且cr≤ri，如圖8c所示。

4）如果ri＜ci、接收循環是有序的“谷-峰”值，則必須滿足：rr≤ri且cr≥ci，如圖8d所示。

Tecware軟件已經將雨流矩陣模擬算法內置于軟件中，對“from-to”形式的雨流矩陣提供了雨流編輯和雨流重構功能。將 2.2小節中得到的外推后“from-to”形式的雨流矩陣文件輸入到 Tecware軟件中，加載Rainflow Reconstruction模塊進行時域載荷重構，重構的時域載荷譜如圖9所示。對外推數據進行雨流重構后，共獲得約5.11×104s的時域載荷數據。其中，x軸載荷振幅在?0.229～0.216 m/s2，y軸載荷振幅在?0.232～0.230 m/s2，z軸載荷振幅在?0.341～0.328 m/s2。

圖8 雨流重構規則Fig.8 Rainflow reconstruction rules: a) ri＞ci, the receive cycle is an ordered "valley-peak" value, b) ri＞ci,the receive cycle as "peak-to-valley" value, c) ri＜ci, the receive cycle as "peak-to-valley" value,d＝ri＜ci, the receive cycle is an ordered "valley-peak" value

圖9 外推時域載荷譜Fig.9 Extrapolated time domain load spectrum

3 有限元仿真與疲勞損傷分析

3.1 有限元仿真

依據建立的振動載荷譜，對固體發動機藥柱在值班條件下的振動疲勞展開分析。進行疲勞計算時，通常需要具備三種要素：關鍵部位的應力幅信息、材料的疲勞性能信息以及疲勞損傷計算模型。固體發動機關鍵部位的應力幅可以通過有限元仿真獲得；材料疲勞性能的獲取通常需要開展往復拉伸試驗。

首先建立某型固體發動機模型，如圖10所示。模型由殼體、襯層、推進劑和人工脫粘層組成，各部件參數以及邊界條件設置參照文獻[7]。將載荷譜中的振動加速度數據作為邊界條件輸入到固體發動機模型中，依據發動機在實際值班環境下的貯存狀態，對模型同時施加豎直向下的重力加速度g。有限元計算結果如圖11所示，可以看出，在振動過程中，藥柱內部應力較小，越靠近粘接界面，應力越大；藥柱應力較為集中的部位在前后封頭附近，其中，后封頭附近部位的應力值最大。提取后封頭附近危險點處的Mises應力幅，圖12給出了部分時間段內的應力時程曲線，在該時間段內，危險點處的應力主要在0.014～0.022 MPa波動。

使用雨流計數法對危險點處的應力幅計算結果進行統計[14]，得到的統計結果如圖13所示。發動機危險點處歷經的循環應力載荷幅值主要集中在 0～0.001 MPa，少部分大幅值載荷分布在 0.001 MPa～0.006 MPa。

圖10 固體發動機模型Fig.10 Solid rocket motor model: a) motor model,b) 1/2 motor model

圖11 振動過程中藥柱Mises應力分布Fig.11 Mises stress distribution of the propellant column during vibration

圖12 危險點處應力時程曲線Fig.12 Stress time history curve at dangerous point

圖13 雨流計數結果Fig.13 Rainflow counting result

3.2 累積疲勞損傷計算

使用雨流計數法對危險點的應力幅進行統計后，就可以進行發動機藥柱的疲勞損傷計算。研究表明，Miner線性累積損傷理論以通用性較高、對隨機載荷作用的材料損傷預測效果較好的優點而得到了廣泛使用[15]。依據Miner線性累積疲勞損傷原理，材料的累積損傷是每個應力循環造成損傷的線性疊加，計算公式如下：

式中：D為累積損傷；Di為第i個應力循環造成的損傷；Ni為第i個應力循環幅對應的疲勞破壞次數。

文獻[16]針對推進劑啞鈴型試件開展了定應力往復拉伸疲勞試驗，擬合得到了試驗應力σ與對應的試件疲勞破壞次數N的關系式：

將式(7)代入式(6)，得到推進劑應力幅與累積損傷的關系公式：

將雨流統計結果代入式(8)，得到承受一個載荷譜周期振動的固體發動機藥柱損傷值為1.059×10?4。以一個載荷譜周期的累積損傷值作為基元，計算得到發動機連續值班6個月（180 d）造成的藥柱損傷（D）為0.0323。

4 結論

1）應用非參數雨流外推方法和時域載荷重構法能夠有效地建立固體發動機海上值班長周期時域載荷譜，可以實現將分段監測數據中未涵蓋的更加嚴酷的振動載荷進行量化表示，提高了載荷譜的精確性。

2）在振動過程中，藥柱內部應力較小，越靠近粘接界面，應力越大；藥柱應力較為集中的部位在前后封頭附近，后封頭附近危險點處的應力幅主要集中在0.014～0.022 MPa范圍內。

3）以一個載荷譜周期的累積損傷值作為計算基元，固體發動機在海上連續值班6個月時，由振動造成的推進劑藥柱累積損傷（D）為0.0323。