基于θ映射法的燃氣渦輪葉片高溫蠕變變形分析

全昌彪，廖明夫，米 棟，李 堅，劉 揚

(1.西北工業大學動力與能源學院，西安710072；2.中國航發湖南動力機械研究所，湖南株洲412002)

1 引言

高溫工作構件在持續應力作用下，其蠕變變形逐漸累積，最終可導致構件塑性變形過大或蠕變斷裂[1]。高溫合金的蠕變變形與溫度密切相關，工程上當溫度達到材料熔點溫度50%以上時蠕變就不能被忽略[2]。航空燃氣渦輪發動機的強度設計準則中，對熱端部件的蠕變變形、蠕變強度和蠕變/應力斷裂壽命均加以了規定[3]。近年來，航空燃氣渦輪發動機對高效率和高推重比/功重比的追求使得渦輪前溫度不斷提高[4]，從而導致其熱端部件(如燃氣渦輪葉片/渦輪盤)因長時蠕變而導致的結構失效問題日益突出。因此，準確進行蠕變變形計算對航空燃氣渦輪發動機熱端部件的強度設計極具工程意義。

通常，恒定溫度下構件受單軸恒定載荷發生的蠕變行為，可分為蠕變初始減速、蠕變恒速和蠕變加速三個階段[2]。其中，第三階段變形速率迅速上升導致最終失效，比較危險。傳統的蠕變方程，如時間硬化理論和應變硬化理論，適用于蠕變的第一和第二階段的分析計算。但工程實際中，由于溫度、載荷的多變性，有些材料并不能明顯體現出蠕變的三階段劃分，這就增加了蠕變計算分析的難度。要想準確描述蠕變的三階段行為，θ映射法無疑為一種較好的選擇。θ映射法是將蠕變變形描述為時間的指數函數[5]，該函數的兩個和式正好體現了蠕變減速階段和蠕變加速階段的疊加，而蠕變恒速階段則體現了減速階段和加速階段之間的平衡關系。Brown等[6]采用θ映射法進行了蠕變變形計算，指出該方法可完整描述蠕變三階段變形。Hayhurst等[7]應用θ映射法描述鍛造合金鋼1%Cr0.5%Mo0.25%V的蠕變行為，給出了擬合參數，其結果表明θ映射法在模擬高應力水平下的蠕變行為時比較接近實驗數據。Ibanez[8]采用θ映射法對定向凝固合金GTD-111的蠕變行為進行了模擬，其工作表明θ映射法在模擬定向凝固合金的蠕變行為時，其內插性優良，而外推能力相對較差。

從本質上講，θ映射法對材料蠕變力學性能的表征，是一種基于宏觀連續介質力學的唯象方法，其數學模型沒有考慮材料微觀層次的變形機制，因此計算效率相對較高。然而，為實現三維結構的蠕變分析，還需要將宏觀唯象的單軸本構拓展為多軸形式。鑒于此，本文按照Prandt-Reuss[9]塑性流動法則將θ映射法拓展成適合三維有限元分析的多軸形式，并進一步將其編制成UMAT用戶子程序植入到ABAQUS有限元軟件，對某型渦軸發動機渦輪葉片材料DZ408的縱向蠕變行為進行模擬計算。同時，對葉片在該型發動機60 h整機持久試車過程的蠕變行為進行近似模擬，并就計算所得葉尖位移量與3次整機60 h試車后的實測伸長量進行對比。

2 蠕變模型

2.1 單軸蠕變方程

對于各向同性材料，蠕變應變εc通常可表示為構件所受應力σ、溫度T和時間t的函數，即：

如果給定試驗條件(σ，T)時，蠕變應變只是時間t的函數：

Evans和Wilshire認為蠕變過程可描述如下[5]：

式中：θi(i=1～4)是與材料、溫度及應力有關的常數。對于特定材料，θi可表示為應力和溫度的函數：

式中：ai、bi、ci、di為材料特性相關的參數，可根據蠕變試驗曲線進行優化獲取。

2.2 多軸蠕變方程

由于θ映射模型的基本方程只是單軸的形式，而復雜結構一般都處于多軸應力狀態，這就要求將單軸形式的本構模型向多軸形式拓展。多軸蠕變理論主要特點是考慮了變形與時間的關系，且隨著時間的變化其蠕變規律會發生非線性變化。一般認為，多軸應力狀態下蠕變模型滿足以下基本假設[10]：①多軸應力狀態下的蠕變公式必須可以退化成正確的單軸蠕變公式；②蠕變變形前后模型體積不變；③蠕變計算方程不受靜水壓力的影響；④各向同性材料的主應力和主應變的主方向一致。根據基本假設，塑性應力應變理論可類似在蠕變分析中推廣應用。小變形情況下，蠕變應變率張量ε?c服從Prandt-Reuss塑性流動[11]：

式中：λ為塑性乘子，Sij為應力偏張量，δij為Kronecker Delta函數，σkk為3個主應力分量之和。

由塑性理論分析可知，為確定式(5)中λ的值，根據單一曲線假設[12](硬化條件)，多軸的等效蠕變應變率應和單軸的蠕變應變率具備相同的形式，則λ值由下式決定。

對式(3)求導，可得到蠕變曲線的斜率公式：

