基于渦粘性模型的翼型湍流流場熵產(chǎn)率計算

， ， ，

(華中科技大學(xué) 能源與動力工程學(xué)院，武漢 430074)

1 引 言

翼型作為葉片的截面形狀，對流體機(jī)械的工作性能起著至關(guān)重要的影響，因此翼型流場的數(shù)值計算在工程上有重要意義。傳統(tǒng)的翼型氣動性能分析側(cè)重于在熱力學(xué)第一定律的基礎(chǔ)上進(jìn)行氣動力計算以及流動分離特性研究，缺乏對翼型流場中能量轉(zhuǎn)換不可逆特性的研究。如果能夠掌握流場中能量損失的大小及分布，就可以有針對性地改善其性能，從而提高整個流體機(jī)械的工作性能。

熵產(chǎn)表征系統(tǒng)的不可逆性及流動中能量損失的大小，基于熱力學(xué)第二定律的熵分析方法，廣泛應(yīng)用于不可逆現(xiàn)象比較集中的流動與傳熱系統(tǒng)中[1-3]。Alabi等[2]將熵分析引入飛行器一體化設(shè)計中，并采用數(shù)值方法計算了B747-200整個流場的熵產(chǎn)。Li等[3]將熵產(chǎn)作為優(yōu)化目標(biāo)對二維翼型進(jìn)行氣動優(yōu)化設(shè)計，卻只考慮粘性耗散的熵產(chǎn)，未能詳細(xì)分析湍流引起的熵產(chǎn)。Doty等[4,5]采用響應(yīng)面方法對非穩(wěn)態(tài)系統(tǒng)的熵產(chǎn)率進(jìn)行計算。Adeyinka等[6]推導(dǎo)了熵的輸運(yùn)方程。Herwig等[7,8]采用RANS方法對熵的輸運(yùn)方程進(jìn)行了時間平均，并將熵產(chǎn)分為四項。計算湍流熵產(chǎn)時，通常采用兩種方法,文獻(xiàn)[7]用湍動能耗散率代替脈動速度場對熵產(chǎn)進(jìn)行求解，從而建立了湍動能耗散率與熵產(chǎn)的關(guān)系，但壁面處的熵產(chǎn)須采用經(jīng)驗公式單獨求解；文獻(xiàn)[8]采用渦粘性模型，通過湍流粘性系數(shù)來代替脈動速度場對熵產(chǎn)進(jìn)行求解。

本文以NACA0012翼型為研究對象，應(yīng)用RANS渦粘性模型中五種不同的湍流模型，計算不同攻角下翼型流場的熵產(chǎn)率，評估由粘性耗散和湍流耗散所引起的不可逆損失的大小及分布情況。翼型流場計算網(wǎng)格采用代數(shù)插值方法生成，通過雙曲正切函數(shù)控制壁面網(wǎng)格間距?？刂品匠痰碾x散與求解分別采用延遲修正技術(shù)(deferred correction)和強(qiáng)隱式法(SIP)，壓力修正采用非線性的動量插值方法(MIM)。

2 控制方程與數(shù)值方法

直角笛卡爾坐標(biāo)系下，穩(wěn)態(tài)雷諾時均控制方程為

(1)

對流項和擴(kuò)散項的離散采用延遲修正方法[9]。以P為節(jié)點的控制體單元第j個界面上的對流項為Cj，則

(2)

第j個界面上的擴(kuò)散項Dj為

(3)

式中Sj為界面上的面積矢量，其正方向與外法向單位矢量一致。對非正交網(wǎng)格，將Dj分解成正交項和非正交項[10]有

(4)

式中dj為相鄰兩節(jié)點間的距離矢量，第一項為正交項，第二項為非正交項，將非正交項歸入源項，從而避免計算單元頂點上的變量，減少了計算量。

對于得到的五對角代數(shù)方程組采用SIP方法迭代求解。

(5)

(6)

3 湍流流場熵產(chǎn)率計算

根據(jù)能量方程，熵的平衡方程為

(7)

由熱力學(xué)第二定律，進(jìn)一步推導(dǎo)出單位體積流體的熵產(chǎn)率為

(8)

式(8)表明流場的熵產(chǎn)率由傳熱產(chǎn)生的不可逆和流體粘性耗散的不可逆引起，其大小取決于流體的溫度梯度場和速度梯度場。

