姚熊亮, 熊幫虎, 王志凱, 楊娜娜, 張文啟, 萬澤臣
(哈爾濱工程大學 船舶工程學院, 黑龍江 哈爾濱 150001)
艦船在服役期間,會受到嚴酷的海況以及各種水下武器如水雷、魚雷等的攻擊,這些水下武器的攻擊會對艦船彎曲外板、龍骨翼板及艙室設備等造成嚴重損傷。根據水下武器與艦船距離的不同,可將艦船水下爆炸沖擊分為兩大類:接觸爆炸和非接觸爆炸。水下接觸爆炸會造成艦船局部破損,導致艦船失衡;在水下非接觸爆炸的情況下,會對艦船動力系統和設備等造成損傷,導致艦船失去動力。艦用燃氣輪機作為艦船的心臟,不僅要擁有優良的動力性能,還需要擁有較強的抗沖擊性能。目前對于燃氣輪機的抗沖擊研究主要集中在結構設計,對燃氣輪機的強度設計、爆炸沖擊波載荷在燃氣輪機內部的傳遞相關研究較少。為得到水下爆炸沖擊載荷作用于艦船或設備的傳遞特性,謝耀國等基于小波分析,對水下爆炸自由場壓力載荷和壁壓載荷進行時頻域特征分析,得到載荷所包含的頻域信號并統計了載荷在每個頻段的能量。李麗萍等采用小波包對爆炸原始信號進行分析,并建立起沖擊波能量與不同類型炸藥、爆距與毀傷目標固有頻率之間的關系。
目前得到信號能量的方法多為快速傅里葉變換、小波分析等。快速傅里葉變換多為分析信號為周期性的平穩信號,小波分析在處理非平穩的信號時具有很強的優勢而且具有多辨率特性,但選擇小波基有異,處理的結果會也不同,而且有限長的小波基也會造成能量的泄露。爆炸沖擊載荷具有突變快、持續時間短、含多種噪聲干擾的強非線性信號。因此,以上方法不再適合爆炸載荷能量的提取。
本文用經驗模態分解(EMD)方法和與之對應的希爾伯特- 黃變換(HHT),對爆炸沖擊波載荷在燃氣輪機中的傳遞特性進行分析。Jin等利用希爾伯特- 黃變換法對節理梁的非線性特性進行研究,并通過對比加速度信號和位移信號驗證了數據的正確性。周鵬運用經驗模態分解法對比分析了水下爆炸沖擊波載荷在單體船和三體船在型寬及型深方向上的分布規律和傳遞特性,最后得出水下爆炸沖擊波信號中高頻和低頻與艦船沖擊響應的關系,表明EMD方法可以分析載荷分布規律和傳遞路徑。相較于艦船研究,現如今對艦船設備的載荷傳遞特性研究較少,故本文擬采用經驗模態分解法對爆炸沖擊載荷在燃氣輪機中的傳遞特性分析,并對燃氣輪機進行力學模型簡化,以驗證傳遞特性的正確性。本文研究內容對燃氣輪機抗沖擊設計具有參考價值。
1.1.1 燃氣輪機關鍵部件類型確定
隨著我國海軍的迅速發展,高性能動力裝備已顯得愈發重要,燃氣輪機在艦船中應用可提高經濟性,在一定的燃料供應下可以提高船舶的航程,其燃氣輪機在艦船上的應用已經成為現代化艦船的重要標志。
燃氣輪機基座是燃氣輪機強有力的支承構件,它的抗沖擊特性是國內外科研人員研究的重點。根據韓璐等的試驗結果,基座與機匣連接處位置在受到沖擊后易發生斷裂,是抗沖擊的薄弱部位。機匣作為燃氣輪機與艦船相連接的主要承力構件,承受由燃氣輪機工作時產生的各種氣體載荷、沖擊載荷、渦動載荷以及慣性載荷等的作用,故機匣不僅僅要滿足支承功能,其剛度、強度和可靠性以及穩定性也需要滿足需要。轉子及軸在外載荷和內載荷雙重的沖擊下,極易發生徑向跳動,產生彎曲現象。這些情況的發生極易使燃氣輪機發生故障進而使艦船喪失動力系統。而在仿真計算時,由于傳動軸系統零部件多且結構復雜,在建模、設置邊界條件及接觸的過程中,會面臨計算復雜與工作量大的情況。綜合上述分析,本文明確燃氣輪機關鍵部件,主要包括:燃氣輪機基座、燃氣輪機機匣及轉子軸系。根據以上關鍵部件,建立準確可靠的數值模型對于計算燃氣輪機在水下爆炸沖擊載荷的作用下的動力響應并基于此研究載荷傳遞路徑及系統動力學特性而言是非常關鍵的。
1.1.2 關鍵部件力學等效方法
燃氣輪機作為一件復雜的船用設備,其內部由大量部件構成,而各個部件間存在復雜的連接方式,同時構件有復雜的曲面,還存在減振器等具有非線性因素的構件,按實際結構建立燃氣輪機模型非常困難,即使建立完整的燃機模型,對燃氣輪機進行有限元網格劃分的質量也不理想。在進行計算時,這樣低質量的有限元網格會導致計算難以收斂,且對不需要關注的部件也進行精密計算,造成計算資源的浪費。故需要對燃氣輪機實際結構簡化,考慮燃氣輪機重要部件,對不需要關注的部件可以作為質量點加載在有限元模型上,根據實際結構的力學特點對燃氣輪機做適當的簡化,并建立易于有限元網格劃分及計算的幾何模型。
簡化幾何模型應遵循兩條基本原則:1)保留模型的主要尺寸與基本特征;2)能夠使得劃分的網格質量有明顯改善。模型簡化后只可以近似反映工程的實際情況,模型正確是有限元計算結果正確的關鍵,故要充分考慮到結構的實際情況。此外,在建立有限元模型時,在滿足一定精度的同時不要劃分過多的單元,減少計算時間,即要具有良好的經濟性能。
根據燃氣輪機轉子系統的三維模型可知,轉子葉片、輪穀和部分軸段需要簡化。運用轉動慣量互等的原理,將葉片簡化成圓柱,空心圓柱轉動慣量參數表示如圖1所示,具體計算公式如下:

