起動(dòng)過程渦輪葉片瞬態(tài)溫度場(chǎng)的準(zhǔn)穩(wěn)態(tài)處理研究

常 勝，杜 強(qiáng)，王 沛，柳 光，劉 軍，肖向濤

(1. 中國科學(xué)院 工程熱物理研究所，北京 100190；2. 中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院，北京 100049；3. 中國科學(xué)院 輕型動(dòng)力創(chuàng)新研究院，北京 100190；4. 中國科學(xué)院大學(xué) 航空宇航學(xué)院，北京 100049)

燃?xì)鉁u輪發(fā)動(dòng)機(jī)起動(dòng)、停車等過程(也稱過渡態(tài))中來流溫度和壓力的急劇變化，將會(huì)引起燃?xì)馀c渦輪葉片換熱的大幅度動(dòng)態(tài)變化，同時(shí)在金屬導(dǎo)熱影響下，葉片的溫度場(chǎng)也會(huì)處于非穩(wěn)定狀態(tài)。過渡態(tài)下異常的溫度分布會(huì)導(dǎo)致渦輪轉(zhuǎn)子承受較大的局部熱應(yīng)力及熱變形，進(jìn)而使葉片發(fā)生熱疲勞現(xiàn)象，嚴(yán)重影響渦輪部件工作壽命，因此有必要開展發(fā)動(dòng)機(jī)過渡過程渦輪葉片瞬態(tài)溫度場(chǎng)的研究，為葉片可靠性分析提供依據(jù)。

過渡過程涉及流動(dòng)穩(wěn)定和熱穩(wěn)定兩類問題，通常，建立系統(tǒng)的熱穩(wěn)定(即建立穩(wěn)定的溫度場(chǎng))所需的時(shí)間要遠(yuǎn)遠(yuǎn)大于建立流動(dòng)穩(wěn)定(即建立穩(wěn)定的速度場(chǎng))所需要的時(shí)間[1]。當(dāng)高溫燃?xì)饬鬟^渦輪葉片通道后，其熱擾動(dòng)穿透固體的深度有限，此時(shí)葉片的熱響應(yīng)將慢于氣動(dòng)響應(yīng)，應(yīng)考慮金屬熱慣性帶來的溫度滯后效應(yīng)[2]。過渡態(tài)葉片溫度的滯后會(huì)引起熱變形的滯后，隨之對(duì)葉尖間隙產(chǎn)生影響，進(jìn)一步引起渦輪效率的下降[3]，嚴(yán)重時(shí)會(huì)導(dǎo)致葉尖發(fā)生刮磨，因此須要全面地預(yù)測(cè)過渡態(tài)下金屬葉片的溫度場(chǎng)。由于發(fā)動(dòng)機(jī)整個(gè)飛行包線涉及眾多過渡過程，快速準(zhǔn)確的溫度預(yù)估方法就顯得尤為重要。

作為目前燃?xì)鉁u輪發(fā)動(dòng)機(jī)設(shè)計(jì)的關(guān)鍵技術(shù)，流熱耦合數(shù)值模擬在渦輪換熱領(lǐng)域得到了廣泛的應(yīng)用。Takahashi[4]應(yīng)用商業(yè)化軟件Fluent 5對(duì)燃?xì)鉁u輪發(fā)動(dòng)機(jī)起動(dòng)和停車過程中第一級(jí)渦輪轉(zhuǎn)子葉片進(jìn)行了瞬態(tài)流熱耦合計(jì)算，總結(jié)了來流變化條件下渦輪葉片的換熱規(guī)律：前緣和尾緣溫度變化最劇烈，葉根處溫度梯度增加較快，但對(duì)葉片內(nèi)部的溫度并沒有過多的關(guān)注。Rossette[5]在Takahashi的工作基礎(chǔ)上，研究了發(fā)動(dòng)機(jī)起動(dòng)過程中帶熱障涂層的高溫合金轉(zhuǎn)子葉片的對(duì)流換熱和導(dǎo)熱，分析了葉片中截面上熱障涂層和葉片基體的溫度梯度隨時(shí)間的變化規(guī)律，此外，將非定常項(xiàng)引入固體域的流熱耦合模型可以用來評(píng)判影響固體熱疲勞機(jī)理的熱慣性。為實(shí)現(xiàn)有限元分析(Finite Element Analysis，F(xiàn)EA)和計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)的高效耦合，Sun[6]將“SC89”插件程序應(yīng)用于發(fā)動(dòng)機(jī)四級(jí)低壓渦輪盤腔的流熱耦合計(jì)算中，數(shù)值模擬采用英國羅-羅公司開發(fā)的軟件Hydra，使用程序調(diào)用一組典型時(shí)刻CFD算例的熱邊界條件，為有限元分析軟件SC03提供輸入，在確保FEA和CFD計(jì)算之間溫度和熱流連續(xù)的前提下，實(shí)現(xiàn)二者的快速耦合。該過程相當(dāng)于采取了一種“準(zhǔn)穩(wěn)態(tài)處理”的方法，對(duì)非穩(wěn)態(tài)中典型時(shí)刻的流場(chǎng)進(jìn)行一組穩(wěn)態(tài)數(shù)值模擬，利用中間程序調(diào)用各工況下的流動(dòng)參數(shù)，再進(jìn)行固體的有限元分析，達(dá)到縮短計(jì)算時(shí)間的目的。因?yàn)榱鲌?chǎng)響應(yīng)的時(shí)間尺度遠(yuǎn)小于固體熱傳導(dǎo)的時(shí)間尺度，所以相對(duì)固體的熱慣性，流動(dòng)的瞬態(tài)效應(yīng)可以忽略。

