施能博， 張溫宇， 孫志新， 林 立

(福州大學石油化工學院， 福建 福州 350116)

0 引言

燃氣輪機的效率提升主要通過提高渦輪進口溫度實現(xiàn)， 目前先進燃氣輪機的渦輪進口溫度已遠超過渦輪材料的耐受溫度， 因此必須對高溫部件進行有效的冷卻， 以保證各部件的最高溫度與溫度梯度在允許范圍內(nèi).

傳統(tǒng)關(guān)于內(nèi)部冷卻空氣系統(tǒng)的研究多專注于設(shè)計工況或非設(shè)計工況的穩(wěn)態(tài)分析， 隨著研發(fā)水平的不斷提高， 非定常效應(yīng)逐漸成為研究的關(guān)注點. Zimmer等[1]考察了將非定常沖擊射流引入燃氣輪機透平葉片冷卻的效果， 結(jié)果表明斯特勞哈數(shù)小于一定值的情況下， 非定常射流沖擊效果弱于穩(wěn)態(tài)射流. Scherhag等[2]通過實驗和數(shù)值研究， 考察了脈沖射流沖擊對透平冷卻應(yīng)用的可行性及其對各類參數(shù)的敏感程度，脈沖來源為動葉造成的周期性非定常壓力分布， 研究主要圍繞流場的時變特征展開， 實際燃機中級壓比更大， 因而能造成更強的脈沖射流. Hwang等[3]采用數(shù)值模擬進行了穩(wěn)態(tài)和非穩(wěn)態(tài)情況下氣熱耦合分析， 研究重點關(guān)注了非定常氣熱耦合分析對熱負荷分析的影響及葉片壽命評估的影響. Schmidt等[4]也同樣考察了非定常流動下的氣熱耦合問題， 開發(fā)了一個理論模型用于表征非定常的對流換熱， 并用一個透平靜葉的CFD模擬驗證了該模型， 研究還表明， 非定常壁溫對于時均的熱負荷影響較小. 另一方面， 非定常流動強化換熱的基礎(chǔ)研究也日漸成熟， 通過合理的設(shè)計非定常流動可以在不增加冷氣用量的基礎(chǔ)上強化換熱， 這也將為燃氣輪機設(shè)計水平的進一步提高提供新的思路.

總體而言， 目前關(guān)于燃機內(nèi)冷卻空氣系統(tǒng)的非定常效應(yīng)影響的研究并不多見， 氣熱耦合情況下冷卻系統(tǒng)的非定常現(xiàn)象對于高溫部件溫度分布的影響尚不明確. 本研究針對非定常邊界下的一維導熱問題， 通過求解無量綱方程， 揭示固體溫度在邊界發(fā)生周期性擾動情況下的變化特征， 及其中的關(guān)鍵影響因素. 基于無量綱數(shù)的一維分析結(jié)果， 可以方便地應(yīng)用于工業(yè)設(shè)計, 以快速分析材料及幾何尺寸對問題的影響.

1 數(shù)學模型

圖1 燃氣輪機透平葉片冷卻結(jié)構(gòu)[5] Fig.1 Turbine blade cooling structure ofthe gas turbine

圖2 一維傳熱示意圖Fig.2 Schematic view of a one-dimensional heat transfer

燃氣輪機高溫部件的冷卻， 主要有內(nèi)部對流冷卻和外部氣膜冷卻兩種形式[5]， 典型的渦輪葉片冷卻系統(tǒng)如圖1所示. 考慮沿葉片壁厚方向的傳熱問題， 其包含3個傳熱過程： 主流燃氣側(cè)對流換熱(研究不考慮輻射換熱與氣膜冷卻)、 葉片內(nèi)部固體導熱和冷氣側(cè)對流換熱. 為了簡化問題并嘗試從理論上找出控制各參變量變化的因素， 定性分析各控制因素對數(shù)學模型結(jié)果的影響， 又因為對于局部等溫線與葉片幾何邊界近似平行的區(qū)域， 可認為局部傳熱是近似一維的， 所以研究僅考慮沿葉片壁厚方向的一維導熱， 如圖2所示， 并做如下假設(shè).

