基于梁理論的渦輪冷卻葉片蠕變計算

李錦紅，張 勇，張如剛，章 勝

(中國航發湖南動力機械研究所，湖南株洲412002)

1 引言

隨著航空發動機技術的發展，發動機的工作溫度越來越高。特別是其中的燃氣渦輪轉子葉片[1]，在工作中承受很大的熱負荷和機械負荷，蠕變損傷成為渦輪葉片的一種主要失效形式，短時強度的設計理念已無法滿足渦輪轉子葉片的使用要求。因此，對渦輪轉子葉片進行準確的蠕變分析，對其設計具有重要意義。

目前，國內外比較廣泛應用的是利用損傷原理，將材料不同工況下的蠕變曲線轉化為不同溫度下的等效等時應力應變曲線進行三維彈塑性有限元計算[2-4]。該方法可利用有限元分析軟件較為簡便地模擬冷卻葉片的復雜結構，但僅能計算蠕變行為的前兩個階段，且每個溫度下只能有一條等效等時應力應變曲線，這與材料的實際蠕變行為不符。部分學者引用多軸時間硬化理論[5]，依據Prandt-Reuss塑性流動法則將單軸的θ模型方程擴展為多軸形式，并將其編寫成UMAT用戶子程序，嵌入到Abaqus有限元軟件，對渦輪冷卻葉片進行三維蠕變計算。該方法引用θ參數可較準確地描述材料的蠕變行為，但未考慮蠕變對應力的影響，且其計算極其繁瑣。

本文以梁理論為基礎，將三維葉片簡化為二維截面的疊加，采用累積損傷理論及材料的拉遜-米勒曲線方程計算各截面在不同工況下的等效應力，然后結合給定工況的持久時間和溫度，利用實測材料的蠕變曲線進行線性插值計算截面的平均蠕變應變，最后對各截面及各工況的蠕變伸長量進行疊加。該方法直接利用材料的蠕變曲線插值，不僅可計算蠕變行為三個階段[6]的蠕變應變，而且還可考慮同溫度下不同應力對蠕變行為的影響。

2 基本原理與方法

高溫下金屬材料的組織結構對材料的蠕變特性影響很敏感，蠕變力學特性比較復雜，同一溫度下存在多條隨時間變化的蠕變應變曲線(圖1)，當溫度一定時蠕變應變εc為應力σ和時間t的函數：

圖1 同一溫度下材料的蠕變曲線Fig.1 Creep curve at the same temperature

一般認為，多維應力狀態下的蠕變公式必須能退化成正確的一維蠕變計算公式[1]，因而需要將圖1中的多維蠕變曲線轉化為一維蠕變曲線。

高溫下，在恒定載荷作用時間tΣ內，材料的應力并非一定值，為此引入等效應力σequiv的概念。恒定載荷在作用時間內的等效應力表示為：

式中：m為隨溫度變化的持久強度指數。同一溫度下，對金屬材料有常數，其中σr為持久強度(MPa)，t為σr對應的持久時間(h)。

σequiv在tΣ時段是恒定的，其引起的損傷相當于在相同時間內隨時間變化的實際應力造成的損傷。當溫度一定時，蠕變應變與應力及時間的關系可轉化為：

式中：q為給定試驗條件(溫度T，tΣ)對材料造成的損傷，其表達式為

由拉遜-米勒曲線可知，σequiv與其在某一溫度下對應的持久時間t的關系為：

式中：a1、a2、a3、a4可由材料的拉遜-米勒曲線(圖2)方程得到。

圖2 拉遜-米勒曲線Fig.2 The Larson-miller curve

3 等效應力計算

按照平截面的假設，葉片橫截面任意點的縱向變形，是由彈性變形εe、塑性變形εp、εc和因葉片截面不均勻溫度場而引起的熱變形一起疊加的結果，即截面的應變為：

式中：α為線膨脹系數，ΔT為溫升。

由梁理論可知，在平行于葉片截面的x-y平面上，葉片截面的應變為：

式中：F為葉片截面上的法向力合力，Mx為相對于y軸的彎矩，My為相對于x軸的彎矩，S為葉片截面面積。

考慮應力松弛后，應力為：