國內(nèi)在渦輪葉片瞬態(tài)溫度場(chǎng)預(yù)測(cè)方面也開展了相關(guān)研究，錢惠華[7]針對(duì)某型渦扇發(fā)動(dòng)機(jī)高壓渦輪導(dǎo)向葉片的裂紋故障，采用自編程序計(jì)算了葉片各部分的換熱系數(shù)，結(jié)合有限元分析軟件ANSYS模擬了起動(dòng)過程50 s內(nèi)的瞬態(tài)溫度場(chǎng)，發(fā)現(xiàn)導(dǎo)向器在各瞬時(shí)的溫度最高點(diǎn)均在葉片尾緣中部，該區(qū)域溫度在9 s時(shí)已基本穩(wěn)定，認(rèn)為這是尾緣的換熱系數(shù)比較大的緣故。徐磊[8]應(yīng)用一維解析方法和數(shù)值方法分別對(duì)航空發(fā)動(dòng)機(jī)整個(gè)飛行循環(huán)內(nèi)薄壁類和軸對(duì)稱緣盤類熱端部件進(jìn)行了熱狀態(tài)計(jì)算分析，并依次研究了各零部件(葉片、輪盤及機(jī)匣)的熱響應(yīng)，為工程設(shè)計(jì)提供了一種過渡態(tài)葉尖間隙的分析方法。關(guān)鵬[9]對(duì)渦輪導(dǎo)向葉片進(jìn)行瞬態(tài)流熱耦合計(jì)算，研究了來流溫度階躍條件下的葉片外表面及內(nèi)部冷卻通道的溫度場(chǎng)，然后將溫度結(jié)果導(dǎo)入ANSYS求解器，計(jì)算了葉片的熱應(yīng)力分布與模態(tài)振型變化。屠秋野[10]以雙軸混排渦扇發(fā)動(dòng)機(jī)為研究對(duì)象，通過一種非定常熱交換模型分析了高壓轉(zhuǎn)子各部件的吸熱和放熱特點(diǎn)，得到葉片的時(shí)間常數(shù)[11]最小的結(jié)論。王伯鑫[12]提出一種葉片對(duì)流冷卻二維模型，并結(jié)合CFX對(duì)美國通用電氣公司E3發(fā)動(dòng)機(jī)高壓渦輪第一級(jí)導(dǎo)葉進(jìn)行數(shù)值模擬，將葉片表面不同位置的平均換熱系數(shù)和近壁面燃?xì)鉁囟容斎朐摾鋮s模型，實(shí)現(xiàn)了葉片壁面溫度分布變化的快速求解。

流熱耦合的計(jì)算方法雖然可以通過設(shè)置流固交界面的方式將燃?xì)鈧?cè)的對(duì)流換熱、金屬固體內(nèi)的導(dǎo)熱一并考慮，且對(duì)傳熱過程的模擬與實(shí)際情形較為接近，但就起動(dòng)過程渦輪葉片換熱計(jì)算而言，模擬全程變化的速度場(chǎng)和溫度場(chǎng)需要耗費(fèi)大量的計(jì)算資源，且計(jì)算周期長，不易實(shí)現(xiàn)葉片瞬態(tài)溫度場(chǎng)的快速求解。國內(nèi)目前對(duì)渦輪葉片熱分析中的邊界條件通常根據(jù)經(jīng)驗(yàn)公式或者約化模型預(yù)估，因此得到的結(jié)果與燃?xì)鉁u輪發(fā)動(dòng)機(jī)的真實(shí)情況存在一些偏差，且在燃?xì)鉁u輪發(fā)動(dòng)機(jī)過渡態(tài)葉片換熱問題中，對(duì)固體溫度滯后效應(yīng)的關(guān)注還比較少。本文將借鑒上文中提到的準(zhǔn)穩(wěn)態(tài)處理的思想，以期對(duì)葉片的瞬態(tài)溫度場(chǎng)進(jìn)行快速求解，掌握其熱響應(yīng)規(guī)律，為渦輪葉片的熱分析提供支撐，進(jìn)而縮短發(fā)動(dòng)機(jī)熱端部件的研制周期。

在準(zhǔn)穩(wěn)態(tài)處理過程中，為準(zhǔn)確計(jì)算渦輪葉片的溫度場(chǎng)，必須確定燃?xì)馀c葉片之間的換熱關(guān)系，這就涉及渦輪葉型外部邊界層的流動(dòng)與換熱問題。目前邊界層流動(dòng)換熱計(jì)算的方法和程序有很多種，但從葉片型面邊界層計(jì)算效果看，還是STAN5[13-14]比較適用。20世紀(jì)90年代，西北工業(yè)大學(xué)的劉松齡、朱惠人將STAN5程序引入國內(nèi)，并在Hylton[14]的工作基礎(chǔ)上，對(duì)計(jì)算程序添加了駐點(diǎn)模型、轉(zhuǎn)捩模型等內(nèi)容，提高了STAN5程序在渦輪葉片流動(dòng)換熱計(jì)算中的準(zhǔn)確性。近年來，劉松齡對(duì)STAN5程序作了進(jìn)一步的改進(jìn)，包括對(duì)駐點(diǎn)附近計(jì)算模型的改進(jìn)、加入Mayle[15]建立的轉(zhuǎn)捩模型、針對(duì)逆壓梯度和不合理速度分布造成計(jì)算失敗的現(xiàn)象制定科學(xué)的修型應(yīng)對(duì)措施等，大大改善了程序?qū)Χ喾N流動(dòng)狀態(tài)下葉片換熱系數(shù)求解的適用性，改進(jìn)后的版本命名為“NPUSTAN7Z”。