圖1 空心圓柱轉動慣量參數表示Fig.1 Parameter representation of moment of inertia of hollow cylinder

(1)
式中:為輪盤葉片數;為輪盤單個葉片質量;為輪盤葉片長度;為葉片重心到轉軸中心距離;為設計圓盤外環總質量;為設計圓盤外環外徑;為設計圓盤外環內徑。
滾動軸承通過鼠籠彈性支承與軸承座彈性連接,滾動軸承- 鼠籠組合支承如圖2所示,帶有軸的滾動軸承固定在鼠籠內壁上,鼠籠固定在軸承座上,其中與分別為滾動軸承剛度和鼠籠剛度,為滾動軸承和鼠籠的等效質量,為彈性支承相應的剛度,為軸承和鼠籠等效軸承的位移,為軸徑的位移。

圖2 滾動軸承- 鼠籠組合支承Fig.2 Combined bearing of rolling bearing and squirrel cage
由圖2(b)簡化為圖2(c)的彈性支承模型,其等效質量運動微分方程可以表示為

(2)
設方程的解為
=ei
(3)
式中:為運動幅值;i為復數表示的虛數;為轉子渦動頻率。則

(4)
建立了軸承座中心坐標與軸徑中心之間的關系后,與圖2(c)聯立,彈性支承相應的剛度可以求得,即

(5)
不是一個常數,與轉子的渦動角速度有關。滾動軸承剛度在10N/m級別,鼠籠剛度在10N/m級別,質量較小,在低速轉動時,總剛度的變化不大,因此可以把滾動軸承和鼠籠的剛度等效為各部件串聯后的剛度,串聯后的剛度系數為

