基于θ映射法的燃?xì)鉁u輪葉片高溫蠕變變形分析

全昌彪，廖明夫，米 棟，李 堅(jiān)，劉 揚(yáng)

(1.西北工業(yè)大學(xué)動(dòng)力與能源學(xué)院，西安710072；2.中國(guó)航發(fā)湖南動(dòng)力機(jī)械研究所，湖南株洲412002)

1 引言

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

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

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

2 蠕變模型

2.1 單軸蠕變方程

對(duì)于各向同性材料，蠕變應(yīng)變?chǔ)與通常可表示為構(gòu)件所受應(yīng)力σ、溫度T和時(shí)間t的函數(shù)，即：

如果給定試驗(yàn)條件(σ，T)時(shí)，蠕變應(yīng)變只是時(shí)間t的函數(shù)：

Evans和Wilshire認(rèn)為蠕變過(guò)程可描述如下[5]：

式中：θi(i=1～4)是與材料、溫度及應(yīng)力有關(guān)的常數(shù)。對(duì)于特定材料，θi可表示為應(yīng)力和溫度的函數(shù)：

式中：ai、bi、ci、di為材料特性相關(guān)的參數(shù)，可根據(jù)蠕變?cè)囼?yàn)曲線進(jìn)行優(yōu)化獲取。

2.2 多軸蠕變方程

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

式中：λ為塑性乘子，Sij為應(yīng)力偏張量，δij為Kronecker Delta函數(shù)，σkk為3個(gè)主應(yīng)力分量之和。

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

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

式中：為等效應(yīng)力，其表達(dá)式如下：

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

3 蠕變模型參數(shù)的確定

DZ408材料在同一溫度下進(jìn)行了多種應(yīng)力水平下的蠕變?cè)囼?yàn)，表1匯總了其蠕變?cè)囼?yàn)條件。本文根據(jù)其蠕變?cè)囼?yàn)數(shù)據(jù)，采用最小二乘法按照式(3)和式(4)擬合出其中的蠕變參數(shù)，結(jié)果見(jiàn)表2和表3。

表1 DZ408材料蠕變曲線試驗(yàn)條件Table 1 Creep test conditions of DZ408

表2 θ參數(shù)擬合結(jié)果Table 2 Fitting results of parameterθ

表3 a、b、c、d參數(shù)擬合結(jié)果Table 3 Fitting results of parametera,b,c,d

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

4 算例分析

對(duì)某型燃?xì)鉁u輪葉片60 h整機(jī)持久試車過(guò)程的蠕變行為進(jìn)行了近似模擬。圖2為簡(jiǎn)化處理后的葉片有限元網(wǎng)格模型。網(wǎng)格采用四節(jié)點(diǎn)四面體單元(C3D4）劃分，共計(jì)148 582個(gè)單元，40 188個(gè)節(jié)點(diǎn)。

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

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

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

表5 各分析步位移計(jì)算結(jié)果Table 5 Displacement calculation results of each procedure

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

葉片在整個(gè)蠕變過(guò)程中，葉身應(yīng)力重新分布。圖5給出了葉身中部位置N1341和N696節(jié)點(diǎn)當(dāng)量應(yīng)力隨蠕變時(shí)間推移而發(fā)生的應(yīng)力松弛現(xiàn)象。由圖可看出，應(yīng)力松弛現(xiàn)象使得葉身高應(yīng)力區(qū)的應(yīng)力有所降低，低應(yīng)力區(qū)的應(yīng)力有所增加，整個(gè)葉片的應(yīng)力分布趨向均勻化，這有益于葉片的使用壽命。

5 計(jì)算變形與實(shí)測(cè)伸長(zhǎng)量對(duì)比

卸載后葉尖的最大徑向位移為0.086 94 mm。發(fā)動(dòng)機(jī)在3次60 h持久試車后，渦輪葉片平均伸長(zhǎng)量的實(shí)測(cè)值分別為0.058、0.083、0.080 mm，平均值為0.074 mm。計(jì)算結(jié)果與實(shí)測(cè)結(jié)果平均值相對(duì)誤差約18.02%，其中有兩次不到10.00%，對(duì)比結(jié)果見(jiàn)表6。

表6 葉片伸長(zhǎng)量計(jì)算值和實(shí)測(cè)值對(duì)比Table 6 Comparison between blade elongation calculation results and test results

6 結(jié)論

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

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

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

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