本文針對(duì)燃?xì)鉁u輪發(fā)動(dòng)機(jī)起動(dòng)過程渦輪葉片瞬態(tài)換熱問題，以某型燃?xì)鉁u輪發(fā)動(dòng)機(jī)第一級(jí)渦輪葉片為研究對(duì)象，結(jié)合商業(yè)化軟件ANSYS CFX 2019 R2和外換熱程序NPUSTAN7Z，探索了一種對(duì)過渡態(tài)進(jìn)行準(zhǔn)穩(wěn)態(tài)處理[6]的方法。準(zhǔn)穩(wěn)態(tài)處理采用非定常計(jì)算的方法求解三維Navier-Stokes(N-S)方程，提取有關(guān)流動(dòng)參數(shù)輸入外換熱程序，獲得求解固體域的第三類邊界條件，最后進(jìn)行非穩(wěn)態(tài)導(dǎo)熱計(jì)算，得到過渡態(tài)下葉片的瞬態(tài)溫度場(chǎng)。同時(shí)還對(duì)過渡態(tài)渦輪葉片進(jìn)行了瞬態(tài)流熱耦合計(jì)算，即在流體域和固體域中同時(shí)求解能量方程，流固交界面施加通量守恒條件，僅在流體域中求解輸運(yùn)方程，對(duì)起動(dòng)過程進(jìn)行了總時(shí)間3 s的非穩(wěn)態(tài)計(jì)算，最終也得到過渡態(tài)下葉片的瞬態(tài)溫度場(chǎng)。分析了葉片熱慣性帶來的溫度滯后效應(yīng)，通過對(duì)比準(zhǔn)穩(wěn)態(tài)處理和瞬態(tài)耦合計(jì)算兩種方法得到的渦輪葉片對(duì)來流熱響應(yīng)的差異，驗(yàn)證了準(zhǔn)穩(wěn)態(tài)處理的方法對(duì)過渡態(tài)渦輪葉片換熱特性預(yù)測(cè)的準(zhǔn)確性，最終實(shí)現(xiàn)對(duì)葉片瞬態(tài)溫度場(chǎng)的快速求解。

1 計(jì)算模型及方法

1.1 計(jì)算模型及邊界條件

以某型燃?xì)鉁u輪發(fā)動(dòng)機(jī)在地面全程加速過程中五級(jí)低壓渦輪第一級(jí)葉片為研究對(duì)象，約化后建立渦輪計(jì)算模型，包括3個(gè)導(dǎo)葉、5個(gè)動(dòng)葉， 如圖1所示。本文主要研究轉(zhuǎn)子葉片的瞬態(tài)溫度場(chǎng)，其幾何參數(shù)如表1所示，葉頂間隙取葉身高度的0.7%。

(a) 約化后模型

(b) 動(dòng)葉進(jìn)口速度三角形示意圖1 計(jì)算模型

表1 低壓渦輪轉(zhuǎn)子葉片幾何參數(shù)

燃?xì)夤べ|(zhì)設(shè)置為理想可壓縮氣體，其粘性系數(shù)μ和熱導(dǎo)率k隨溫度的變化滿足薩瑟蘭(Sutherland)公式[16]，設(shè)置定壓比熱cp為溫度的表達(dá)式[17]。葉片材料選擇K417G高溫合金，其導(dǎo)熱系數(shù)、比熱容等參數(shù)隨溫度變化[18]。

對(duì)燃?xì)鉁u輪發(fā)動(dòng)機(jī)地面全程加速的試驗(yàn)數(shù)據(jù)進(jìn)行?；?，模化范圍對(duì)應(yīng)試車過程的中間14 s，以模化后來流條件變化率最大的3 s(6.2～9.2)為研究區(qū)間，截取第一級(jí)低壓渦輪進(jìn)口總溫Tt、進(jìn)口總壓pt、出口靜壓p和轉(zhuǎn)速值N，以起始時(shí)刻第6.2 s為參考對(duì)各變量進(jìn)行無量綱處理，得到隨時(shí)間變化的邊界條件相對(duì)值C/C6.2(C代表Tt、pt、p、N),如圖2所示。

圖2 起動(dòng)過程3 s內(nèi)邊界條件的變化曲線

燃?xì)饨?jīng)過渦輪過渡段后，氣流角沿徑向非均勻分布，根據(jù)測(cè)試數(shù)據(jù)，進(jìn)口各速度分量沿徑向的變化如圖3所示,進(jìn)口湍流度取5%。

圖3 主流進(jìn)口沿徑向的速度分量

1.2 瞬態(tài)流熱耦合邊界設(shè)置

瞬態(tài)流熱耦合計(jì)算采用商業(yè)化軟件ANSYS CFX 2019 R2進(jìn)行數(shù)值模擬，計(jì)算域包含導(dǎo)葉流體域、動(dòng)葉流體域、導(dǎo)葉固體域和動(dòng)葉固體域四個(gè)計(jì)算域，而非定常計(jì)算域僅包含導(dǎo)葉流體域和動(dòng)葉流體域兩部分，相對(duì)流熱耦合計(jì)算，非定常計(jì)算不含固體域，圖4(a)、4(b)分別為瞬態(tài)流熱耦合計(jì)算域和非定常計(jì)算域。

流熱耦合計(jì)算中進(jìn)出口邊界條件按照?qǐng)D2與圖3給出的數(shù)據(jù)設(shè)置，以起始時(shí)刻的邊界條件進(jìn)行非定常流熱耦合計(jì)算，將計(jì)算結(jié)果作為瞬態(tài)流熱耦合和準(zhǔn)穩(wěn)態(tài)計(jì)算的初場(chǎng)。轉(zhuǎn)靜交界面設(shè)置為Transient Rotor Stator混合模型，葉片表面(含葉尖)設(shè)置為通量守恒的流固交界面，除交界面外所有壁面均滿足無滑移條件?？紤]葉根榫槽內(nèi)充滿冷卻氣體，3 s加速過程葉根設(shè)置為800 K恒定的溫度，如圖4(c)所示，輪轂、機(jī)匣(上、下面)為絕熱條件，兩側(cè)(左、右面)設(shè)置周期交界面。

(a) 耦合計(jì)算域

(b) 非定常計(jì)算域

(c) 各邊界位置

1.3 準(zhǔn)穩(wěn)態(tài)處理計(jì)算過程

準(zhǔn)穩(wěn)態(tài)處理的計(jì)算方法結(jié)合了商業(yè)化軟件ANSYS CFX 2019 R2和外換熱程序NPUSTAN7Z，分以下幾個(gè)步驟來開展：CFX求解流體域、數(shù)據(jù)預(yù)處理、表面換熱系數(shù)計(jì)算、三維節(jié)點(diǎn)插值等，其流程如圖5所示。