(6)
根據文獻[19],可選取滾動軸承剛度為9882 1×10N/m,鼠籠支承剛度為1382 1×10N/m,將剛度值代入剛度公式,可得鼠籠彈性支承剛度為1211 8×10N/m。
當艦船遭受到水下爆炸載荷沖擊作用時,設備基礎會受到加速度沖擊。在對設備進行瞬態沖擊響應分析時,通常采用實測的加速度時歷信號或者標準加速度時歷曲線作為設備的沖擊輸入。本文燃氣輪機使用的沖擊加速度信號是根據德國軍用標準 BV 043/85,使用組合三角波的形式對其基礎進行加載。根據設備安裝位置不同,選擇相應的抗沖擊標準。由于設備的沖擊環境極其復雜,設備的沖擊載荷會以響應譜的形式表達。如果把沖擊信號以頻域表示,則如圖3所示設計譜,左右拐點分別為左頻率和右頻率,低于左頻率則為低頻響應,此時以激起艦船及設備位移為主,在左右頻率中間則以激起艦船及設備速度為主,高于右頻率則以激起艦船及設備加速度為主。典型沖擊響應譜如圖3所示,典型三折線譜如圖4所示,低頻段的直線代表等位移段,用譜位移來表示,中頻段的直線為等速度段,用譜速度來表示,高頻段的直線代表等加速度段,用譜加速度來表示,三條直線相交,便形成了三折線設計譜。

圖3 典型沖擊譜曲線Fig.3 Typical impact spectrum curve

圖4 典型設計沖擊譜Fig.4 Impact spectrum of typical design
根據規范規定及設備的安裝情況,可查得設備的譜位移值、譜速度值和譜加速度值。當艦載設備的質量大于5 t時,需對參數進行折減計算:

(7)
式中:為隔離安裝的設備質量(t);為質量常數,恒等于5 t;表示折減后的譜加速度值(g);表示折減后的譜速度值(m/s)。
在使用時域分析方法對燃氣輪機進行試驗時,在試驗中測得的時間加速度曲線具有很大的隨機性,對結果的分析也會產生很大的影響。故選用德國軍用標準BV 04385規范將規定譜值轉化為雙三角波時間歷程曲線,如圖5所示。圖5中,軸表示加速度(m/s)、軸表示時間(s),、分別表示與時刻的加速度,、分別表示正三角波積分面積和負三角波積分面積,為正三角波最大值時間、為正負三角波交界時間點、為負三角波取值最大時間點,為負三角波作用時間結束點。

圖5 正負三角波歷程曲線Fig.5 Positive and negative triangul wave history curves
其轉化關系如下:
1)第1個正三角形的加速度峰值,約為最大加速度的06倍;
2)第1個三角形的面積約為最大速度的075倍;
3)第2個負三角形面積應與第1個三角形面積大小相等,從而使基礎的最終速度為0 m/s;
4)對加速度歷程兩次積分、得到位移,此位移比設計譜的最大相對位移要大105倍;
5)再選擇和,滿足=04和-=06(-)。
根據以上關系,可得到組合三角波各參數的計算公式為

(8)
式中:為沖擊譜譜位移。
通過(8)式將沖擊譜的譜值轉化為組合三角波時間歷程曲線,得到的計算輸入載荷如表1所示。

表1 組合三角波時歷曲線Table 1 Time history curve of combined triangle wave
選取某型燃氣輪機為研究對象,燃氣輪機的結構支承座分別為前支承、后支承和中間定位支承,其有限元模型如圖6所示。

圖6 燃氣輪機支承有限元模型Fig.6 Finite element model of gas turbine support
模型包括支承系統結構、燃氣輪機機匣和轉子系統。因為有些原因,只給出支承系統結構。
因為沖擊波載荷在箱裝體下部近似為均勻分布,故選取燃氣輪機前基座處剖面位置進行載荷傳遞計算,在選擇測點位置時避免焊接位置等應力集中點。各測點位置的分布具體如圖7所示。使用相同的工況進行計算,載荷加載在箱裝體正下方,載荷加載符合德國軍用標準BV 085規定。