蠕變應變εc隨時間的變化規律可用蠕變變化率描述。此處引入Lemaitre應變等價原理，單軸損傷蠕變變化率為：

由蠕變試驗可看到蠕變模型在變形前后體積不變[1]。本文假設多軸蠕變與單軸蠕變規律相當，則多軸狀態下的蠕變變化率為：

將載荷作用時間離散為若干時間步，每一時間步長為 Δti(i=0,…,k)。

(1) 當i=0時，εic=0，根據式(7)和式(8)可計算該時間步下的等效應力σ(Δti)。(2) 當i=i+1時，由式(5)可知σ(Δti)在給定溫度條件下的持久時間為t，造成的損傷q=Δti+1/t，根據式(10)，蠕變應變，根據式(7)和式(8)可計算該時間步下的等效應力σ(Δti+1)。

通過以上計算可得到每個時間步下的等效應力，根據式(2)可得到載荷作用時間段的等效應力。

4 計算實例

以某航空發動機燃氣渦輪工作葉片為例，用上述方法計算渦輪工作葉片持久試車試驗后的蠕變變形。根據飛行包線，該持久試車可劃分為5個工況。本文忽略加卸載過程的加減速，認為每個工況內的載荷恒定。

4.1 計算模型

該渦輪工作葉片為鑄件，各工況下的持久試車時間見表1。

表1 持久試車時間Table 1 Time of endurance test

將渦輪葉片沿徑向平均分為n段，第n段對第n截面提供離心拉力及彎矩。其中第1截面為葉根截面，第n+1截面為葉尖截面，如圖3所示。

圖3 渦輪葉片截面劃分Fig.3 Turbine blade section division

本次計算中n=9，考慮同一截面中相鄰節點之間的影響，取整個截面的平均載荷進行計算，結果見表2。

4.2 渦輪葉片

根據文獻[7]，當T≥0.5Tm(Tm為材料熔點溫度)時，蠕變變形達到不可忽略的程度。本文忽略50%及以下熔點溫度的蠕變變形。通過MATLAB對給定的多軸蠕變曲線進行多維擬合，再將表2數據帶入進行多維線性插值，計算得到該渦輪葉片各截面各狀態下的平均蠕變應變見表3。

表2 各截面的平均溫度及等效應力Table 2 Average temperature and equivalent stress of all sections

表3 各截面不同狀態下的蠕變應變Table 3 Creep elongation on different sections under different conditions

經過5個工況后，該渦輪工作葉片各截面的平均蠕變應變見表4。

表4 各截面蠕變應變Table 4 Creep elongation on each section

本文認為，葉尖截面的蠕變應變為0，相鄰截面的平均蠕變應變為對應葉片段的蠕變應變，則相鄰截面對應段的蠕變伸長量為：

式中：hn為截面n的徑向高度(mm)。

該渦輪葉片的蠕變伸長量為：

根據以上公式，計算得到該葉片葉尖的平均蠕變伸長量為0.37 mm。

4.3 與試車結果的對比

表5 計算結果與持久試車結果的對比Table 5 Comparison between calculation results and endurance test results

表5為本文計算結果與發動機持久試車結果對比。由表中可知，持久試車渦輪葉片的蠕變伸長量分散性較大，本文計算結果與持久試車結果的平均值相對誤差為14.0%。誤差原因主要為：試驗結果為葉尖變形的最大值，其結果包括渦輪葉片工作過程中的塑性變形，本文提出的方法未考慮塑性變形，且計算結果為葉尖截面的平均蠕變伸長量。

5 結論

提出了一種基于梁理論的蠕變計算方法，并用該方法對某航空發動機燃氣渦輪工作葉片持久試車試驗在經過5個工況后的平均蠕變伸長量進行了計算。主要結論為：

(1)本文提出方法在溫度梯度較小且單純考慮蠕變伸長量時具有參考價值。

(2)持久試車渦輪葉片的蠕變伸長量分散性較大，本文提出的計算方法依賴于實測蠕變曲線的置信度，其計算結果在實測蠕變曲線的置信區間內精度更高。

(3)表3計算結果為線性插值得到，部分狀態存在外插，對計算精度有影響。

(4)本文提出的方法已應用于多種型號發動機研制中，被證明為適用性較好、工程應用性較強。