圖5 準(zhǔn)穩(wěn)態(tài)處理流程圖

1.3.1 非定常計(jì)算邊界設(shè)置

對(duì)過渡態(tài)進(jìn)行準(zhǔn)穩(wěn)態(tài)處理，選取起動(dòng)過程研究區(qū)間3 s內(nèi)的7個(gè)典型工況，各工況點(diǎn)間隔約0.5 s，分別對(duì)應(yīng)的時(shí)刻是A:第6.2 s、B:第6.699 5 s、C: 第7.200 5 s、D:第7.7 s、E:第8.199 5 s、F:第8.700 5 s、G:第9.2 s，以上述時(shí)刻的邊界條件進(jìn)行非定常數(shù)值模擬。非定常計(jì)算域如上文描述，設(shè)置所有壁面為無滑移絕熱條件，轉(zhuǎn)靜交界面設(shè)置為Transient Rotor Stator混合模型，左右兩側(cè)設(shè)置周期交界面。

1.3.2 計(jì)算表面換熱系數(shù)

表面換熱系數(shù)的計(jì)算是準(zhǔn)穩(wěn)態(tài)處理中的關(guān)鍵步驟。采用非定常的方法求解時(shí)均流場(chǎng)后，提取相關(guān)流動(dòng)參數(shù)輸入到外換熱程序中，程序采用流函數(shù)坐標(biāo)變換法，求解邊界層微分方程，計(jì)算葉片表面換熱系數(shù)，作為邊界條件，進(jìn)一步再進(jìn)行葉片瞬態(tài)溫度場(chǎng)的求解。

本文提取轉(zhuǎn)子葉片表面10%、20%、…、90%共9個(gè)展向高度上的型線數(shù)據(jù)輸入NPUSTAN7Z程序求解邊界層方程，如圖6所示為不同時(shí)刻外換熱程序計(jì)算的部分展向高度上沿型面分布的表面換熱系數(shù)，橫軸表示無量綱的型面弧長，0表示前緣，區(qū)間(-1,0)表示壓力面，區(qū)間(0,1)表示吸力面。由于前緣受到主流燃?xì)獾臎_擊，邊界層非常薄，因此具有較高的換熱系數(shù)。沿壓力面和吸力面型線至尾緣邊界層逐漸增厚，相應(yīng)地?fù)Q熱系數(shù)較低。隨著發(fā)動(dòng)機(jī)起動(dòng)過程中來流條件的改變，表面換熱系數(shù)也逐漸增加。

1.3.3 三維節(jié)點(diǎn)插值

進(jìn)行渦輪葉片導(dǎo)熱計(jì)算時(shí)須要在網(wǎng)格節(jié)點(diǎn)上加載邊界條件，而葉片表面的網(wǎng)格節(jié)點(diǎn)與外換熱程序計(jì)算的9條已知數(shù)據(jù)型線的坐標(biāo)并不一致，不能滿足三維溫度場(chǎng)計(jì)算的需要[19]，進(jìn)一步還應(yīng)將9條型線上已知的換熱系數(shù)插值到固體域的網(wǎng)格節(jié)點(diǎn)上。

(a) 70%展向高度

(b) 50%展向高度

(c) 30%展向高度圖6 隨時(shí)間變化的表面換熱系數(shù)

本文采用MATLAB編寫了三維彎扭葉片表面節(jié)點(diǎn)數(shù)據(jù)的插值程序，采用K-近鄰算法[20]，按照歐式距離，尋找距葉片表面網(wǎng)格節(jié)點(diǎn)最近的已知型線上的3～4個(gè)坐標(biāo)點(diǎn)，以坐標(biāo)點(diǎn)間的距離確定權(quán)重，將相應(yīng)坐標(biāo)點(diǎn)上的換熱系數(shù)插值到葉片表面節(jié)點(diǎn)上。尋找近鄰點(diǎn)時(shí)，兩型線中間的節(jié)點(diǎn)分別尋找上、下型線上距離最近的2點(diǎn)，葉根及葉頂區(qū)域?qū)ふ艺瓜?0%、展向90%型線上的3點(diǎn)，這樣可以避免已知型線附近的網(wǎng)格節(jié)點(diǎn)只尋找了該條型線上的數(shù)據(jù)，造成局部換熱系數(shù)偏差過大。如圖7為初始時(shí)刻9條型線上的換熱系數(shù)插值出的三維葉片表面所有節(jié)點(diǎn)的換熱系數(shù)分布圖。

圖7 初始時(shí)刻插值后的葉片表面換熱系數(shù)分布

1.3.4 求解葉片瞬態(tài)溫度場(chǎng)

依次對(duì)準(zhǔn)穩(wěn)態(tài)處理的工況A至G進(jìn)行求解流體域、計(jì)算表面換熱系數(shù)、進(jìn)行節(jié)點(diǎn)插值，得到7個(gè)時(shí)刻下非穩(wěn)態(tài)導(dǎo)熱計(jì)算離散點(diǎn)的第三類邊界條件，接著對(duì)節(jié)點(diǎn)上相鄰兩個(gè)時(shí)刻的表面換熱系數(shù)和絕熱壁溫取平均，獲得中間6個(gè)時(shí)刻的邊界條件，計(jì)算步驟如表2所示。3 s過渡過程的瞬態(tài)溫度場(chǎng)計(jì)算分12個(gè)區(qū)間進(jìn)行，時(shí)間步長設(shè)置為7.5×10-4s，初次計(jì)算時(shí)與耦合計(jì)算設(shè)置相同的初場(chǎng)，之后每次計(jì)算均以上次的結(jié)果為初場(chǎng)，轉(zhuǎn)子葉根處給定800 K均溫，最終得到過渡態(tài)葉片內(nèi)部的瞬態(tài)溫度場(chǎng)。

表2 葉片溫度場(chǎng)計(jì)算步驟

1.4 網(wǎng)格劃分