圖7 測點分布位置示意圖Fig.7 Distribution of measuring points
應力波在艦船及設備中傳遞引起結構的響應,設備中不同位置處的響應即可反映沖擊波載荷的傳遞情況。因為設備自身的屬性以及設備基座等的強非線性特征,設備的沖擊響應信號也具有強非線性特征,傳統的快速傅里葉變換對此等信號的處理就變得不那么有優勢。故本文運用HHT對燃氣輪機的加速度響應信號進行時頻域分析。該方法主要包含完整集合經驗模態疊加法(CEEMD法)以及希爾伯特變換。
EMD是Huang等提出的一種時頻域分析方法,其優勢在于處理和分析非平穩信號。該方法已被用到設備診斷、生物醫學、結構分析等領域的研究中。EMD方法是根據經驗把復雜的數據轉換成有限個簡單數據即本征函數(IMF)的疊加。對于時域信號()首先判別出其極大值、極小值,分別用三次樣條曲線組成的包絡線將其連接起來,確保樣條曲線均經過這些極值點,取包絡線均值點,則()與之差記為,即
()-=
(9)
若滿足本征值模態函數,則記為1,否則需要進一步分解。重復進行上述步驟,對進行新一輪分解,以此類推,直至分解到滿足判別式:

(10)
式中:表示對原函數的第次分解。記第1個本征函數1,將函數分解為若干個本征模態及單調的殘差:
()=1+2+…++()
(11)
式中:為滿足(10)式時的取值;()為殘差。
對仿真數據()做EMD,其分解過程如圖8所示,其中表示信號振動幅值。

圖8 仿真信號EMD過程Fig.8 Simulating the signal EMD process
限于篇幅,下面以箱裝體底板測點做EMD及能量分析:
1)提取時間加速度時域信號,并畫出圖9所示時間加速度曲線。

圖9 時間- 加速度曲線Fig.9 Time-acceleration curve
2)對步驟1中圖9所示的時間加速度信號做EMD,分解后的IMF各分量如圖10所示,得到的時間加速度信號分解為10個IMF分量,并依次按順序排列。其中1、2、3、4具有較高的頻率,隨著分解的慢慢遞進,得到的IMF分量的頻率逐漸減小,而波長呈現增長趨勢。

圖10 加速度信號的EMDFig.10 EMD of acceleration signal
分解出的本征函數大都具有清晰的物理意義。為了解載荷沿裝備基座到設備內部關鍵部件的載荷傳遞情況,對各個本征函數進行能量分析,其計算公式為

(12)
式中:為第個IMF分所對應的能量;()為第個IMF分量;為各個IMF分量各個離散點的幅值;為各個IMF分量離散點的數量。對本征函數能量計算后得到的結果如圖11所示,并且能量的分布隨著頻率的降低能量有變大的趨勢。經過快速傅里葉變換得到每個函數的頻域信號,可以得到5~10的中心頻率分別為329 Hz、179 Hz、139 Hz、59 Hz、29 Hz、19 Hz。這與艦船設備在水下爆炸載荷尤其是沖擊波載荷作用下呈現低頻響應的理論相符合。10是分解后頻率更小的分量,可能是信號固有的,也可能是其他情況下引起的。

圖11 本征函數能量分布圖Fig.11 Energy distribution of eigenfunction
在研究載荷沿燃氣輪機由基座向內部轉子時,分別選取前中基座各測點向渦輪內部傳遞。選取測點1、2、3、4、5,代表載荷向渦輪軸與轉子軸載荷的傳遞特征。對各測點加速度數值進行經驗模態分解,同時分析其各個IMF分量的主頻率m及其能量占比,計算結果如表2所示。

