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

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

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

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

本文利用某固體發(fā)動機實測振動載荷數(shù)據(jù)，通過非參數(shù)雨流外推方法獲得了長周期條件下發(fā)動機循環(huán)載荷的幅值-頻次信息，應(yīng)用雨流重構(gòu)法構(gòu)建了固體發(fā)動機振動時的域載荷譜。

1 數(shù)據(jù)采集與預(yù)處理

1.1 數(shù)據(jù)采集

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

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

表1 監(jiān)測數(shù)據(jù)Tab.1 Monitored data

1.2 數(shù)據(jù)預(yù)處理

1.2.1 消除趨勢項

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

1.2.2 數(shù)字濾波

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

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

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

1.2.3 消除奇異點

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

圖1 原始信號與數(shù)據(jù)處理后的信號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 數(shù)據(jù)拼接

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

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

2 載荷譜編制

以拼接載荷數(shù)據(jù)為基礎(chǔ)，進行振動載荷譜的編制。載荷譜的編制主要包含三個內(nèi)容：載荷頻次統(tǒng)計、載荷外推和載荷重構(gòu)。載荷外推是載荷譜編制流程中的重要一環(huán)，是指運用數(shù)理統(tǒng)計的方法，將測得的短時載荷擴展至全壽命周期載荷。雨流矩陣外推法能夠?qū)崿F(xiàn)頻次和載荷的雙向外推，其中，雨流矩陣非參數(shù)外推法可以避免參數(shù)法難以用已知分布擬合雨流矩陣的問題，對于隨機載荷的外推效果更好。因此，采用雨流矩陣非參數(shù)外推法進行振動載荷的外推。

2.1 載荷頻次統(tǒng)計

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

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

2.2 載荷外推

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

1）計算初始估計。

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

式中：n為數(shù)據(jù)點個數(shù)；h為初始帶寬，一般通過經(jīng)驗公式獲得。

2）對雨流矩陣進行自適應(yīng)核密度估計。

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

式中：λi為自適應(yīng)帶寬因子；h*為自適應(yīng)帶寬。

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

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

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

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

經(jīng)驗理論認(rèn)為，經(jīng)歷106次循環(huán)載荷，對包括極少發(fā)生的高幅值載荷之內(nèi)的全部載荷具有足夠代表性[12]。將載荷累計頻次外推至106次，外推因子k為：

圖5 四種形式的核函數(shù)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可知，非參數(shù)雨流外推實現(xiàn)了載荷頻次與幅值的雙重外推，外推后的x、y、z三軸的極值載荷分別提高了51%、83%、115%，實現(xiàn)了分段監(jiān)測數(shù)據(jù)中未涵蓋的更加嚴(yán)酷的振動載荷的量化表示。

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

2.3 載荷重構(gòu)

對雨流外推后的載荷進行時域外推的方法稱為雨流矩陣模擬法。該方法基于以下基本原則：新序列必須呈現(xiàn)與先前序列相同的雨流計數(shù)結(jié)果，載荷重構(gòu)過程相當(dāng)于載荷循環(huán)提取的逆過程，將提取的載荷循環(huán)插入到時間序列中。如圖8所示，將幅值劃分為等間隔的幾行，從下到上行數(shù)依次遞增，載荷起始點所在位置稱為行數(shù)（row），載荷轉(zhuǎn)折點所在位置稱為列數(shù)（col.）。一般稱插入到時間序列的載荷循環(huán)為插入循環(huán)，ri、ci分別是循環(huán)載荷的行數(shù)與列數(shù)；接收插入循環(huán)的時間序列為接收循環(huán)，rr、cr分別是其行數(shù)與列數(shù)。雨流重構(gòu)主要有如下4種規(guī)則[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、接收循環(huán)是有序的“谷-峰”值，則必須滿足：ri≤ci且cr≥ri，如圖8a所示。

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

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

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

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

圖8 雨流重構(gòu)規(guī)則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 有限元仿真

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

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

使用雨流計數(shù)法對危險點處的應(yīng)力幅計算結(jié)果進行統(tǒng)計[14]，得到的統(tǒng)計結(jié)果如圖13所示。發(fā)動機危險點處歷經(jīng)的循環(huán)應(yīng)力載荷幅值主要集中在 0～0.001 MPa，少部分大幅值載荷分布在 0.001 MPa～0.006 MPa。

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

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

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

圖13 雨流計數(shù)結(jié)果Fig.13 Rainflow counting result

3.2 累積疲勞損傷計算

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

式中：D為累積損傷；Di為第i個應(yīng)力循環(huán)造成的損傷；Ni為第i個應(yīng)力循環(huán)幅對應(yīng)的疲勞破壞次數(shù)。

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

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

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

4 結(jié)論

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

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

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