所有計(jì)算域均采用ICEM CFD劃分結(jié)構(gòu)化網(wǎng)格，為減少插值帶來的誤差，流體域和固體域的交界面上設(shè)置了相互匹配(相同數(shù)目及分布規(guī)律)的網(wǎng)格，在固體域周邊布置O型網(wǎng)格，保證了沿葉片表面網(wǎng)格的正交性。

劃分了350 萬、709 萬、1 001 萬三種不同數(shù)目的網(wǎng)格進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證，各網(wǎng)格劃分方案如表3所示。按照初始時(shí)刻第6.2 s的邊界條件做非定常流熱耦合計(jì)算，通過對(duì)比三種方案得到的瞬態(tài)初場(chǎng)的溫度、壓力值，綜合考慮計(jì)算精度和所需成本，網(wǎng)格劃分選擇方案2，流熱耦合與非定常計(jì)算采用同樣的流體域網(wǎng)格。

表3 網(wǎng)格劃分方案

渦輪動(dòng)葉的網(wǎng)格細(xì)節(jié)如圖8所示,壁面第一層網(wǎng)格的尺度為0.001 mm，使得壁面無量綱單位y+小于1，滿足湍流模型的要求。為使流固交界面上網(wǎng)格具有較好的匹配性，固體域也采用O型網(wǎng)格，壁面以內(nèi)第一層網(wǎng)格尺度設(shè)置為0.008 mm，沿葉片型線的網(wǎng)格數(shù)目和分布規(guī)律相同，即交界面兩側(cè)網(wǎng)格密度相同。葉頂間隙內(nèi)沿徑向布置25層網(wǎng)格，兩端網(wǎng)格尺度設(shè)置為0.001 mm，輪轂壁面采取同樣的處理方式。

圖8 動(dòng)葉流體域與固體域的網(wǎng)格細(xì)節(jié)

2 數(shù)值方法

2.1 控制方程

2.1.1 CFX

數(shù)值模擬采用商業(yè)化軟件ANSYS CFX 2019 R2，其核心求解器CFX Solver中，固定的笛卡爾坐標(biāo)系(i，j，k)下，微分形式的連續(xù)方程、動(dòng)量方程和能量方程：

(1)

(2)

(3)

式中：SM為動(dòng)量源項(xiàng)；SE為能量源項(xiàng)；h*為比總焓。且有

(4)

因?yàn)樵诠腆w區(qū)域內(nèi)沒有流動(dòng)存在，可以將對(duì)流項(xiàng)和擴(kuò)散項(xiàng)取消，并且在固體域內(nèi)只有熱傳導(dǎo)傳熱方式存在，所以在固體區(qū)域中能量守恒方程簡(jiǎn)化為：

(5)

CFX軟件對(duì)流體域和固體域均采用有限體積法求解，在流體域中求解三維可壓縮雷諾時(shí)均N-S方程，在固體域中僅求解能量方程。瞬態(tài)流熱耦合計(jì)算中，固體域和流體域滿足能量方程的同時(shí)，在流固交界面上施加通量守恒的條件。

耦合計(jì)算和非定常計(jì)算中，湍流模型均選擇基于雷諾平均方法的二方程模型——剪切應(yīng)力輸運(yùn)(Shear Stress Transport, SST)格式的k-ω模型，求解過程中，控制方程對(duì)流項(xiàng)均采用高精度差分格式(High Resolution Scheme)，時(shí)間項(xiàng)采用隱式歐拉二階向后差分格式(Second Order Backward Euler Scheme)。

2.1.2 NPUSTAN7Z

外換熱程序NPUSTAN7Z是一個(gè)邊界層有限差分的計(jì)算程序，通過求解壁面邊界層微分方程，得到渦輪葉片表面的換熱系數(shù)。

邊界層內(nèi)連續(xù)方程：

(6)

x方向的動(dòng)量方程：

(7)

式中：τ為邊界層內(nèi)流體的粘性應(yīng)力。

y方向的動(dòng)量方程：

(8)

總焓形式的能量方程：

(9)

式中：q為單位面積的熱流量；X為單位質(zhì)量流體的徹體力，對(duì)于氣體?？梢月匀?；I為流體的總焓。

外換熱程序中采用的湍流模型有混合長度模型、單方程湍流模型，同時(shí)考慮了湍流度修正。其中較為常用的是普朗特(Prandtl)雙方程混合長度模型。

2.2 計(jì)算結(jié)果驗(yàn)證

非定常計(jì)算的時(shí)間步長設(shè)置為1×10-5s，此時(shí)各典型時(shí)刻動(dòng)葉流體域單個(gè)扇區(qū)(6.372°)通過周期內(nèi)包含49～73個(gè)時(shí)間步，能夠保證計(jì)算結(jié)果的準(zhǔn)確性。如圖9為7個(gè)工況計(jì)算的動(dòng)葉出口的溫度監(jiān)控曲線，受到上游尾跡影響，當(dāng)動(dòng)葉出口的溫度表現(xiàn)出明顯的非定常特性時(shí)，證明計(jì)算已收斂。每次計(jì)算收斂后，對(duì)200步的流場(chǎng)結(jié)果取時(shí)均值，用于外換熱程序的輸入。

圖9 非定常計(jì)算動(dòng)葉出口監(jiān)控點(diǎn)的溫度

圖10 耦合計(jì)算葉片前緣監(jiān)控點(diǎn)的溫度

瞬態(tài)流熱耦合計(jì)算總時(shí)間較長，需要消耗大量的計(jì)算資源，選擇了三種時(shí)間步長1.5×10-2s、1.5×10-3s、1.5×10-4s進(jìn)行無關(guān)性驗(yàn)證，對(duì)起動(dòng)過程進(jìn)行總時(shí)間0.5 s的瞬態(tài)流熱耦合計(jì)算，考察轉(zhuǎn)子葉片中截面附近前緣一點(diǎn)的溫度值，如圖10所示，可以看到隨著時(shí)間步長的增加，監(jiān)控點(diǎn)的溫度隨時(shí)間的推移逐漸起伏變化，這與圖2中來流條件的平順變化規(guī)律不符，當(dāng)時(shí)間步長為1.5×10-3s 和1.5×10-4s時(shí)，溫度變化趨勢(shì)更接近真實(shí)情況，但綜合考慮耦合計(jì)算的經(jīng)濟(jì)性，最終確定1.5×10-3s的時(shí)間步長。