1) 固體導熱系數(shù)為常數(shù). 該假設(shè)使得計算結(jié)果中, 固體內(nèi)部溫度時均值的空間分布由非線性變?yōu)榫€性， 而局部溫度波動造成的導熱系數(shù)波動幅值一般不超過2%， 因此該假設(shè)對于溫度幅頻特性的影響可以忽略.

2) 熱流體側(cè)邊界條件為定常狀態(tài). 研究目的是研究冷流體側(cè)的邊界條件波動對傳熱過程的影響， 因此熱流體側(cè)邊界條件設(shè)為定常.

由此建立數(shù)學模型， 考察冷氣側(cè)對流換熱系數(shù)周期性變化對固體內(nèi)溫度分布的影響.

微分方程及邊界條件如式(1)～(3)所示， 式中各變量含義參考文后符號列表[6].

(1)

(2)

其中:

(3)

則該問題的解可表示為t=f(x,τ,a,λ,δ,h1,h2,tf1,tf2).

為了確定該問題的關(guān)鍵影響因素， 對上述方程進行無量綱化， 令

則微分方程及邊界條件可整理為：

(4)

(5)

其中

(6)

2 解的特征

由于一維非穩(wěn)態(tài)導熱方程式(4)關(guān)于Θ是線性的， 則在周期性的邊界條件下， 該方程的解可表示為

(7)

(8)

(9)

采用分離變量法可以對式(9)進行求解， 假設(shè)解的形式如式(10)所示:

Θ′=f(Fo)g(X)

(10)

代入式(9)求解， 整理后可求得

(11)

式中:A、B、C和n均為常量. 特別地， 若給定周期性變化的第一類邊界條件， 則有

(12)

式中:ωFo, n為Θ關(guān)于Fo的波動角速度 (rad·s-1); 系數(shù)Tfw, n與Tbw, n根據(jù)邊界條件確定， 分別影響向X正負方向傳播的波動幅值.

(13)

本文所研究問題即為第三類邊界條件中對流傳熱系數(shù)周期性波動而流體溫度保持常量的情形， 因此雖然方程(4)的解具有式(12)的形式， 但其中系數(shù)A、B無法直接通過邊界條件確定， 因此采用數(shù)值求解的方法進行分析.

此外， 對流傳熱系數(shù)隨時間周期性變化帶來的邊界條件非線性， 也將對傳熱強化或削弱有所貢獻. 為便于理解， 此處以量綱參量表示， 如式(14)所示， 時均的傳熱量并不等于時均的對流傳熱系數(shù)與時均溫差的乘積， 而是多了一項協(xié)方差項， 協(xié)方差項的正負將使得非定常情況下傳熱量大于或小于定常情況. 注意到式(14)中，h和(Tw-Tf)的乘積是兩個周期性信號的乘積， 將產(chǎn)生額外的頻率分量， 亦即， 即使邊界傳熱系數(shù)的波動僅有一個頻率分量， 由此產(chǎn)生的固體溫度波動及傳熱量波動也將包含多個頻率分量.

(14)

3 數(shù)值求解及結(jié)果分析

3.1 數(shù)值求解參數(shù)設(shè)置

由于以燃氣輪機高溫部件冷卻為研究背景， 因此第2節(jié)中數(shù)學模型的相關(guān)系數(shù)， 應(yīng)根據(jù)相應(yīng)的實際情況選取， 以確保各無量綱參數(shù)處于合理范圍. 表1給出了Mark II翼型實驗的相關(guān)參數(shù)[7]， 其中壁厚及對流傳熱系數(shù)的數(shù)據(jù)取自葉片前緣滯止區(qū)， 并忽略固體導熱系數(shù)隨溫度變化， 以此為基礎(chǔ)進行計算. 特別地， 冷氣側(cè)對流傳熱系數(shù)的波動幅值取為時均傳熱系數(shù)值的20%， 波動周期為1 s. 表2給出了數(shù)值計算所需的無量綱系數(shù).