當(dāng)為湍流流動時，采用RANS方法對式(8)進(jìn)行時均運(yùn)算，有

(9)

根據(jù)Boussinesq的渦粘性假設(shè)，采用渦粘性模型EVMs(Eddy Viscosity Models)將湍流的脈動值和時均值聯(lián)系起來，

(10)

式中 湍流粘性系數(shù)μt和湍流Prandtl數(shù)Prt是空間坐標(biāo)函數(shù)，取決于流動狀態(tài)而不是物性參數(shù)。因此，湍流的流場熵產(chǎn)率可以表示為

(11)

對式(11)進(jìn)行體積分，可得流場總熵產(chǎn)率為

(12)

由上述分析可知，引入Boussinesq假設(shè)后，計算湍流流動的流場熵產(chǎn)率關(guān)鍵在于如何確定湍流粘性系數(shù)μt。當(dāng)選擇不同的湍流模型時，μt與湍流時均參數(shù)聯(lián)系起來的關(guān)系式不同。

4 計算結(jié)果與分析

本文采用文獻(xiàn)[15]的實驗設(shè)置和實驗數(shù)據(jù)，以NACA0012翼型繞流為研究目標(biāo)，采用非正交同位網(wǎng)格下的SIMPLE算法計算雷諾數(shù)為2.88×106時，從0°攻角～16.5°攻角流場的熵產(chǎn)率。

翼型流場C型計算網(wǎng)格采用代數(shù)方法生成，在確定邊界上的節(jié)點分布后，通過代數(shù)插值得到計算區(qū)域內(nèi)的節(jié)點分布。為了控制翼型近壁面的節(jié)點分布，采用雙曲正切關(guān)系式(13)，

(13)

為了保證壁面處第一個節(jié)點布置在對數(shù)分布律成立的范圍內(nèi)，讓控制參數(shù)α滿足給定第一個網(wǎng)格間距的已知條件，即

(14)

式中M為節(jié)點個數(shù)，Δs為第一層網(wǎng)格間距，Δq=(sM-s1)/(M-1)，s1和sM分別為起始點與終止點。

翼型的計算域和計算網(wǎng)格如圖1所示，翼型弦長為c，入口邊界距機(jī)翼后緣12.5c，為一半圓，出口邊界距機(jī)翼后緣21c，為21c×25c的矩形，網(wǎng)格大小為240×60。

本文給定入口速度、溫度和壓力，出口邊界條件按局部單向化處理；出口速度由內(nèi)點外推得到，法向速度滿足總體質(zhì)量守恒條件；壁面邊界為無滑移邊界，當(dāng)選擇k-ε和RNGk-ε湍流模型時，采用壁面函數(shù)法來確定壁面上流速的有效擴(kuò)散系數(shù)。

圖2為NACA0012翼型在0°,10°和15°攻角下，翼型壓力系數(shù)分布與實驗數(shù)據(jù)的對比，可以看出,5種湍流模型的計算結(jié)果基本相同，且都和實驗結(jié)果基本吻合，說明網(wǎng)格劃分與數(shù)值方法可信。

圖3為五種湍流模型在0°攻角～16.5°攻角下的流場熵產(chǎn)率曲線，由于不同湍流模型對于湍流粘性系數(shù)計算不同，故熵產(chǎn)率的計算也存在差異。可以看出，流場熵產(chǎn)率隨著翼型攻角的增大不斷增大。當(dāng)攻角小于14°時，熵產(chǎn)率隨攻角平穩(wěn)增大，但隨著攻角進(jìn)一步增大，流場熵產(chǎn)率顯著增大。

圖4為五種湍流模型在不同攻角下的翼型壁面的阻力系數(shù)，可以看出，不同湍流模型計算得到的阻力系數(shù)變化趨勢與相應(yīng)湍流模型計算得到的熵產(chǎn)率變化趨勢一致。圖5給出不同湍流模型下，翼型流場熵產(chǎn)率與翼型阻力系數(shù)的關(guān)系，可以看出，流場熵產(chǎn)率與阻力系數(shù)存在線性關(guān)系。文獻(xiàn)[3]指出，對于不可壓流動，翼型阻力與流場熵產(chǎn)率存在如下線性關(guān)系，

(15)

進(jìn)一步驗證了流場熵產(chǎn)率計算的正確性。

