流固耦合作用下航空發動機燃燒室熱疲勞研究

張俊紅, 戴胡偉, 魯鑫, 王杰, 馬梁, 孫宇博

(1.天津大學內燃機燃燒學國家重點實驗室, 300354, 天津; 2.天津大學仁愛學院機械工程系, 301636, 天津;3.中國民航大學天津市民用航空器適航與維修重點實驗室, 300300, 天津)

航空發動機燃燒室是航空發動機的關鍵熱端部件,其強度和可靠性對發動機運行安全性有直接的影響。目前針對航空發動機燃燒室的研究主要分為流場研究和疲勞壽命研究。

在流場研究方面,由于民航運行安全規定,難以監測航空發動機工作時燃燒室內流場分布的真實情況。對燃燒室流場的研究大多以仿真為主。王寶官等采用二維紊流模型對有氣膜冷卻的燃燒室內流場及壁溫進行數值分析,得到速度場、壁面溫度、輻射通量隨進口參數的變化規律[1]。雷雨冰等采用三維紊流模型通過多區耦合法對某環形燃燒室化學反應流場進行數值模擬[2],較全面地反映燃燒室的氣流流動、換熱、燃燒等現象;劉常春等采用SST(shear-stress transport)湍流模型研究了孔陣排列和偏轉角對火焰筒冷卻效果的影響[3];Fureby采用大渦模擬的方法對環形燃燒室三維兩相流場進行模擬,獲得燃燒室內流場分布及溫度、壓力的波動情況[4]。然而,諸多學者的研究重點大多集中于燃燒室的流動特性,對燃燒室內流體域與結構之間的耦合傳熱特性及其對燃燒室疲勞可靠性研究較少。

在疲勞壽命研究方面,Barrett等通過實驗研究了高溫環境下哈氏合金低周疲勞性能[5];Meyer等通過對比試驗研究了焊縫對哈氏合金疲勞特性的影響[6];易慧通過對燃燒室施加簡化的溫度載荷獲得燃燒室應力應變,并對燃燒室進行壽命預測[7]。實驗研究及對燃燒室基體施加的簡化載的荷仿真研究不能準確地反應燃燒室復雜多變的工作環境,壽命預測結果準確性較低。目前,對燃燒室工作過程進行CFD仿真準確獲得燃燒室溫度及應變載荷,并對燃燒室進行疲勞可靠性研究很少。

本文建立了某航空發動機燃燒室湍流燃燒流固耦合模型,對慢車、起飛、爬升、巡航、下降5種典型工況下燃燒室流場進行了數值模擬,獲得了燃燒室的溫度分布。在驗證流固耦合計算準確性的基礎上,基于流場計算結果通過非線性靜力學分析獲得了燃燒室基體的應力和應變分布。

1 基本物理方程

1.1 流固耦合計算方程

本文采用簡化的聯合概率密度函數(PDF)模型模擬發動機燃燒過程,非預混燃燒PDF模型的平均混合分數方程及混合分數方差的守恒方程分別為

(1)