表2 燃氣輪機加速度測點IMF分量主頻率及其能量占比Table 2 Main frequency and energy ratio of IMF component at acceleration measurement points on the gas turbine
由表2中對艦載燃氣輪機沿前基座傳遞路線的結果中可以看出:加速度載荷在箱裝體底部所激起的高頻響應占總能量較高,約為總能量的0.98;隨著載荷沿燃氣輪機向內部轉子的傳遞,受阻尼和材料固有屬性的影響,載荷所激起高頻響應逐漸減低;載荷從燃氣輪機底部向內部傳遞的過程中,雖然總的能量在衰減,但是高頻部分的占比逐漸降低,燃氣輪機轉子軸系相較于基座及箱裝體底部的響應卻很小,這是因為燃氣輪機的前幾階模態都以低階模態為主,爆炸沖擊波載荷所引起設備振動的部分為低頻部分,而爆炸沖擊波載荷在燃氣輪機的傳遞過程中主要以高頻部分存在。
結合載荷在燃氣輪機中的載荷傳遞特征,可以給出沖擊波載荷在燃氣輪機中傳遞路徑,其示意圖如圖12所示。結合云圖給出燃氣輪機的傳遞路徑。沖擊波載荷以加速度的形式加載在箱裝體下部的減振元件上,并沿箱裝體內底板向機匣下部的每個基座傳遞,載荷到達燃氣輪機機匣與基座連接處位置后還需要以彎曲波的形式沿機匣長度方向傳遞并沿軸承向轉子傳遞,造成載荷在燃氣輪機內部分布的不均勻。為驗證上述載荷的傳遞規律,取燃氣輪機前基座附近處的剖面,給出在同一工況下不同時刻的應力云圖分布,如圖13所示。

圖12 燃氣輪機載荷傳遞路徑Fig.12 Gas turbine load transfer path

圖13 燃氣輪機典型剖面應力云圖Fig.13 Typical stress cloud picture of gas turbine section
從圖12和圖13中可以發現:爆炸沖擊波載荷作用在燃氣輪機上的范圍廣,作用強度大;在=1 ms 時,沖擊波首先作用在箱裝體底板,造成箱裝體底板的塑性應變,應力最大值達到345 MPa;隨著爆炸沖擊波載荷向燃氣輪機內部的傳遞,載荷作用在基座支承及基座;在=2.5 ms時基座處出現較大的應力,應力最大處主要出現在前部兩邊基座的中間轉接處位置,最大值為150 MPa,比箱裝體底板小,未達到材料的屈服極限,也未出現明顯的塑性變形,這是因為基座作為支承機匣和內部轉子的強力構件出現了應力集中的現象,在沖擊環境中的考核中需要著重注意此類強力構件;在3.75 ms時,機匣及轉子開始出現應力,隨著載荷的傳遞機匣和轉子的應力開始變大,相較于基座處應力值出現了明顯的衰減,但也未出現明顯的塑性變形,這是因為基座的載荷傳遞到燃氣輪機內部所引起的。
通過圖13中的應力云圖變化情況中可以發現,爆炸沖擊波載荷在燃氣輪機中的傳遞過程與理論結果一致,在傳遞過程中能量在衰減,應力響應也在衰減。總之,燃氣輪機結構的變化導致沖擊波載荷傳遞路徑和分布的不一致,為燃氣輪機在水下爆炸載荷沖擊作用下的不同響應提供了條件。
當使用多自由度系統對燃氣輪機進行力學分析時,可以將燃氣輪機簡化成如下3自由度力學模型,它包含基座、機匣、轉子三個質量點和剛度為、、,阻尼為、、的3自由度系統。在多數情況下,阻尼對系統的影響可忽略不計,當沖擊的激勵力與系統的固有頻率相近時,必須考慮阻尼的影響。

圖14 燃氣輪機力學模型Fig.14 Mechanical model of gas turbine
對如上3自由度系統進行力學分析并建立系統的運動微分方程,運用牛頓第二運動規律建立的方程為

(13)


(14)
在對設備進行抗沖擊考核時,只有設備基座會遭受到沖擊載荷,因此()=()=0 N。對系統的各參數及爆炸沖擊載荷確認后即可對系統的響應、能量傳遞及載荷傳遞規律進行分析。其中=5 130 kg,=170 kg,=4 860 kg,=518×10N/m,=17×10N/m,=49×10N/m,系統的阻尼比均為001。在瞬時沖擊激勵下得到、和的沖擊響應。
為驗證多自由度系統的有效性,首先假設=0 kg,=0 N/m,則系統為2自由度系統,且結構的初始值

(15)
則此方程為求解初始值已知的常微分方程。通常,用1階顯示的微分方程組來描述:

(16)

(16)式實際為運動方程的狀態空間表達形式,即(·)為
[,()]=A+
(17)
式中:

(18)
為單位矩陣。
本文用4階定步長Runge-Kutta法求解微分方程,先定義4個附加向量:

(19)
式中:為計算步長,一般取常數。
則下一個步長的狀態變量為

(20)
在已知初始值的情況下,可用(20)式一次次迭代,即可求解每時刻結構的響應。
在有限元軟件中提取機匣處的垂向加速度響應,將信號經過處理后作為2自由度系統的基礎激勵,其中=5 130 kg,=518×10N/m,=170 kg,=17×10N/m,阻尼比設為001,對設備的響應進行理論分析,用MATLAB軟件SimuLink中的仿真模塊對2自由度系統進行仿真,求解示意圖如圖15所示。圖15中,、分別為基座、機匣的加速度,、分別為基座、機匣的位移。

圖15 SIMULINK 2自由度系統求解框圖Fig.15 Block diagram of SIMULINK 2-DoF system solution
在Abaqus軟件中建立相同2自由度系統的彈簧質量點模型,沖擊載荷力施加在處,得到和處的位移響應。圖16為Abaqus 2自由度計算模型。

圖16 Abaqus 2自由度計算模型Fig.16 Abaqus 2-DoF calculation model
然后將理論計算結果與有限元仿真結果進行對比,其對比結果如圖17所示。

圖17 數值仿真結果與理論計算結果對比Fig.17 Comparison between numerical simulation results and theoretical calculation results
從圖17中可以看出,理論計算結果與有限元軟件仿真結果峰值相同,頻率一致,少量誤差可能來自有限元軟件求解器的原因,從而說明了燃氣輪機力學模型理論分析的正確性。
故可以使用MATLAB軟件中的SimuLink仿真功能對(13)式進行求解,用變步長4階、5級Runge-Kutta公式的ODE45算法,求解框圖如圖18所示,可得到、和的位移、速度及加速度響應。

圖18 SIMULINK 3自由度系統求解框圖Fig.18 Block diagram of SIMULINK 3-DoF system solution
采用與數值模型相同的信號處理方法對、和處得到的加速度信號進行分析,其分析結果如表3所示。

表3 燃氣輪機質量點加速度測點主頻率及 能量占比Table 3 Main frequency and energy ratio of gas turbine mass points as acceleration measurement points
在研究力學模型載荷傳遞的規律時,分別選取3個質量點為測點。對測得的數據進行經驗模態分解,記錄經驗模態分解后的主頻率以及能量占比情況如表3所示,基本呈現與有限元模型相似的規律。在載荷由基座向內部轉子傳遞的過程中,首先在基座部位激起的高頻響應占比較高為0.71,隨著載荷的向上傳遞,機匣的高頻響應為0.62,在轉子處的高頻響應占比為0.54,高頻響應呈現衰減的趨勢。分析結果與有限元模型分析相一致,表明可以用此力學模型來預先對燃氣輪機的響應進行預示。
本文對燃氣輪機沖擊載荷下的加速度響應進行EMD,以高頻信號的分布作為燃氣輪機響應的依據。得出以下主要結論:
1)載荷沿燃氣輪機底部向內部轉子傳遞,在傳遞的過程中,載荷隨傳遞距離的增加逐漸衰減,使得燃氣輪機各部件的沖擊響應也在逐漸減小。載荷在向燃氣輪機內部的傳遞過程中,高頻部分的占比逐漸減少,與燃氣輪機的應力響應變化相呼應。
2)載荷在傳遞的過程中,前部兩邊基座的中間轉接處位置處會產生較大的應力,尚未達到材料的屈服極限以及出現明顯的塑性變形,基座作為支承機匣和內部轉子的強力構件出現了應力集中的現象,在沖擊環境中的考核中需要著重注意此類強力構件。
3)通過將燃氣輪機簡化成3自由度力學預報模型,在驗證力學模型的有效性后,對力學模型的加速度響應進行EMD,得出與有限元模型相類似的規律。故燃氣輪機力學簡化模型能夠很好地預測燃氣輪機的載荷傳遞情況,可以成為快速工程化預報力學模型。