傳統(tǒng)翼型阻力的數(shù)值計算，通常采用物面壓力、摩擦力積分法、遠(yuǎn)場邊界面積分法和尾跡動量損失積分法[16]。式(15)從熱力學(xué)第二定律能量不可逆的角度提供了一種翼型阻力的新計算方法，通過對流場熵產(chǎn)率的積分，得到翼型阻力。

流場熵產(chǎn)由粘性耗散和湍流耗散兩部分組成。圖6為粘性耗散的熵產(chǎn)率隨攻角的變化曲線，可以看出，隨著攻角增大，粘性耗散的熵產(chǎn)率會呈增大趨勢，但隨著攻角進(jìn)一步增大至一定值時，粘性耗散的熵產(chǎn)率會隨攻角的增加而下降。圖7給出14°攻角下，不同湍流模型計算的翼型流線圖?？梢钥闯?14°攻角下，k-ε、RNGk-ε和S -A湍流模型計算流場仍為附著流，未出現(xiàn)分離渦結(jié)構(gòu)；k-ω和SSTk-ω湍流模型下，翼型尾緣出現(xiàn)明顯的分離渦結(jié)構(gòu)。

圖1 NACA0012翼型流場C型網(wǎng)格

Fig.1 C type mesh of NACA0012 flow field

圖2 不同攻角下壓力系數(shù)分布與實驗值比較

Fig.2 Comparison of the pressure coefficient distribution obtained from calculation and experiment at different angle of attack

圖3 不同攻角下翼型流場熵產(chǎn)率

Fig.3 Entropy generation rate vs.angle of attack

圖4 不同攻角下翼型阻力系數(shù)

Fig.4 Drag coefficient vs.angle of attack

因此，當(dāng)攻角增大到一定程度出現(xiàn)分離渦結(jié)構(gòu)時，由粘性耗散引起的熵產(chǎn)會下降，但是由湍流耗散引起的熵產(chǎn)會顯著上升。這是由于出現(xiàn)流動分離時，流場壁面附近速度梯度會變小，邊界層變成紊流，使得粘性耗散的熵產(chǎn)降低，湍流耗散的熵產(chǎn)顯著增加。因此，可以通過粘性耗散的熵產(chǎn)率曲線判定翼型產(chǎn)生分離渦時的攻角。

圖8為10°攻角下，不同湍流模型計算的流場熵產(chǎn)率云圖?？梢钥闯?，翼型流場的熵產(chǎn)主要發(fā)生在翼型前緣、近壁面邊界層處以及尾流處。由式(11)可知，影響流場熵產(chǎn)率的主要因素是湍流粘性系數(shù)與速度梯度。翼型前緣的熵產(chǎn)主要由氣流折轉(zhuǎn)的速度梯度產(chǎn)生的阻力造成；近壁面邊界層的熵產(chǎn)要遠(yuǎn)高于主流區(qū)，其本質(zhì)是一種摩擦損失；湍尾流場由于尾流的高湍流度導(dǎo)致較大的湍流粘性系數(shù)。

圖5 翼型熵產(chǎn)率與阻力系數(shù)關(guān)系

Fig.5 Drag coefficient vs.entropy generation rate

圖6 不同攻角下翼型流場粘性耗散的熵產(chǎn)率

Fig.6 Viscous entropy generation rate vs.attack angle

圖7 14°攻角下不同湍流模型流線圖

Fig.7 Streamlines for attack angle 14°at different turbulence models

圖8 10°攻角下不同湍流模型翼型熵產(chǎn)率云圖

Fig.8 Contours of entropy generation rate for attack angle 10° for different turbulence models

5 結(jié) 論

本文基于非正交同位網(wǎng)格下的SIMPLE算法，根據(jù)渦粘性模型，對比分析了5種不同的湍流模型下翼型的流場熵產(chǎn)率，得到如下結(jié)論。

(1) 采用代數(shù)方法生成翼型流場的C型網(wǎng)格，網(wǎng)格生成速度快，質(zhì)量便于控制；基于非正交網(wǎng)格下的SIMPLE算法，采用延遲修正技術(shù)和強(qiáng)隱式迭代方法，具有變量存儲量小和計算效率高的優(yōu)點。

(2) 翼型流場熵產(chǎn)率由粘性耗散、湍流耗散、溫差傳熱和湍流溫度脈動四部分組成，其大小取決于流體的溫度梯度場和速度梯度場，其中湍流耗散是引起熵產(chǎn)的主要因素。翼型流場的熵產(chǎn)主要發(fā)生在翼型前緣區(qū)、近壁面區(qū)和翼型尾流區(qū)域。