3 結(jié)果及分析

3.1 過渡態(tài)葉片表面瞬態(tài)溫度分布

通過準(zhǔn)穩(wěn)態(tài)計(jì)算與瞬態(tài)耦合計(jì)算兩種方法，得到轉(zhuǎn)子葉片的瞬態(tài)溫度場(chǎng)，如圖11為壓力面和吸力

(a) 壓力面瞬態(tài)溫度場(chǎng)

(b) 吸力面瞬態(tài)溫度場(chǎng)圖11 耦合計(jì)算和準(zhǔn)穩(wěn)態(tài)計(jì)算的瞬態(tài)溫度場(chǎng)

面溫度隨時(shí)間的變化云圖。在發(fā)動(dòng)機(jī)加速3 s時(shí)間內(nèi)，隨著來流溫度、壓力和轉(zhuǎn)速的增加，轉(zhuǎn)子葉片表面溫度也隨之上升，并且呈現(xiàn)出前緣和尾緣溫度升高快于葉中、由葉尖至葉根溫度逐漸降低、葉根附近存在較大的溫度梯度等特點(diǎn)。

(a) 70%徑向高度

(b) 50%徑向高度

(c) 30%徑向高度圖12 不同徑向高度型面上的溫度變化規(guī)律

對(duì)型面弧長坐標(biāo)進(jìn)行無量綱化，得到葉片表面不同徑向高度型線上溫度的變化曲線，如圖12所示。隨著時(shí)間的推移，葉片表面不同徑向高度上的溫度逐漸升高，前緣和尾緣相對(duì)中弦位置溫度升高更快。由于耦合計(jì)算的結(jié)果受到流場(chǎng)的瞬態(tài)效應(yīng)(尾跡、二次流等)影響，沿弧長的溫度分布相比準(zhǔn)穩(wěn)態(tài)的計(jì)算結(jié)果波動(dòng)較大。兩種計(jì)算方法得到的表面溫度分布存在一定差異，末尾時(shí)刻溫度差較為顯著，在70%、50%、30%徑向高度上，二者計(jì)算的型線平均溫差為6.7 K、10.1 K、8.5 K，前緣附近的差異尤為顯著，末尾時(shí)刻準(zhǔn)穩(wěn)態(tài)的結(jié)果高于耦合計(jì)算14.4 K、18.4 K、16.3 K。這源于兩種計(jì)算方法耦合方式的不同，準(zhǔn)穩(wěn)態(tài)通過求解邊界層方程得到表面換熱系數(shù)進(jìn)而計(jì)算溫度場(chǎng)，流熱耦合計(jì)算采用設(shè)置交界面上通量守恒的方式得到壁面溫度。但總體上二者溫度變化趨勢(shì)接近，說明準(zhǔn)穩(wěn)態(tài)處理方法具有一定的準(zhǔn)確性。

3.2 葉片內(nèi)部瞬態(tài)溫度場(chǎng)

3.2.1 不同區(qū)域的瞬態(tài)溫度

為詳細(xì)了解葉片內(nèi)部的溫度變化趨勢(shì)，查看葉片30%、50%、70%無量綱徑向高度截面上的溫度值，圖13為7個(gè)典型時(shí)刻兩種計(jì)算方法求解的葉片內(nèi)部的溫度場(chǎng)云圖。由圖可知，葉片內(nèi)溫度隨著時(shí)間的推移而增加，非穩(wěn)態(tài)導(dǎo)熱對(duì)葉片內(nèi)溫度場(chǎng)產(chǎn)生一定影響。隨著來流條件的改變，熱擾動(dòng)由壁面逐漸向內(nèi)傳遞，葉片內(nèi)部的熱響應(yīng)慢于葉片表面。此外，轉(zhuǎn)子葉片榫槽內(nèi)被冷卻空氣環(huán)繞，相對(duì)葉身溫度較低，因此葉片內(nèi)部存在由葉身向葉根的導(dǎo)熱，由葉尖至葉根溫度逐漸降低，進(jìn)而造成葉片內(nèi)部沿徑向的熱響應(yīng)不一致，葉片中下部對(duì)燃?xì)獾臒犴憫?yīng)相對(duì)較慢。

圖13 葉片不同徑向截面溫度云圖

葉片內(nèi)部的等值線方向接近于橫向，表明葉片內(nèi)部存在縱向熱傳導(dǎo)，熱量由前緣、尾緣向中弦位置傳遞，沿中弧線方向存在由中段指向前緣、尾緣的溫度梯度，最終在葉片中段形成低溫區(qū)。準(zhǔn)穩(wěn)態(tài)和耦合計(jì)算得到的葉片內(nèi)部溫度變化趨勢(shì)仍然是接近的，其中準(zhǔn)穩(wěn)態(tài)處理預(yù)測(cè)的溫度場(chǎng)略高。

圖14為葉片內(nèi)部9條特征弧線(如圖15中從壓力面穿過葉片至吸力面的圓弧段)上的溫度隨時(shí)間的變化曲線，反映了葉片從壓力側(cè)到吸力側(cè)的橫向溫度變化規(guī)律，從左至右，從上至下依次為10%、50%、90%無量綱軸向位置和70%、50%、30%無量綱徑向高度上的溫度值，其中橫軸表示無量綱的周向弧長，-1表示壓力側(cè)，(-1，1)表示葉片內(nèi)部，1表示吸力側(cè)。10%、50%、90%軸向位置可以用來表征葉片前緣、葉片中弦位置、葉片尾緣的換熱特性。