數(shù)值求解采用Matlab軟件提供的pdepe函數(shù)基于線上法(method of lines)對一維拋物型和橢圓型方程進行求解， 即采用近似方法對偏微分方程的空間微分項進行離散， 從而將偏微分方程轉(zhuǎn)化為不再顯含空間獨立變量的常微分方程組， 具體計算方法參考文獻 [8-10]. pdepe函數(shù)計算過程中, 對時間項的積分采用自動時間步長和格式， 對空間項離散采用基于網(wǎng)格尺寸的二階精度格式. 計算中時間長度的選取保證固體內(nèi)部溫度分布已進入穩(wěn)定的周期性變化階段， 空間網(wǎng)格數(shù)的選擇經(jīng)過了網(wǎng)格無關(guān)性測試. 計算結(jié)果表明空間網(wǎng)格數(shù)取128和256時計算結(jié)果幾乎相同， 因此后續(xù)計算分析基于空間網(wǎng)格數(shù)128進行.

表1 Mark II翼型實驗參數(shù)

表2 數(shù)值計算所需無量綱系數(shù)

3.2 固體溫度變化過程及時均值

圖3 邊界溫度變化過程(以時的穩(wěn)態(tài)解為初值)Fig.3 Changing processes of boundary temperature(stationary solution to the as the initial value)

圖4 固體溫度時均值分布Fig.4 Distribution of the time-average of solid temperature

圖5給出了固體兩側(cè)邊界溫度穩(wěn)定周期性波動階段的傅里葉分析結(jié)果. 由圖5可知， 冷氣側(cè)溫度波動幅值的峰值對應(yīng)的頻率與該側(cè)畢渥數(shù)Bi2關(guān)于傅里葉數(shù)Fo的波動頻率fBi2相符， 而在此頻率之外， 仍有一些幅值微小的頻率分量， 其根本原因在于式(4)～(6)構(gòu)成的數(shù)學模型關(guān)于Θ的非線性， 這與圖3中的溫度變化趨勢一致， 也與第2節(jié)中的分析相符. 而固體高溫側(cè)邊界溫度變化的傅里葉分析結(jié)果中， 對應(yīng)fBi2頻率處仍有明顯的峰值， 但其幅值已與其余頻率分量的幅值大致相當， 因此圖5中高溫側(cè)溫度隨時間的變化并不像低溫側(cè)一樣在fBi2頻率下呈現(xiàn)明顯的正弦波動.

圖5 固體邊界溫度波動頻譜分析Fig.5 Spectral analysis of the solid boundary temperature fluctuations

圖6 固體內(nèi)部各處溫度波動幅值的分布 Fig.6 Distribution of temperature fluctuation amplitude in solid

由于固體溫度波動均在fBi2頻率下存在最大幅值， 因此對于溫度波動幅值的分析僅考慮該頻率. 圖6顯示了固體溫度穩(wěn)定周期性變化階段各處溫度波動在fBi2下的幅值分布. 由圖6可知， 固體冷氣側(cè)邊界由于受冷氣對流傳熱系數(shù)周期性變化的直接影響， 其溫度波動幅值最大， 而與冷氣側(cè)距離越遠， 溫度波動幅值越小. 由微分方程(4)的通解形式(12)可知， 溫度波動幅值在固體內(nèi)部以指數(shù)規(guī)律衰減(由于高溫側(cè)近似穩(wěn)態(tài)， 因此僅考慮擾動向X負方向傳播)， 將fBi2代入式(12)， 計算得溫度波動幅值衰減的指數(shù)為6.2， 與圖6擬合結(jié)果較為接近， 考慮到固體內(nèi)溫度波動是多個頻率信號的疊加(如圖5)， 僅由邊界擾動頻率fBi2來估計幅值衰減速率帶來的偏差是可以接受的.

進一步地， 可以根據(jù)溫度波動幅值的衰減定義穿透深度Xp， 參考文獻[11]定義， 將穿透深度定義為溫度波動幅值衰減為邊界溫度波動幅值10%的位置深度. 依此定義， 算例結(jié)果如圖7所示， 穿透深度位置在X=0.625. 同樣的， 根據(jù)式(12)， 可以推導得穿透深度的表達式如式(15)所示(由于高溫側(cè)近似穩(wěn)態(tài)， 因此僅考慮擾動向X負方向傳播)， 將fBi2代入(忽略其余頻率分量)， 可估算得穿透深度位置為Xp=0.628 5， 與數(shù)值計算結(jié)果基本一致.