(3) 驗證了翼型流場的熵產(chǎn)率與翼型阻力系數(shù)為線性相關(guān)，不同的湍流模型計算流場熵產(chǎn)率由于湍流粘性系數(shù)計算的不同，存在差異?；跓崃W(xué)第二定律的熵分析方法提供了一種翼型阻力的新計算方法。

(4) 當(dāng)出現(xiàn)分離渦結(jié)構(gòu)時，粘性耗散的熵產(chǎn)率下降，湍流耗散的熵產(chǎn)率顯著上升。粘性耗散的熵產(chǎn)率曲線可以用來判定翼型初始產(chǎn)生分離渦時的攻角。

:

[1] Baker M，Riggins D，Huang J，et al.Second-law methods for the analysis & design of hypersonic vehicles[A].39t hAIAA Thermophysics Conference[C].2007.

[2] Alabi K，Ladeinde F，von Spakovsky M，et al.Asses-sing CFD modeling of entropy generation for the air frame subsystem in an integrated aircraft design/synthesis procedure[A].AIAA Aerospace Sciences Meeting and Exhibit[C].2006.

[3] Li H，Stewart J，F(xiàn)igliola R.Exergy based design methodology for airfoil shape optimization and wing analysis [A].25t hInternational Congress of the Aeronautical Science [C].2006.

[4] Doty J，Camberos J，Yerkes K.Approximate approach for direct calculation of unsteady entropy generation rate for engineering applications[A].50t hAIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition[C].2012.

[5] Doty J，Camberos J，Yerkes K.Comparison of entropy generation rate calculation methods for engineering applications[A].43r dAIAA Thermophysics Confe -rence[C].2012.

[6] Adeyinka O B，Naterer G F.Modeling of entropy production in turbulent flows[J].JournalofFluidsEngineering，2004，126(6):893-899.

[7] Herwig H，Kock F.Local entropy production in turbulent shear flows: A tool for evaluating heat transfer performance[J].JournalofThermalScience，2006，15(2):159-167.

[8] Moore J，Moore J G.Entropy production rates from viscous flow calculations:Part I—A turbulent boun-dary layer flow[A].ASME 1983 International Gas Turbine Conference and Exhibit[C].1983.

[9] Khosla P K，Rubin S G.A diagonally dominant second-order accurate implicit scheme [J].Computers&Fluids，1974，2(2):207-209.

[10] Mathur S R，Murthy J Y.A pressure-based method for unstructured meshes[J].NumericalHeatTransfer,1997,31(2):195-215.

[11] Stone H L.Iterative solution of implicit appro-ximations of multidimensional partial differential equations[J].SIAMJournalonNumericalAnalysis,1968,5(3):530-558.

[12] Rhie C M，Chow W L.Numerical study of the turbulent flow past an airfoil with trailing edge separation[J].AIAAJournal，1983，21(11):1525-1532.

[13] 劉 波，李忠媛，張 濤.一種基于三角形非結(jié)構(gòu)化網(wǎng)格SIMPLE算法的程序設(shè)計[J].計算力學(xué)學(xué)報，2015，32(6):813-819.(LIU B o，LI Zhong-yuan，ZHANG Tao.A computer programming of SIMPLE algorithm based on triangle unstructured grid[J].ChineseJournalofComputationalMechanics，2015，32(6):813-819.(in Chinese))

[14] 柏 威，鄂學(xué)全.基于非結(jié)構(gòu)化同位網(wǎng)格的SIMPLE算法[J].計算力學(xué)學(xué)報，2003，20(6):702-710.(BAI Wei，E Xue -quan.Implement of SIMPLE algorithm based on unstructured colocated meshes[J].ChineseJournalofComputationalMechanics，2003，20(6):702-710.(in Chinese))

[15] Gregory N，Reilly.Low-speed aerodynamic characte -ristics of NACA0012 aerofoil section，including the effects of upper-surface roughness simulating hoar frost[J].Cheminform，1970，23(48):6697-6700.

[16] Vinh H，van Dam C P，Yen D，et al.Drag prediction algorithms for Navier-Stokes solutions about airfoils[A].13t hAIAA Applied Aerodynamics Conference[C].1995.