Cgμt(

(2)

(3)

考慮流體與結構之間的相互作用,在作用處邊界需滿足流體與固體溫度、熱流量、位移、應力相等,滿足以下方程

(4)

式中:T為溫度;q為熱流量;d為位移;τ為應力;下標f、s分別表示流體、固體。

1.2 非線性靜力學分析方程

假定構件溫度變化為ΔT,則應力、應變和溫度的關系式可表示為

在彈性區域

σ=De(ε-εT)

(5)

在塑性區域

σ=Dep(ε-εT)

(6)

εT=αT

(7)

式中:De為彈性矩陣;Dep為相應的彈塑性矩陣;εT為熱應變增量;α為熱膨脹列陣。

根據虛位移原理,可建立節點載荷有限元方程的增量表達式

在彈性區域中的單元

(8)

在塑性區域中的單元

(9)

總體載荷

R=∑Re

(10)

基本平衡方程式可寫為

Kδ=R

(11)

式中:K為構件剛度矩陣;δ為結構節點位移列向量。

1.3 熱疲勞分析方程

在結構疲勞壽命計算中,Manson-Coffin公式和線性損傷理論應用廣泛。Manson在獨立研究熱疲勞問題的過程中提出一種以塑性應變幅為參量的疲勞壽命描述法[8],即

(12)

線性累積損傷理論由Miner提出[9],零部件的損傷變量為

(13)

式中:ni為第i段載荷下的循環數;Ni為第i段載荷下結構失效的循環數;k為載荷數。當D為1時,可認為零部件失效,發生疲勞破壞[10]。

2 幾何建模和力學性能實驗

2.1 幾何建模

本文以某航空發動機環形燃燒室為研究對象。該環形燃燒室火焰筒前端沿周向均布20個旋流器,每個旋流器中心有一個雙油路離心噴嘴;火焰筒沿周向均布80個主燃孔、120個摻混孔;火焰筒內、外壁分別沿軸向布置4、5道冷卻氣膜,共10 540個氣膜冷卻孔。

為得到與燃燒室實體吻合度較高的計算模型,采用HEXAGON公司的ROMER絕對關節臂外接激光型非接觸式三坐標掃描儀獲取環形燃燒室截面點云數據及整體點云數據。通過對后期坐標數據的逆向處理實現燃燒室實體模型的建立,燃燒室點云數據及三維模型如圖1所示。

(a)截面點云數據 (b)整體點云數據 (c)燃燒室三維模型圖1 燃燒室點云數據及三維模型

2.2 哈氏合金拉伸試驗

燃燒室基體材料為哈氏合金X,為獲得哈氏合金X準確力學性能參數,進行拉伸試驗。依據GB/T228.1—2010[11]設計矩形截面拉伸試樣,尺寸和試驗前、后試樣如圖2所示。試驗在INSTRON萬能力學試驗機上進行,選用4個平行試樣依據GB/T228.1—2010規定的方法,通過控制試驗機橫梁位移速率控制試件的應變速率,橫梁位移速率為1 mm/min,通過0.5級全自動軸向引伸計測量試樣軸向變形,實驗數據采樣頻率為10 Hz。

(a)拉伸試驗矩形截面試樣尺寸

(b)試驗前試樣

(c)試驗后試樣圖2 哈氏合金X力學性能試驗試樣

實驗測得哈氏合金X彈性模量為278 MPa,屈服強度為463 MPa,抗拉強度為980 MPa。繪制哈氏合金X拉伸試驗工程應力-應變曲線,依據GB/T228.1—2010用多項式回歸方法對工程應力-應變曲線進行光滑處理,哈氏合金X工程應力-應變曲線如圖3所示。

圖3 哈氏合金X工程應力-應變曲線

在拉伸試驗中,由于試樣的橫截面積減小,其真實應力比工程應力大,試件的真實應力為

σa=σ(1+ε)

(14)

(15)

式中:ε為試件應變。

3 流固耦合模擬

3.1 數學模型

利用FLUENT系統對典型工況下燃燒室內的湍流燃燒進行模擬,通過編寫UDF(user define function)程序實現迭代計算過程中流體域與固體域之間的數據交換,使得流體域與固體域耦合面上的位移、熱流量、溫度、應力等相等。選用模擬旋流能力強且收斂性較好的Realizablek-ε雙方程模型模擬燃燒室內的湍流流動[12],選用P1輻射模型模擬燃燒室內的輻射傳熱,選用壓力霧化噴嘴模型模擬燃油噴射過程,采用隨機軌道模型進行兩相流計算,用拉格朗日法跟蹤離散液滴在流場中的運動和運輸,用歐拉法描述氣相守恒,通過在氣相守恒方程中加入相應的源項來考慮液滴對氣相的影響;選用PDF非預混燃燒模型模擬湍流擴散燃燒,通過在能量守恒中加入相應的源項來考慮燃燒對流場的影響,采用SIMPLE算法和二階迎風格式進行求解。

3.2 網格劃分

燃燒室具有以2個旋流器頭部為周期的周向周期性結構特征。為減少計算量,選取2個旋流器頭部,即1/10扇形段作為計算域,如圖4a所示。燃燒室流體計算域形狀復雜,為更準確地描述流體域形狀,采用變形能力強的四面體網格對流體域進行網格劃分,基本網格尺寸大小為5 mm。氣膜冷卻孔對燃燒室壁面溫度分布有極大的影響,其孔徑僅有1.5 mm,為保證計算準確性,對氣膜冷卻孔部分網格進行局部加密,此處網格大小0.3 mm。氣膜冷卻孔處網格尺寸與其他區域網格尺寸差別較大,為防止尺寸差別過大導致網格畸形引起CFD計算不準確,設定網格尺寸增長率不超過20%。為穩定氣流、改善流場計算的收斂性,添加進、出口流道,設置進、出口流道長度為30 mm,以減小進出口邊界壓力反射對流場計算域的影響。

燃燒室基體為薄壁件,壁厚小于2 mm,其形狀對流場影響較大,為更好地反應燃燒室形狀對流場的影響,采用基本尺寸為1 mm具有二階精度四面體單元對燃燒室基體進行網格劃分,對燃燒室基體氣膜冷卻孔、摻混孔孔邊等位置網格進行局部加密,網格尺寸0.3 mm。燃燒室基體表面噴涂有一層厚約0.2 mm的熱障涂層,熱障涂層起隔熱作用,可降低燃燒室基體的工作溫度[13],本文計算考慮熱障涂層的影響,采用基礎尺寸0.2 mm的四面體網格對其進行網格劃分。對流體域及固體域進行網格劃分,獲得2 339 390個流體網格單元及2 305 363個固體網格單元,流固耦合計算模型如圖4所示。

(a)1/10扇形段計算域

(b)流體域網格

(c)固體域網格圖4 流固耦合計算模型

3.3 結果分析

本文基于發動機試車數據,將該發動機完整的工作循環簡化為慢車、最大、爬升、巡航、下降5個典型工況,各工況進氣壓力、進氣溫度、噴油量如表1所示。流場入口采用壓力入口邊界條件,出口截面采用自由流邊界條件,壁面采用無滑移邊界條件,通過壁面函數法確定近壁處湍動能k及湍動能耗散率ε,燃燒室兩側面為旋轉周期性邊界條件,旋轉角度為36°。

表1 典型工況下發動機運行參數

在熱障涂層服役過程中,環境溫度高于熱障涂層相變溫度時,熱障涂層會發生相變和燒結[14]隨著熱障涂層服役時間的增加,不同相變程度的位置宏觀樣貌會有所不同。為驗證流固耦合仿真結果的準確性,將熱障涂層外表面溫度場與實際服役燃燒室熱障涂層宏觀樣貌進行對比。由于各個工況下熱障涂層表面溫度值大小有一定差異,但其分布規律相似,選用巡航狀態下熱障涂層溫度場與實際服役燃燒室熱障涂層宏觀樣貌進行對比,如圖5所示。

(a)溫度場 (b)燃燒室圖5 熱障涂層溫度場與實際服役燃燒室宏觀樣貌對比

不同工況下燃燒室基體表面溫度分布如圖6所示,慢車狀態、最大狀態、爬升狀態、巡航狀態、下降狀態燃燒室基體表面最高溫度分別為1 081.60、1 387.39、1 329.93、1 243.26、1 147.82 K,由于燃燒室內環壁摻混孔下游區域靠近輻射高溫區,該位置整體溫度較高,最高溫度點也位于此區域;在摻混孔正下方小部分區域,從摻混孔處進入燃燒室內部的氣流對此區域有一定的冷卻作用,使得該區域溫度較低;在摻混孔下游區域,溫度分布不均勻,高溫區與低溫區最大溫差可達237 K(最大狀態下),這易使得該區域熱應力較大導致疲勞破壞,燃燒室基體實際開裂位置也多位于此處。

(a)慢車狀態 (b)最大狀態

(c)爬升狀態 (d)巡航狀態

(e)下降狀態圖6 不同工況下燃燒室基體表面溫度分布

4 非線性靜力學分析

4.1 載荷及邊界條件

燃燒室承受冷熱氣流沖擊受熱不均勻,各處變形不一致相互約束產生熱應力[15]。將流固耦合計算得到的溫度作為載荷;將哈氏合金X力學性能拉伸試驗結果作為燃燒室基體材料屬性;按照燃燒室實際裝配情況,約束燃燒室基體底面軸向自由度,約束周期性對稱面周向自由度;通過有限單元法對燃燒室基體基本平衡方程進行求解,獲得各節點的應變分布。

圖7 主燃孔及摻混孔編號

4.2 結果分析

為便于后續分析,對燃燒室基體壁面部分主燃孔摻混孔進行編號,如圖7所示。各工況下燃燒室基體塑性變形如圖8所示。由圖8可知,各工況下燃燒室基體塑性應變主要發生在主燃孔、摻混孔及氣膜冷卻孔下游區域,隨著發動機負荷的增大,燃燒室基體熱負荷增加,燃燒室整體塑性變形逐漸增大,燃燒室各工況下最大應變點位置及最大應變如表2所示。

表2 各工況下燃燒室最大應變點位置及最大應變

(a)慢車狀態 (b)最大狀態

(c)爬升狀態 (d)巡航狀態

(e)下降狀態圖8 燃燒室塑性應變

應變較大位置與實際裂紋對比如圖9所示,各應變較大位置對應實際服役燃燒室基體位置均出現了不同程度的裂紋;其中,最大狀態下及下降狀態下,燃燒室基體最大應變位置與實際服役燃燒室基體裂紋位置基本吻合;慢車狀態及巡航狀態下,燃燒室基體應變最大位置對應實際服役燃燒室基體開裂位置偏差約為1.3 mm、1.5 mm,這是由于慢車狀態及巡航狀態下,燃燒室基體最大應變位置出現在1號主燃孔及3號主燃孔處,該區域靠近燃燒室噴油嘴,燃油濃度較高,燃燒過程中易形成碳煙,隨著航空發動機服役時間的增加,碳煙附著在燃燒室壁面使得該區域局部傳熱特性發生改變,使得最大應變位置與實際裂紋位置有較小的偏差。總體而言,模擬結果的最大應變位置與燃燒室實際失效位置對應良好,表明有限元計算結果可靠,可作為下一步計算的依據。

(a)慢車狀態下1號孔應變 (b)1號孔裂紋

(c)最大狀態下2號孔應變 (d)2號孔裂紋

(e)巡航狀態下3號孔應變 (f)3號孔裂紋

(g)下降狀態下4號孔應變 (h)4號孔裂紋圖9 應變較大位置與實際裂紋對比

5 疲勞壽命預測

5.1 危險點的選取及載荷譜的確定

根據應變分析結果及實際服役航空發動機燃燒室基體裂紋分布情況,對1號、3號主燃孔及2號、4號摻混孔4個易失效的位置進行壽命預測。

本文研究對象為某民用飛機上使用的航空發動機,從其飛行記錄中摘取飛機發動機19個起落循環的數據,獲得其典型工作循環,如圖10所示[16]。由于機動狀態負荷較小,加速、反推力、減速過程時間很短,這幾個階段對疲勞壽命影響較小,因此忽略這幾個階段的影響,將發動機典型工作循環簡化為“慢車—起飛—爬升—巡航—下降—慢車”。

圖10 某航空發動機典型工作循環

5.2 危險點疲勞壽命預測

該機型的設計維修手冊中要求對燃燒室可靠度的為99.5%,采用文獻[17]所述的方法,設定樣件存活率99.5%,置信區間95%,對Hong等的試驗數據[18]進行擬合,得到873、1 033、1 143 K環境下哈氏合金X的Manson-Coffin公式,分別為

εp1=0.072 6(2Nf)-0.340 87

(16)

εp2=0.051 9(2Nf)-0.351 59

(17)

εp=0.275(2Nf)-0.590 44

(18)

根據流場仿真結果,危險點不同工況下工作溫度范圍在856～1 191 K之間,分別將不同工況下危險點應變代入上式,求得危險點在873、1033、1 143 K下疲勞壽命,通過Lagrange插值法獲得危險點在實際工作溫度下疲勞壽命。通過線性累積損傷理論計算得到在典型工作循環下各危險點位置的損傷累積值,當D達到1時,可認為零部件失效,發生疲勞破壞。計算結果表明1號、3號主燃孔及2號、4號摻混孔壽命循環數為13 754、9 381、7 126、11 693。

6 結 論

本文基于流固耦合方法對某航空發動機燃燒室典型工況下流場進行了仿真,獲得了不同工況下燃燒室基體溫度,通過非線性靜力學分析得到了不同工況下燃燒室基體塑性應變分布,對危險點進行了壽命預測,得到以下結論。

(1)通過考慮燃燒室基體及熱障涂層的典型工況燃燒室CFD仿真發現,在慢車狀態下,燃燒室基體整體溫度較低,隨著發動機負荷上升,燃燒室基體整體溫度逐漸升高。由于靠近燃燒高溫區域,燃燒室摻混孔下游大部分區域溫度較高;由于摻混氣流對基體的冷卻作用使得摻混孔下游小部分區域燃燒室基體溫度較低,摻混孔下游區域溫度分布不均勻。

(2)通過非線性有限元仿真發現,不同工況下燃燒室基體塑性應變分布規律相似,隨著發動機負荷的上升,燃燒室基體塑性應變值逐漸增大;塑性應變主要出現于主燃孔及摻混孔下游區域,各工況最大應變位置對應實際服役燃燒室基體位置均有不同程度的裂紋,仿真最大應變位置與實際燃燒室位置偏差小于1.5 mm。

(3)考慮溫度的影響,通過Manson-Coffin公式及線性累積損傷理論計算得到在典型工作循環下危險點疲勞壽命,壽命最低點為2號摻混孔處,最低壽命循環數為7 126。

參考文獻:

[1] 王寶官, 李永康, 胡正義. 帶氣膜冷卻的火焰筒壁溫的數值分析 [J]. 航空學報, 1995, 16(4): 415-421.

WANG Baoguan, LI Yongkang, HU Zhengyi. Numerical analysis of the wall-temperature with film-cooling in the combustion chamber [J]. Acta Aeronautica ET Astronautica Sinica, 1995, 16(4): 415-421.

[2] 雷雨冰, 胡好生, 趙堅行. 環形回流燃燒室兩相反應流場的數值研究 [J]. 航空動力學報, 2001, 16(4): 350-354.

LEI Yubing, HU Haosheng, ZHAO Jianxing. Numerical study of two-phase reacting flow in annual reversed combustor [J]. Journal of Aerospace Power, 2001, 16(4): 350-354.

[3] 劉常春, 吉洪湖, 楊芳芳, 等. 孔陣排列和偏轉角對多斜孔壁火焰筒冷卻效果的影響研究 [J]. 推進技術, 2013, 34(10): 1369-1375.

LIU Changchun, JI Honghu, YANG Fangfang, et al. Numerical study on cooling effects of hole arrangement and deflected angle of inclined multi-hole on annular flame tube [J]. Journal of Propulsion Technology, 2013, 34(10): 1369-1375.

[4] FUREBY C. LES of a multi-burner annular gas turbine combustor [J]. Flow, Turbulence and Combustion, 2010, 84(3): 543-564.

[5] BARRETT P R, AHMED R, MENON M, et al. Isothermal low-cycle fatigue and fatigue-creep of Haynes 230 [J]. International Journal of Solids & Structures, 2016, 88/89: 146-164.

[6] MEYER-OLBERSLEBEN F, KASIK N, ILSCHNER B, et al. The thermal fatigue behavior of the combustor alloys In 617 and HAYNES 230 before and after welding [J]. Metallurgical and Materials Transactions: A, 1999, 30(4): 981-989.

[7] 易慧. 環形燃燒室火焰筒強度壽命技術研究 [D]. 南京: 南京航空航天大學, 2008: 5-6.

[8] MANSON S S. Fatigue: A complex subject some simple approximations [J]. Experimental Mechanics, 1965, 5(7): 193-226.

[9] MINER M A. Cumulative damage in fatigue [J]. Journal of Applied Mechanics, 1945, 12(3): 159-164.

[10] 付曦, 張俊紅, 寇海軍, 等. 復雜載荷下軸流壓氣機葉片疲勞損傷數值研究 [J]. 西安交通大學學報, 2017, 51(5): 149-155.

FU Xi, ZHANG Junhong, KOU Haijun, et al. Numerical study on the fatigue damage of compressor blade under complex loads [J]. Journal of Xi’an Jiaotong University, 2017, 51(5): 149-155.

[11] 全國鋼標準化技術委員會. 金屬材料拉伸試驗第1部分: 室溫試驗方法GB/T228.1—2010 [S]. 北京: 中國標準出版社, 2010.

[12] JONES W P, LAUNDER B E. The prediction of laminarization with a two-equation model of turbulence [J]. International Journal of Heat & Mass Transfer, 1972, 15(2): 301-314.

[13] PADTURE N P, GELL M, JORDAN E H. TBCs for gas-turbine engine applications [J]. Science, 2002, 296(5566): 280-284.

[14] 劉懷菲. 二元稀土氧化物復合穩定氧化鋯熱障涂層材料的制備及性能研究 [D]. 長沙: 中南大學, 2011: 47-50.

[15] 嚴宗達, 王洪禮. 熱應力 [M]. 北京: 高等教育出版社, 1993: 76-82.

[16] 付娜. 某航空發動機渦輪盤和葉片的強度分析與壽命計算 [D]. 西安: 西北工業大學, 2006: 12-13.

[17] 傅惠民, 劉成瑞.ε-N曲線和P-ε-N曲線整體推斷方法 [J]. 航空動力學報, 2006, 21(6): 957-961.

FU Huimin, LIU Chengrui. Integral influence method forε-NandP-ε-Ncurves [J]. Journal of Aerospace Power, 2006, 21(6): 957-961.

[18] HONG H U, KIM I S, CHOI B G, et al. Effects of temperature and strain range on fatigue cracking behavior in hastelloy X [J]. Materials Letters, 2008, 62(28): 4351-4353.

[本刊相關文獻鏈接]

季家東,葛培琪,畢文波,等.采用不同管排組合的換熱器彈性管束殼程流體誘導振動響應.2018,52(3):69-75.[doi:10.7652/xjtuxb201803010]

仲繼澤,謝志強,沈渡,等.基于虛功原理的流固耦合面力和位移傳遞方法.2018,52(3):160-167.[doi:10.7652/xjtuxb 201803022]

文鍵,李科,劉育策,等.利用流固耦合分析的板翅式換熱器鋸齒型翅片多目標優化.2018,52(2):130-135.[doi:10.7652/xjtuxb201802020]

仲繼澤,謝志強,沈渡,等.基于空間分布彈性模量的快速動網格方法.2018,52(2):136-139.[doi:10.7652/xjtuxb2018 02021]

仲繼澤,徐自力.流固單向耦合的能量法及機翼顫振預測.2017,51(1):109-114.[doi:10.7652/xjtuxb201701017]

范增華,榮偉彬,王樂鋒,等.壓電微噴輔助液滴的多物理場耦合與實驗.2016,50(11):56-61.[doi:10.7652/xjtuxb2016 11009]

仲繼澤,徐自力,陶磊.基于虛擬彈性體的快速動網格方法.2016,50(10):132-138.[doi:10.7652/xjtuxb201610020]

姜濤,黃偉,王安麟.多路閥閥芯節流槽拓撲結構組合的神經網絡模型.2016,50(6):36-41.[doi:10.7652/xjtuxb2016 06006]

高炎,晏鑫,李軍.燃氣透平葉片尾緣開縫結構冷卻性能的數值研究.2016,50(3):29-37.[doi:10.7652/xjtuxb201603005]

郭濤,管志成,孫光普,等.調頻振子-液體聯合水平減振的流固耦合機理研究.2016,50(1):28-33.[doi:10.7652/xjtuxb 201601005]

宋明毅,吳偉烽,李直.汽車空調壓縮機氣閥運動規律模擬.2015,49(12):144-150.[doi:10.7652/xjtuxb201512023]

季家東,葛培琪,畢文波.換熱器內彈性管束流體組合誘導振動響應的數值分析.2015,49(9):24-29.[doi:10.7652/xjtuxb201509005]

曹文瑾,孫中國,席光.系泊型浮體運動的無網格法數值模型.2015,49(3):62-66.[doi:10.7652/xjtuxb201503011]