(15)

穿透深度的意義在于揭示邊界周期性變化的影響范圍， 前述算例的固體壁厚要大于冷氣側(cè)邊界對流傳熱系數(shù)波動的穿透深度， 因此， 可近似認為高溫側(cè)壁溫處于定常狀態(tài). 對于葉片中局部壁厚較薄的地方， 冷氣側(cè)對流傳熱系數(shù)波動對應(yīng)的穿透深度可能大于固體壁厚， 此時高溫側(cè)壁溫也將呈現(xiàn)明顯的周期性變化.

3.3 邊界條件非線性產(chǎn)生的換熱削弱

上文第2節(jié)中分析了對流傳熱系數(shù)的周期性波動造成的邊界條件非線性將對傳熱產(chǎn)生強化或削弱的效果， 如式(14)所示. 圖8給出了數(shù)值計算結(jié)果中固體低溫側(cè)邊界無量綱溫度(Θ)、 無量綱傳熱量(QΘ=Bi2Θ)和畢渥數(shù)(Bi2)的變化過程， 由圖可知三者波動的主頻率是一致的(即fBi2)， 而三者波動存在明顯的相位差. 當邊界對流傳熱系數(shù)增大時， 將增大邊界的傳熱量， 從而導致固體溫度降低， 對流換熱溫差的減小又對換熱量造成影響， 反之亦然. 總體而言， 對流傳熱系數(shù)增大的階段大致對應(yīng)傳熱量增大的階段和固體溫度降低的階段， 由圖8可知，Bi2的峰值時刻、QΘ的峰值時刻和Θ的谷值時刻三者并不嚴格重合， 但較為接近.

圖7 某一時刻各位置下溫度波動幅值與邊界波動幅值的比值Fig.7 Each location at a certain moment under temperature fluctuation amplitude and the ratio of boundary wave amplitude

圖8 固體低溫側(cè)邊界上無量綱溫度、 無量綱傳熱量和畢渥數(shù)隨時間變化規(guī)律Fig.8 Dimensionless temperature, dimensionless heat transfer and Biot number change with time on the boundary of low temperature side

圖9 固體邊界上傳熱量的時均值Fig.9 Time-average of heat transfer on the boundary of low temperature side

同時， 由圖9可知， 穩(wěn)態(tài)情況下， 經(jīng)過固體的熱量要略高于非定常情況， 這與圖4中穩(wěn)態(tài)情況下固體溫度要略低于非定常情況是一致的. 這也意味著， 與穩(wěn)態(tài)情況相比， 冷氣側(cè)對流傳熱系數(shù)的周期性波動削弱了整個系統(tǒng)的換熱， 而這種削弱效應(yīng)是由于固體溫度場在周期性邊界下的響應(yīng)特征造成的， 而與流場細節(jié)無關(guān).

3.4 無量綱參數(shù)分析

圖10 兩側(cè)邊界隨ωFo的變化情況 changes with ωFo on both sides of the boundary

圖11 低溫側(cè)換熱削弱程度隨ωFo的變化情況Fig.11 Extent of weaken heat transfer change with ωFo on the low temperature side

圖12 高溫側(cè)邊界隨畢渥數(shù)Bi1的變化 changes with Bi1 on the high temperature side

圖13 高溫側(cè)換熱削弱程度隨Bi1的變化情況Fig.13 Extent of weaken heat transfer change with Bi1on the high temperature side

圖14 低溫側(cè)邊界隨低溫側(cè)畢渥數(shù)時均值Bi2av的變化 with Bi2av on the low temperature side

圖15 低溫側(cè)換熱削弱程度隨Bi2av的變化情況Fig.15 Extent of weaken heat transfer change with Bi2avon the low temperature side

4 結(jié)語

對第三類邊界發(fā)生周期性波動情況下的固體一維非定常導熱問題進行研究. 首先對數(shù)學模型進行無量綱分析， 推導了第一類邊界下的解析解形式， 而后根據(jù)燃機葉片工作條件進行數(shù)值實驗. 結(jié)果表明：