在第6.699 5 s時(shí)刻，各個(gè)位置上的溫度均已偏離初始值，因此經(jīng)過0.5 s后熱擾動(dòng)已穿透整個(gè)葉片。各特征位置上的熱響應(yīng)逐漸加快，如圖中70%徑向高度上前緣附近在開始的0.5 s內(nèi)，耦合計(jì)算(圖14(a)中實(shí)線橢圓圈出)的溫度平均增加了1.5 K，末尾的0.5 s內(nèi)溫度平均增加了10.7 K。又如50%徑向高度上中弦位置在開始的0.5 s內(nèi)，耦合計(jì)算(圖14(b)中虛線橢圓圈出)的溫度平均增加了0.1 K，末尾的0.5 s內(nèi)溫度平均增加了7.9 K。這是由于起動(dòng)過程中來流雷諾數(shù)的增加帶來葉片表面換熱系數(shù)的增加，同時(shí)受到來流溫度升高的影響，導(dǎo)致越來越多的熱量傳遞到葉片內(nèi)部，由此熱響應(yīng)的也越來越快。

(a) 70%徑向高度

(b) 50%徑向高度

(c) 30%徑向高度圖14 葉片內(nèi)部不同特征弧線上的溫度值

圖15 葉片無量綱相對(duì)位置示意圖

觀察葉片各處的溫度變化規(guī)律可以發(fā)現(xiàn)，葉片的溫度從初場(chǎng)開始升高，之后約0.5 s時(shí)間內(nèi)，各軸向位置幾乎不存在橫向溫差，隨著時(shí)間的推移，葉片壓力側(cè)、吸力側(cè)溫度差開始增加。如70%徑向高度上前緣附近在第6.2 s和第6.699 5 s時(shí)刻的橫向溫差為0 K，在第8.700 5 s和第9.2 s時(shí)刻，耦合計(jì)算的橫向溫差變?yōu)?.0 K和3.0 K。又如50%徑向高度上中弦位置在第6.2 s和第6.699 5 s時(shí)刻的橫向溫差也是0 K，在第8.700 5 s和第9.2 s時(shí)刻，耦合計(jì)算的橫向溫差變?yōu)?.7 K和1.9 K，其余位置也具有同樣的現(xiàn)象。起動(dòng)過程中，由于渦輪葉片吸力面和壓力面流動(dòng)狀態(tài)的不一致，使得兩側(cè)的換熱特性不同，導(dǎo)致葉片內(nèi)部存在一定的橫向溫度梯度，且溫度梯度逐漸增加。

在50%軸向位置上，接近末尾時(shí)刻可以觀察到溫度呈中間低兩端高的現(xiàn)象，說明中弧線附近區(qū)域?qū)θ細(xì)鉁囟鹊臒犴憫?yīng)較慢，由于前緣、尾緣的橫向尺寸較小，熱擾動(dòng)極短的時(shí)間內(nèi)就能將其穿透，因此沒有形成明顯的低溫區(qū)。對(duì)比不同徑向高度前緣、中弦及尾緣溫度隨時(shí)間的變化曲線，發(fā)現(xiàn)三者中中弦位置溫度隨時(shí)間變化最慢，相對(duì)來說，局部熱容量較小的前緣、尾緣對(duì)燃?xì)鉁囟茸兓臒犴憫?yīng)比葉片中弦位置更快。準(zhǔn)穩(wěn)態(tài)計(jì)算的溫度略高，但在溫度變化趨勢(shì)上與流熱耦合的結(jié)果比較一致，對(duì)得比較好。

3.2.2 葉片內(nèi)部平均溫度及其滯后效應(yīng)

為定量分析葉片內(nèi)部的溫度變化趨勢(shì)，對(duì)30%、50%、70%徑向截面的溫度取平均值，得到不同截面的溫度隨時(shí)間的變化曲線，如圖16所示。圖中給出了過渡態(tài)加速3 s過程內(nèi)來流截面(位于葉片上游軸向弦長30%處，如圖7所示)平均的溫度變化，耦合計(jì)算的來流溫度升高了93.0 K(非定常計(jì)算初、末兩個(gè)算例來流溫升93.4 K)，30%、50%、70%徑向高度上溫度分別升高了27.3 K、27.2 K、33.9 K，準(zhǔn)穩(wěn)態(tài)計(jì)算的溫升分別是35.5 K、37.1 K、38.7 K，二者計(jì)算的溫升存在一定差異，但計(jì)算結(jié)果共同反映了葉片中上部相對(duì)中下部熱響應(yīng)更快的現(xiàn)象，這是由葉片中下部熱傳導(dǎo)帶走更多的熱量導(dǎo)致的。由于金屬葉片本身具有熱慣性，葉片內(nèi)部不同截面上的平均溫度隨時(shí)間的變化率低于來流溫度隨時(shí)間的變化率，說明葉片的溫升慢于來流的溫升，因此葉片內(nèi)部存在溫度滯后效應(yīng)。

圖16 不同徑向截面平均溫度的變化曲線

準(zhǔn)穩(wěn)態(tài)計(jì)算與瞬態(tài)流熱耦合計(jì)算的結(jié)果總體上吻合較好，需要說明的是，來流溫度初始值的不同是因?yàn)槎叩挠?jì)算域不同，準(zhǔn)穩(wěn)態(tài)初溫取來流截面時(shí)均溫度的平均值，流熱耦合初溫取來流截面瞬態(tài)溫度的平均值，二者計(jì)算域及邊界條件不一致，所以初始時(shí)刻存在偏差。

以初始時(shí)刻的溫度為參考值，定義燃?xì)鉁厣俜直群徒饘贉厣俜直龋?/p>

(10)

(11)

式中：tg0為燃?xì)獬跏紩r(shí)刻溫度；tg為燃?xì)鉁囟龋沪g為燃?xì)鉁囟认鄬?duì)初始時(shí)刻溫度之差；ts0為葉片初始時(shí)刻溫度，葉片內(nèi)部各處初始溫度不盡相同；ts為葉片溫度；Δts為葉片溫度相對(duì)初始時(shí)刻溫度之差。