式中：為等效應力，其表達式如下：

將式(7)、式(8)代入式(5)，則可得到多軸蠕變應變率的表達式：

3 蠕變模型參數的確定

DZ408材料在同一溫度下進行了多種應力水平下的蠕變試驗，表1匯總了其蠕變試驗條件。本文根據其蠕變試驗數據，采用最小二乘法按照式(3)和式(4)擬合出其中的蠕變參數，結果見表2和表3。

表1 DZ408材料蠕變曲線試驗條件Table 1 Creep test conditions of DZ408

表2 θ參數擬合結果Table 2 Fitting results of parameterθ

表3 a、b、c、d參數擬合結果Table 3 Fitting results of parametera,b,c,d

為獲取蠕變曲線的擬合精度，驗證模型參數的準確性，將2.2節多軸蠕變方程編制成UMAT用戶程序嵌入ABAQUS有限元軟件，采用體積胞元模型對DZ408材料縱向蠕變行為進行了數值模擬。圖1給出了有限元計算曲線與試驗曲線的對比，圖中孤立點為試驗曲線，實線為計算曲線?？煽闯?，計算曲線與試驗曲線較為吻合，驗證了UMAT用戶程序的正確性和蠕變模型方程的有效性。同時，由圖1(a)可知，應力水平為540 MPa的計算曲線落在應力水平為530 MPa和550 MPa的兩條試驗曲線中間，說明蠕變模型參數能夠適應內插。此外，從圖1中蠕變曲線的特征還可看出：在較低的溫度和相對較高的應力水平(圖1(a))下，材料的蠕變曲線基本上沒出現第三階段；而在較高的溫度和相對較低的應力水平(圖1(c)和圖1(d))下，材料在較短時間內就出現了第三階段。由此表明，蠕變變形首先依賴于溫度，其次是應力[13]。

4 算例分析

對某型燃氣渦輪葉片60 h整機持久試車過程的蠕變行為進行了近似模擬。圖2為簡化處理后的葉片有限元網格模型。網格采用四節點四面體單元(C3D4）劃分，共計148 582個單元，40 188個節點。

計算主要考慮了試車過程中不同載荷譜下的離心載荷和溫度載荷。為模擬載荷的瞬態變化，在分析文件中定義了如表4所示的離心載荷和發動機動力渦輪前溫度(T45)的歷程，渦輪葉片的溫度分布如圖3所示。

表4 各分析步計算載荷Table 4 Calculation load of each analytical procedure

表5示出了給定計算載荷下，葉片在60 h蠕變過程中葉尖徑向位移計算結果。卸載后葉尖殘余變形量為0.086 94 mm。

表5 各分析步位移計算結果Table 5 Displacement calculation results of each procedure

從葉身根部至葉尖依次選取4個代表節點N913、N1978、N1833、N476，提取出其徑向位移隨時間的變化關系，見圖4。從圖中可看到，葉身的蠕變位移-時間曲線有明顯的蠕變三階段特征，從初始加載時的蠕變減速階段，很快進入蠕變恒速階段，接近60 h時又逐漸表現出蠕變加速階段的特征。由此可知，渦輪葉片在較高的溫度和離心載荷等的作用下，其變形與載荷不再是簡單的一一對應關系，在載荷作用下還會隨著時間的推移而逐漸增加。

葉片在整個蠕變過程中，葉身應力重新分布。圖5給出了葉身中部位置N1341和N696節點當量應力隨蠕變時間推移而發生的應力松弛現象。由圖可看出，應力松弛現象使得葉身高應力區的應力有所降低，低應力區的應力有所增加，整個葉片的應力分布趨向均勻化，這有益于葉片的使用壽命。

5 計算變形與實測伸長量對比

卸載后葉尖的最大徑向位移為0.086 94 mm。發動機在3次60 h持久試車后，渦輪葉片平均伸長量的實測值分別為0.058、0.083、0.080 mm，平均值為0.074 mm。計算結果與實測結果平均值相對誤差約18.02%，其中有兩次不到10.00%，對比結果見表6。

表6 葉片伸長量計算值和實測值對比Table 6 Comparison between blade elongation calculation results and test results

6 結論

(1)推導出的多軸形式的θ映射法蠕變本構模型，能較好地模擬高溫合金DZ408材料的縱向蠕變特性。

(2)渦輪葉片發生蠕變后，其應力將重新分布，葉身的變形也會隨著工作時間的推進而逐漸增加，在葉片蠕變變形設計和葉尖間隙設計時應予以考慮。

(3)采用多軸形式的θ映射法蠕變本構模型對某渦輪葉片進行三維蠕變分析得到的60 h葉片徑向變形量與3次整機60 h持久試車后實測伸長量的平均值較為接近，平均值相對誤差約18.02%，其中有兩次誤差低于10.00%。

(4)受材料僅有縱向蠕變試驗曲線的限制，文中未考慮DZ408材料蠕變行為本身的各向異性，而是將其當作各向同性處理，這也是計算精度影響因素之一。此外，計算的葉片伸長量比3次實測平均值都要大，表明θ映射法在預測葉片的蠕變變形用于指導設計偏安全，可通過后續研究中持續積累試驗數據與計算模擬逐步形成適用于工程應用的經驗修正系數。