根據(jù)這一發(fā)現(xiàn)，科學家研制出了“鱟試劑”，這種試劑能準確、快速地檢測出革蘭氏陰性菌的內(nèi)毒素，可以廣泛應(yīng)用于醫(yī)藥衛(wèi)生、食品安全和水質(zhì)檢測等和人類生活息息相關(guān)的領(lǐng)域。

1) 影響固體內(nèi)部無量綱溫度分布的主要因素包括： 無量綱位置X， 無量綱時間Fo， 以及固體兩側(cè)的畢渥數(shù)Bi1和Bi2；

2) 第三類邊界條件的非線性會引起對流傳熱系數(shù)波動時壁面?zhèn)鳠嵯魅酰?其根源在于邊界對流傳熱系數(shù)和溫度變化之間的相位差；

3) 邊界對流傳熱系數(shù)周期性變化對應(yīng)的穿透深度主要受邊界畢渥數(shù)變化頻率fBi2或波動角速度ωFo的影響；

4) 關(guān)于Fo的波動角速度參量ωFo越小， 穿透深度越大， 傳熱削弱越明顯； 擾動側(cè)畢渥數(shù)Bi2存在某一取值， 使得擾動側(cè)無量綱壁溫波動幅值最大.

參考文獻：

[1] ZIMMER V J, RUTLEDGE J L, KNIERIEM C,etal. The influence of coolant unsteadiness on impingement heat transfer[C]//ASME Turbo Expo. Duesseldorf: ASME, 2014.

[2] SCHERHAG C, GEIERMANN J P, WARTZEK F,etal. Blade triggered excitation of periodically unsteady impinging jets for efficient turbine liner segment cooling[J]. Journal of Turbomachinery, 2015, 138(5): S15.

[3] HWANG S, SON C, SEO D,etal. Comparative study on steady and unsteady conjugate heat transfer analysis of a high pressure turbine blade[J]. Applied Thermal Engineering, 2016, 99: 765-775.

[4] SCHMIDT M, STARKE C. Coupled heat-transfer simulations of turbines in consideration of unsteady flows[J]. International Journal of Thermal Sciences, 2015, 96: 305-318.

[5] 韓介勤， 杜達， 艾卡德. 燃氣輪機換熱和冷卻技術(shù)[M]. 程代京， 謝永慧, 譯. 西安: 西安交通大學出版社， 2005.

[6] 楊世銘， 陶文銓. 傳熱學[M]. 4版. 北京: 高等教育出版社， 2006: 41-44.

[7] HYLTON L D, MIHELC M S, TURNER E R,etal. Analytical and experimental evaluation of the heat transfer distribution over the surfaces of turbine vanes[C]//Aas/Division of Dynamical Astronomy Meeting. Indianapolis: Final Report Detroit Diesel Allison, 1983.

[8] SKEEL R D, BERZINS M. A method for the spatial discretization of parabolic equations in one space variable[J]. SIAM Journal on Scientific and Statistical Computing, 1990, 11(1): 1-32.

[9] HOWARD P. Partial differential equations in MATLAB 7.0[M]. Maryland: University of Maryland, 2005.

[10] 董平， 黃洪雁， 馮國泰. 高壓燃氣渦輪徑向內(nèi)冷葉片氣熱耦合的數(shù)值分析[J]. 航空動力學報, 2008, 23(2): 201-207.

[11] BERGMAN T L, LAVINE A S, INCROPERA F P,etal. Fundamentals of heat and mass transfer[M]. New Jersey: John Wiley & Sons, 2011.

[12] JAWORSKI A J, PICCOLO A. Heat transfer processes in parallel-plate heat exchangers of thermoacoustic devices-numerical and experimental approaches[J]. Applied Thermal Engineering, 2012, 42(2): 145-153.

[13] INCROPERA F P, DEWITT D P, BERGMAN T L,etal. Fundamentals of heat and mass transfer [Hardcover][J]. Staff General Research Papers, 2006, 27(1/2): 222.

[14] SCHLICHTING. Boundary layer theory[M]. New York: Springer, 2000: 569-578.

[15] ZUDIN Y B. Analysis of the processes of heat transfer with periodic intensity with allowance for temperature fluctuations in the heat carrier[J]. J Eng Phys Thermophys, 2000, 73(2): 243-247.