圖17為不同徑向高度上面平均的溫升隨時(shí)間的變化曲線，圖中給出了過渡態(tài)加速3 s過程內(nèi)來流截面平均的溫升曲線，可以看出，來流溫度升高了11.1%(非定常計(jì)算初、末兩個(gè)算例來流溫升11.3%)，耦合計(jì)算和準(zhǔn)穩(wěn)態(tài)計(jì)算的30%、50%、70%徑向高度面平均的溫升分別是2.7%、2.7%、3.3%和3.5%、3.6%、3.8%，二者計(jì)算的溫升百分比偏差最大處位于葉片中部和下部，分別為0.8%、0.9%，偏差較小，因此認(rèn)為準(zhǔn)穩(wěn)態(tài)處理的方法能夠準(zhǔn)確預(yù)測(cè)過渡態(tài)渦輪葉片的瞬態(tài)溫度場(chǎng)。

圖17 不同徑向截面的平均溫升

葉片中下部、中部、中上部溫度相對(duì)來流溫度分別滯后了8.4%、8.4%、7.8%(耦合)，所以中上部的熱響應(yīng)相對(duì)中下部較快，葉片中下部區(qū)域滯后效應(yīng)較強(qiáng)，考慮是因?yàn)槿~根溫度相對(duì)較低，葉片中下部溫度場(chǎng)受熱傳導(dǎo)的影響更大(觀察準(zhǔn)穩(wěn)態(tài)計(jì)算的結(jié)果，可以得到相同的結(jié)論)。

圖18為葉片內(nèi)部不同位置上的線平均溫度滯后百分比的柱狀圖，受葉片熱慣性的影響，葉片內(nèi)部不同區(qū)域的溫度均會(huì)存在滯后，且溫度滯后程度不一樣。聯(lián)系上文分析，葉片較薄位置(前緣、尾緣)相對(duì)葉片較厚位置(中弦區(qū)域)熱響應(yīng)更快，葉片中下部的溫度滯后效應(yīng)較強(qiáng)。對(duì)比圖18(a)中前緣、中弦和尾緣附近的溫度滯后，可知中弦滯后>前緣滯后>尾緣滯后。對(duì)比圖18(a)中耦合計(jì)算不同徑向高度上的溫度滯后，可知葉片中下部滯后>葉片上部滯后。

葉片內(nèi)局部溫度滯后最大的位置在50%軸向弦長處的中下部區(qū)域，最大滯后達(dá)到9.6%，因?yàn)樵搮^(qū)域靠近葉根，同時(shí)也是葉片中較厚的位置，所以該區(qū)域相對(duì)來流的熱響應(yīng)最慢。

準(zhǔn)穩(wěn)態(tài)計(jì)算結(jié)果的溫升普遍偏高，對(duì)比圖18(a)與圖18(b)不同位置上準(zhǔn)穩(wěn)態(tài)與耦合計(jì)算的滯后百分?jǐn)?shù)偏差，前緣附近70%、50%、30%徑向高度上的偏差分別為0.7%、1.2%、0.8%，中弦位置分別為0、0.9%、0.8%，尾緣附近分別為0.1%、0.2%、0.4%，準(zhǔn)穩(wěn)態(tài)處理的計(jì)算結(jié)果雖然在各局部位置上均低估了金屬葉片的溫度滯后，但相對(duì)耦合計(jì)算的溫度場(chǎng)偏差不大。

(a) 耦合

(b) 準(zhǔn)穩(wěn)態(tài)

4 結(jié)語

本文分別采用瞬態(tài)流熱耦合計(jì)算和準(zhǔn)穩(wěn)態(tài)處理的方法求解過渡態(tài)渦輪葉片的瞬態(tài)溫度場(chǎng)，通過比較二者的計(jì)算結(jié)果，發(fā)現(xiàn)準(zhǔn)穩(wěn)態(tài)處理與瞬態(tài)耦合計(jì)算的結(jié)果吻合較好，僅略微低估了金屬葉片的溫度滯后，因此認(rèn)為準(zhǔn)穩(wěn)態(tài)處理可以實(shí)現(xiàn)過渡態(tài)渦輪葉片換熱特性的快速求解。同時(shí)開發(fā)了三維彎扭葉片的節(jié)點(diǎn)數(shù)據(jù)插值程序，解決了工程實(shí)際中渦輪葉片瞬態(tài)溫度場(chǎng)計(jì)算時(shí)的邊界條件銜接問題。分析葉片的瞬態(tài)溫度場(chǎng)得到如下結(jié)論：

(1) 非穩(wěn)態(tài)導(dǎo)熱會(huì)對(duì)葉片內(nèi)溫度場(chǎng)產(chǎn)生一定影響，葉片內(nèi)部的熱響應(yīng)慢于葉片表面。不同軸向位置溫度滯后程度的不一致，形成了葉片內(nèi)部的縱向溫度梯度，葉片中段存在低溫區(qū)。吸力面和壓力面流動(dòng)狀態(tài)的不一致導(dǎo)致葉片內(nèi)部存在一定的橫向溫度梯度。

(2) 熱擾動(dòng)在較短的時(shí)間內(nèi)已穿透整個(gè)葉片，起動(dòng)過程中，隨著時(shí)間的推移，來流溫度與表面換熱系數(shù)的共同作用使得葉片內(nèi)的熱響應(yīng)逐漸加快。

(3) 葉片本身的熱慣性導(dǎo)致其溫升慢于來流的溫升，在發(fā)動(dòng)機(jī)起動(dòng)過程的3 s內(nèi)，葉片中截面平均溫度相對(duì)來流滯后了8.4%。

(4) 葉片不同軸向位置的厚度不一致，帶來了不同區(qū)域局部熱容量的差異，反映到轉(zhuǎn)子葉片溫度場(chǎng)上的就是中弦滯后>前緣滯后>尾緣滯后。

(5) 由于葉根是低溫區(qū)，離葉根越近熱傳導(dǎo)帶走的熱量越多，所以葉片中下部滯后>葉片上部滯后，葉根附近存在較大的溫度梯度。