田 鋮, 江思雅,2,3, 符 松*
(1.清華大學(xué) 航天航空學(xué)院,北京 100084; 2.中物院高性能數(shù)值模擬軟件中心,北京 100088; 3.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100088)
航空發(fā)動(dòng)機(jī)是飛機(jī)的心臟。現(xiàn)代航空發(fā)動(dòng)機(jī)的核心部件包括壓氣機(jī)、燃燒室和渦輪。其中,壓氣機(jī)承擔(dān)增加空氣壓力的任務(wù),是發(fā)動(dòng)機(jī)熱力循環(huán)必不可少的一部分。我國(guó)正大力發(fā)展更先進(jìn)的航空發(fā)動(dòng)機(jī),這要求壓氣機(jī)具有更好的氣動(dòng)性能。為實(shí)現(xiàn)這一目標(biāo),必須深入理解壓氣機(jī)內(nèi)部流動(dòng)機(jī)理。然而,由于壓氣機(jī)部件高速旋轉(zhuǎn),而且內(nèi)部空間狹小,流動(dòng)實(shí)驗(yàn)測(cè)量受到很多限制。計(jì)算流體力學(xué)CFD(Computational Fluid Dynamics)可以復(fù)現(xiàn)壓氣機(jī)內(nèi)部流動(dòng)細(xì)節(jié),是壓氣機(jī)流動(dòng)研究的有力工具。
在壓氣機(jī)流動(dòng)模擬中,湍流模式起關(guān)鍵作用。目前,工程上最常用的湍流模式是雷諾平均方法RANS(Reynolds-Averaged Navier-Stokes)。然而,壓氣機(jī)內(nèi)部流動(dòng)有強(qiáng)旋轉(zhuǎn)特性,且包含葉尖泄漏渦和端壁分離渦等復(fù)雜渦系。對(duì)于這樣復(fù)雜的流動(dòng),RANS方法在準(zhǔn)確性方面存在缺陷。大渦模擬LES(Large Eddy Simulation)和直接數(shù)值模擬DNS(Direct Numerical Simulation)是更準(zhǔn)確的湍流模擬方法,但計(jì)算成本很高,目前主要應(yīng)用于雷諾數(shù)較低的流動(dòng)。為了兼顧湍流模擬的準(zhǔn)確度和計(jì)算效率,學(xué)者們提出了RANS/LES混合方法。
分離渦模擬DES(Detached Eddy Simulation)是一類RANS/LES混合方法,其基本思想是在附著邊界層使用RANS進(jìn)行模擬,而在大分離區(qū)域使用LES進(jìn)行模擬。通過(guò)混合,DES既克服了RANS對(duì)大分離流動(dòng)模擬不準(zhǔn)確的缺陷,又克服了LES在附著邊界層要求網(wǎng)格過(guò)密的缺陷。DES首先由Spalart等[1]提出,至今已有很大發(fā)展,主要改進(jìn)版本包括DDES(Delayed DES)[2]和IDDES(Improved Delayed DES)[3],本文統(tǒng)稱為DES類方法。采用DES類方法能以適中的計(jì)算成本獲得較為準(zhǔn)確的湍流流場(chǎng)。目前,已有很多學(xué)者將DES類方法用于壓氣機(jī)流動(dòng)模擬,表1列出了其中一部分工作。
應(yīng)用DES類方法模擬壓氣機(jī)流動(dòng)時(shí),需要特別關(guān)注數(shù)值耗散的影響。Marty等[5]對(duì)比了三階和五階重構(gòu)格式的計(jì)算效果,發(fā)現(xiàn)高階重構(gòu)格式能顯著提高小尺度相干結(jié)構(gòu)的分辨率,認(rèn)為這是因?yàn)楦唠A格式的數(shù)值耗散更低。高麗敏等[16]總結(jié)了DES類方法在葉輪機(jī)械流動(dòng)中的應(yīng)用,評(píng)論道:“DES類方法在葉輪機(jī)械中應(yīng)用時(shí)多采用二階或三階精度的格式,雖然在惡劣工況下對(duì)宏觀性能的預(yù)測(cè)能力提高,但仍然存在盲目性?!盚e等[13]在使用DES類方法模擬壓氣機(jī)流動(dòng)時(shí),也強(qiáng)調(diào)了數(shù)值耗散的影響,認(rèn)為全場(chǎng)使用高耗散數(shù)值格式有利于RANS分支的穩(wěn)定性,但會(huì)損害LES分支的湍流解析能力。因此,其推薦使用混合型格式,即在RANS分支使用高耗散的迎風(fēng)格式,而在LES分支使用低耗散的中心格式。
表1 DES類方法在壓氣機(jī)中應(yīng)用總結(jié)
然而,由表1可知,在壓氣機(jī)流動(dòng)DES計(jì)算中,高耗散的二階或三階迎風(fēng)格式仍非常流行,應(yīng)用低耗散數(shù)值格式尚未形成共識(shí)。因此,本文認(rèn)為該類計(jì)算中數(shù)值耗散的作用值得進(jìn)一步研究。
本文用的流動(dòng)求解程序是自主編寫的UNITs,其準(zhǔn)確性在跨聲速軸流壓氣機(jī)[4]和離心壓氣機(jī)[17]中得到了良好驗(yàn)證。
流動(dòng)求解采用有限體積法。粘性項(xiàng)離散采用二階中心格式,時(shí)間推進(jìn)采用雙時(shí)間步隱式LU-SGS(Lower-upper Symmetric Gauss-Seidel)格式。湍流模式采用基于SSTk-ω湍流模式的IDDES方法[3]。
N-S(Navier-Stokes)方程對(duì)流項(xiàng)的離散格式對(duì)數(shù)值耗散影響最大。在有限體積框架下,對(duì)流項(xiàng)的離散過(guò)程通常分為兩步,一是變量重構(gòu),二是通量演化。相應(yīng)的,對(duì)流項(xiàng)離散格式由變量重構(gòu)格式和近似Riemann求解器這兩部分組成。為了對(duì)比不同數(shù)值格式的作用,本文采用了多種不同對(duì)流項(xiàng)離散格式,重構(gòu)格式包括三階MUSCL格式[18]、四階MDCD格式[19]和五階WENO格式[20];近似Riemann求解器包括標(biāo)準(zhǔn)Roe格式[21]和自適應(yīng)耗散Roe格式[22]。
四階MDCD重構(gòu)和自適應(yīng)耗散Roe格式組合成自適應(yīng)耗散格式,與IDDES方法相配合,能夠較準(zhǔn)確模擬湍流,將在下文詳細(xì)介紹。
MDCD(Minimized Dispersion and Controllable Dissipation)是孫振生等[19]提出的先進(jìn)通量重構(gòu)格式。MDCD在每個(gè)方向使用六個(gè)模版點(diǎn)進(jìn)行變量重構(gòu),但只取四階名義精度,從而獲得兩個(gè)自由參量。這兩個(gè)參量分別稱為色散系數(shù)γdisp和耗散系數(shù)γdiss,各自獨(dú)立地控制重構(gòu)格式的色散和耗散特性。經(jīng)過(guò)優(yōu)化,MDCD格式可以得到最優(yōu)色散系數(shù),而耗散系數(shù)大小則根據(jù)算例特點(diǎn)進(jìn)行人為調(diào)整。
為了捕捉激波,孫振生等[19]將WENO格式[20]的子模版加權(quán)思想引入MDCD格式,只是線性權(quán)用MDCD色散和耗散系數(shù)來(lái)表示。從這個(gè)角度來(lái)看,MDCD可以視為WENO線性權(quán)的優(yōu)化技術(shù)。MDCD的變量重構(gòu)計(jì)算式為
(1)
式中W為流動(dòng)變量,下標(biāo)i+1/2為網(wǎng)格界面指標(biāo),上標(biāo)L代表界面左側(cè),上標(biāo)(m)為第m個(gè)子模版重構(gòu)結(jié)果,ωm為第m個(gè)子模版的非線性權(quán)。式(1)只給出了界面左側(cè)變量重構(gòu)的計(jì)算式,而界面右側(cè)變量重構(gòu)是完全對(duì)稱的。
(2)
(3)
(4)
(5)
式中 下標(biāo)i-1,i和i+1等代表網(wǎng)格中心指標(biāo)。式(1)中非線性權(quán)ωm計(jì)算為
(6)
式中Cm為線性權(quán),ISm為光滑因子,具體表達(dá)式見(jiàn)文獻(xiàn)[20]。從式(1~6),MDCD格式和WENO格式完全相同,但其中線性權(quán)Cm的取值并不相同。MDCD格式中,線性權(quán)的計(jì)算式為
(7)
(8)
(9)
(10)
如前所述,γdisp為色散系數(shù),其最優(yōu)值為0.0463783;γdiss為耗散系數(shù),根據(jù)算例特點(diǎn)進(jìn)行調(diào)整,本文取0.012。
(11)
自適應(yīng)耗散函數(shù)φ表達(dá)式為
(12)
詳細(xì)展開后較為復(fù)雜,包含多個(gè)開關(guān)函數(shù),即
A=max{[(lgrid/lturb)/g0-0.5],0}
lgrid=CDESmax(Δx,Δy,Δz)
B0=2Ωmax(Ω,S)/max[(S2+Ω2)/2,10-20]
Ks=max{[(S2+Ω2)/2]1/2,0.1}
CDESε=0.61,CDESω=0.78,Cμ=0.09
(13)
各向同性衰減湍流DIT(Decaying Isotropic Turbulence)是最簡(jiǎn)單的湍流流動(dòng)之一。從系綜平均觀點(diǎn)看,該流動(dòng)的平均速度處處為零,只有脈動(dòng)速度非零,而且湍流統(tǒng)計(jì)量在空間均勻分布。湍動(dòng)能k輸運(yùn)方程的對(duì)流項(xiàng)、生成項(xiàng)和擴(kuò)散項(xiàng)均為零,只余下非定常項(xiàng)和耗散項(xiàng)ε,即
?k/?t=-ε
(14)
因此,DIT非常適合于定量測(cè)試數(shù)值耗散。本文用IDDES方法對(duì)DIT進(jìn)行了模擬,并比較了不同數(shù)值格式的效果。值得注意的是,因?yàn)镈IT不存在壁面,所以全計(jì)算域都激活LES分支,而RANS分支并不生效。
本文所用程序主要求解可壓縮流動(dòng),故選擇初始湍流馬赫數(shù)為0.3的可壓縮DIT作為測(cè)試算例。計(jì)算域是立方體,每條邊長(zhǎng)度均為2π,網(wǎng)格總數(shù)是64×64×64。所有邊界都設(shè)置為周期性邊界。本文計(jì)算結(jié)果與相應(yīng)DNS計(jì)算結(jié)果對(duì)比。DNS工作由閆博文[25]完成。
不同數(shù)值格式預(yù)測(cè)的湍動(dòng)能隨時(shí)間衰減曲線如圖1所示??梢钥闯?數(shù)值耗散從低到高的排列順序?yàn)?(1) 四階中心格式; (2) MDCD重構(gòu)+自適應(yīng)耗散Roe格式; (3) MDCD重構(gòu)+標(biāo)準(zhǔn)Roe格式; (4) 五階WENO重構(gòu)+標(biāo)準(zhǔn)Roe格式; (5) 三階MUSCL重構(gòu)+標(biāo)準(zhǔn)Roe格式。定量地來(lái)看,在t=1時(shí)刻,這五種數(shù)值格式預(yù)測(cè)湍動(dòng)能分別是DNS結(jié)果的89%,84%,69%,50%和32%。即使是耗散最低的中心格式,IDDES預(yù)測(cè)的湍動(dòng)能衰減也快于DNS結(jié)果,這種現(xiàn)象產(chǎn)生的原因稍后將結(jié)合湍流能譜進(jìn)行分析。
圖1 采用不同數(shù)值格式模擬DIT,湍動(dòng)能隨時(shí)間的衰減曲線(橫軸表示時(shí)間,用大渦翻轉(zhuǎn)時(shí)間無(wú)量綱化;縱軸表示湍動(dòng)能,用初始湍動(dòng)能無(wú)量綱化)
針對(duì)t=1時(shí)刻的流場(chǎng),用Fourier變換分析湍動(dòng)能能譜,并和DNS預(yù)測(cè)的能譜作比較,結(jié)果如圖2所示。圖中橫軸是Fourier波數(shù),用K表示,以便與代表湍動(dòng)能的k作區(qū)分。圖2反映出來(lái)的數(shù)值格式耗散高低排列順序,與圖1一致。在K<5的低波數(shù)區(qū)域,各種數(shù)值格式預(yù)測(cè)的結(jié)果差別不大,而且都與DNS結(jié)果吻合較好。然而,在中高波數(shù)區(qū)域,各數(shù)值格式的耗散特性出現(xiàn)明顯區(qū)別。三階MUSCL+標(biāo)準(zhǔn)Roe格式的耗散最強(qiáng),在K>5時(shí)其預(yù)測(cè)的能譜就遠(yuǎn)低于DNS結(jié)果。這意味著,三階MUSCL+標(biāo)準(zhǔn)Roe格式只能捕捉大尺度相干結(jié)構(gòu),且會(huì)耗散掉絕大多數(shù)的小尺度結(jié)構(gòu)。采用高階重構(gòu)后,湍流能譜預(yù)測(cè)結(jié)果得到改善,但在中高波數(shù)區(qū)域的耗散仍然過(guò)高。五階WENO重構(gòu)+標(biāo)準(zhǔn)Roe格式在K>7區(qū)域預(yù)測(cè)的能譜明顯低于DNS結(jié)果,而MDCD重構(gòu)+標(biāo)準(zhǔn)Roe格式在K>9區(qū)域預(yù)測(cè)的能譜明顯低于DNS結(jié)果。這說(shuō)明,如果只改變重構(gòu)格式,而不改變近似Riemann求解器,數(shù)值耗散并不能達(dá)到LES分支要求的理想值。
分析四階中心格式預(yù)測(cè)的湍流能譜。在K<11的區(qū)域,中心格式和DNS結(jié)果總體吻合較好,略微偏高。在11
MDCD重構(gòu)配合自適應(yīng)耗散Roe格式的預(yù)測(cè)結(jié)果最好。在中低波數(shù)區(qū)域,該格式的表現(xiàn)和中心格式類似,基本吻合DNS結(jié)果。而且,該格式有適度數(shù)值耗散,抑制了高波數(shù)區(qū)域的能量累積現(xiàn)象。
圖2 t=1時(shí)刻不同數(shù)值格式預(yù)測(cè)的DIT湍動(dòng)能能譜
總體來(lái)說(shuō),DIT模擬結(jié)果表明,對(duì)于IDDES的LES分支,三階MUSCL+標(biāo)準(zhǔn)Roe格式的耗散過(guò)高,只能分辨大尺度湍流結(jié)構(gòu),無(wú)法完全發(fā)揮LES對(duì)湍流的解析能力。如果僅采用高階的重構(gòu)格式,如WENO或MDCD重構(gòu),可以在一定程度上提高數(shù)值分辨率,但對(duì)中、高波數(shù)區(qū)域的能量耗散仍然過(guò)高。純中心格式能夠比較準(zhǔn)確地預(yù)測(cè)中低波數(shù)的湍流能譜,但在最高波數(shù)區(qū)域會(huì)出現(xiàn)能量累積現(xiàn)象,這很容易引起計(jì)算發(fā)散。MDCD重構(gòu)配合自適應(yīng)耗散Roe格式可以抑制高波數(shù)的能量累積現(xiàn)象,在全波數(shù)范圍內(nèi)都表現(xiàn)出較好的性能。
比較兩種不同數(shù)值格式在跨聲速離心壓氣機(jī)上的IDDES計(jì)算結(jié)果。該跨聲速離心壓氣機(jī)由清華大學(xué)設(shè)計(jì)[26],主要參數(shù)列入表2。本文模擬了離心葉輪和無(wú)葉擴(kuò)壓器內(nèi)部流動(dòng)。表2列出的70500 轉(zhuǎn)/分是60%設(shè)計(jì)轉(zhuǎn)速,也是本文計(jì)算采用的轉(zhuǎn)速。
表2 跨聲速離心壓氣機(jī)主要參數(shù)
本文主要討論數(shù)值格式的影響,為提高計(jì)算效率,計(jì)算域只包括單個(gè)葉片通道及其對(duì)應(yīng)的無(wú)葉擴(kuò)壓器,如圖3所示。總網(wǎng)格數(shù)約340萬(wàn),葉高方向使用91個(gè)網(wǎng)格,其中26個(gè)網(wǎng)格位于葉尖間隙內(nèi)。絕大多數(shù)壁面上的第一層網(wǎng)格大小滿足y+<1,而且Δx+和Δz+均小于60,這符合IDDES湍流模式的一般要求。時(shí)間步長(zhǎng)取為轉(zhuǎn)子旋轉(zhuǎn)周期的1/1440,即每個(gè)葉片通過(guò)周期內(nèi)包括60個(gè)物理時(shí)間步。每個(gè)物理時(shí)間步內(nèi)進(jìn)行20次子迭代,計(jì)算殘差可降低一個(gè)數(shù)量級(jí)。非定常計(jì)算收斂的判斷依據(jù)是:壓氣機(jī)總體性能參數(shù),包括質(zhì)量流量和靜總壓比,時(shí)間平均值基本不再變化。
進(jìn)口邊界給定總溫288 K,總壓101 kPa,速度取為軸向。出口邊界給定均勻分布的靜壓,取值為145 kPa~205 kPa。不同的出口靜壓值對(duì)應(yīng)壓氣機(jī)不同的工作點(diǎn)。環(huán)向取旋轉(zhuǎn)周期性邊界條件。所有固體壁面都給定絕熱無(wú)滑移條件。
圖3 跨聲速離心壓氣機(jī)計(jì)算域和網(wǎng)格
作為一種RANS/LES混合方法,IDDES通過(guò)混合函數(shù)fd來(lái)控制兩個(gè)分支的生效區(qū)域。采用自適應(yīng)耗散格式計(jì)算時(shí),該混合函數(shù)fd的分布如圖4所示。在靠近壁面的區(qū)域,fd趨近于1,RANS分支生效;在遠(yuǎn)離壁面的主流區(qū)域,fd趨近于0,LES分支生效??梢钥闯?fd分布基本符合IDDES方法的設(shè)計(jì)預(yù)期。但同時(shí),圖4也顯示出其模擬壓氣機(jī)流動(dòng)的一個(gè)缺陷,即葉尖間隙區(qū)域內(nèi)fd趨近于1。這說(shuō)明葉尖間隙流動(dòng)全部采用了RANS模式計(jì)算,這很可能會(huì)造成葉尖泄漏渦模擬不準(zhǔn)。目前有學(xué)者[13]正嘗試改進(jìn)DES類方法,以解決該問(wèn)題。
圖4 自適應(yīng)耗散格式計(jì)算時(shí),RANS/LES混合函數(shù)fd分布
圖5 數(shù)值格式的自適應(yīng)耗散函數(shù)分布
CFD預(yù)測(cè)的壓氣機(jī)性能曲線和實(shí)驗(yàn)測(cè)量結(jié)果對(duì)比如圖6所示,圖中橫軸是質(zhì)量流量,縱軸是靜總壓比??梢钥闯?兩種格式的模擬結(jié)果均低估了壓比,而自適應(yīng)耗散格式相比于三階迎風(fēng)格式,結(jié)果與實(shí)驗(yàn)吻合得更好。在流量0.33 kg/s附近,兩種數(shù)值格式結(jié)果的差異最大,自適應(yīng)耗散格式預(yù)測(cè)壓比與實(shí)驗(yàn)結(jié)果誤差為1%,三階迎風(fēng)格式與實(shí)驗(yàn)結(jié)果誤差為7.3%;在流量0.29 kg/s附近,兩種數(shù)值格式差異最小,自適應(yīng)耗散格式預(yù)測(cè)壓比與實(shí)驗(yàn)結(jié)果誤差為4.6%,三階迎風(fēng)格式與實(shí)驗(yàn)結(jié)果誤差為6.8%。在大多數(shù)工作點(diǎn),自適應(yīng)耗散格式預(yù)測(cè)壓比與實(shí)驗(yàn)結(jié)果的相對(duì)誤差約為4%~5%;三階迎風(fēng)格式預(yù)測(cè)的壓比與實(shí)驗(yàn)結(jié)果的相對(duì)誤差約為6%~7%。
在DIT模擬中,本文發(fā)現(xiàn)自適應(yīng)耗散格式可以更準(zhǔn)確地預(yù)測(cè)高波數(shù)能譜,這說(shuō)明其對(duì)湍流小尺度結(jié)構(gòu)的分辨率更高。在跨聲速離心壓氣機(jī)計(jì)算中,自適應(yīng)耗散格式同樣在小尺度結(jié)構(gòu)分辨率方面表現(xiàn)出了優(yōu)勢(shì)。圖7用熵云圖展示了兩種數(shù)值格式預(yù)測(cè)的葉片尾跡流動(dòng),圖8則用渦識(shí)別判據(jù)λ2等值面展示了葉片通道內(nèi)的三維渦結(jié)構(gòu)。圖7和圖8都反映出自適應(yīng)耗散格式比三階迎風(fēng)格式捕捉了更多的小尺度相干結(jié)構(gòu)。在圖7中,自適應(yīng)耗散格式可以捕捉到葉片尾跡的非定常脫落渦,而三階迎風(fēng)格式預(yù)測(cè)的尾跡渦幾乎呈現(xiàn)定常狀態(tài)。在圖8中,自適應(yīng)耗散格式可以捕捉到豐富的小尺度渦結(jié)構(gòu),而三階迎風(fēng)格式只能捕捉大尺度渦結(jié)構(gòu)。
圖6 離心壓氣機(jī)性能曲線
圖7 兩種不同數(shù)值方法預(yù)測(cè)的葉片尾跡(用熵渲染顏色)
圖8 兩種不同數(shù)值方法預(yù)測(cè)的三維渦結(jié)構(gòu)(用λ2等值面識(shí)別渦,用熵渲染顏色)
通過(guò)對(duì)比不同數(shù)值格式的計(jì)算結(jié)果,本文研究了數(shù)值耗散對(duì)各向同性衰減湍流和跨聲速離心壓氣機(jī)分離渦模擬DES的影響。主要結(jié)論如下。
(1) 在各向同性衰減湍流模擬中,三階迎風(fēng)格式的耗散過(guò)高,嚴(yán)重低估了中高波數(shù)的湍流能譜。通過(guò)采用高階重構(gòu),能夠在一定程度上改善能譜預(yù)測(cè)結(jié)果,但中高波數(shù)的耗散仍然偏高。自適應(yīng)耗散格式可以在整個(gè)波數(shù)范圍內(nèi)準(zhǔn)確地預(yù)測(cè)湍流能譜。
(2) 相比于三階迎風(fēng)格式,自適應(yīng)耗散格式更準(zhǔn)確地預(yù)測(cè)了離心壓氣機(jī)的性能參數(shù)。具體來(lái)說(shuō),三階迎風(fēng)格式預(yù)測(cè)的壓比與實(shí)驗(yàn)結(jié)果的相對(duì)誤差約為6%~7%,而自適應(yīng)耗散格式的預(yù)測(cè)誤差則降低至4%~5%。
(3) 在跨聲速離心壓氣機(jī)的IDDES計(jì)算中,三階迎風(fēng)格式會(huì)耗散大部分小尺度相干結(jié)構(gòu),只能捕捉大尺度渦旋。相比之下,自適應(yīng)耗散格式顯著提高了小尺度相干結(jié)構(gòu)的分辨率。
綜合來(lái)看,在采用DES類方法模擬壓氣機(jī)流動(dòng)時(shí),有必要采用數(shù)值耗散較低的離散格式,以充分發(fā)揮DES類方法的精度優(yōu)勢(shì)。由MDCD重構(gòu)和自適應(yīng)耗散Roe格式組合成的自適應(yīng)耗散格式模擬結(jié)果比高耗散迎風(fēng)格式更加準(zhǔn)確。
后續(xù)研究的可能方向是改良自適應(yīng)耗散函數(shù),以使其更加適合于壓氣機(jī)流動(dòng)模擬,進(jìn)一步提升模擬精度。
計(jì)算力學(xué)學(xué)報(bào)2024年1期