田 鋮, 江思雅,2,3, 符 松*

(1.清華大學 航天航空學院,北京 100084; 2.中物院高性能數(shù)值模擬軟件中心,北京 100088; 3.北京應用物理與計算數(shù)學研究所,北京 100088)

1 引 言

航空發(fā)動機是飛機的心臟。現(xiàn)代航空發(fā)動機的核心部件包括壓氣機、燃燒室和渦輪。其中,壓氣機承擔增加空氣壓力的任務,是發(fā)動機熱力循環(huán)必不可少的一部分。我國正大力發(fā)展更先進的航空發(fā)動機,這要求壓氣機具有更好的氣動性能。為實現(xiàn)這一目標,必須深入理解壓氣機內(nèi)部流動機理。然而,由于壓氣機部件高速旋轉,而且內(nèi)部空間狹小,流動實驗測量受到很多限制。計算流體力學CFD(Computational Fluid Dynamics)可以復現(xiàn)壓氣機內(nèi)部流動細節(jié),是壓氣機流動研究的有力工具。

在壓氣機流動模擬中,湍流模式起關鍵作用。目前,工程上最常用的湍流模式是雷諾平均方法RANS(Reynolds-Averaged Navier-Stokes)。然而,壓氣機內(nèi)部流動有強旋轉特性,且包含葉尖泄漏渦和端壁分離渦等復雜渦系。對于這樣復雜的流動,RANS方法在準確性方面存在缺陷。大渦模擬LES(Large Eddy Simulation)和直接數(shù)值模擬DNS(Direct Numerical Simulation)是更準確的湍流模擬方法,但計算成本很高,目前主要應用于雷諾數(shù)較低的流動。為了兼顧湍流模擬的準確度和計算效率,學者們提出了RANS/LES混合方法。

分離渦模擬DES(Detached Eddy Simulation)是一類RANS/LES混合方法,其基本思想是在附著邊界層使用RANS進行模擬,而在大分離區(qū)域使用LES進行模擬。通過混合,DES既克服了RANS對大分離流動模擬不準確的缺陷,又克服了LES在附著邊界層要求網(wǎng)格過密的缺陷。DES首先由Spalart等[1]提出,至今已有很大發(fā)展,主要改進版本包括DDES(Delayed DES)[2]和IDDES(Improved Delayed DES)[3],本文統(tǒng)稱為DES類方法。采用DES類方法能以適中的計算成本獲得較為準確的湍流流場。目前,已有很多學者將DES類方法用于壓氣機流動模擬,表1列出了其中一部分工作。

應用DES類方法模擬壓氣機流動時,需要特別關注數(shù)值耗散的影響。Marty等[5]對比了三階和五階重構格式的計算效果,發(fā)現(xiàn)高階重構格式能顯著提高小尺度相干結構的分辨率,認為這是因為高階格式的數(shù)值耗散更低。高麗敏等[16]總結了DES類方法在葉輪機械流動中的應用,評論道:“DES類方法在葉輪機械中應用時多采用二階或三階精度的格式,雖然在惡劣工況下對宏觀性能的預測能力提高,但仍然存在盲目性?！盚e等[13]在使用DES類方法模擬壓氣機流動時,也強調(diào)了數(shù)值耗散的影響,認為全場使用高耗散數(shù)值格式有利于RANS分支的穩(wěn)定性,但會損害LES分支的湍流解析能力。因此,其推薦使用混合型格式,即在RANS分支使用高耗散的迎風格式,而在LES分支使用低耗散的中心格式。

表1 DES類方法在壓氣機中應用總結

然而,由表1可知,在壓氣機流動DES計算中,高耗散的二階或三階迎風格式仍非常流行,應用低耗散數(shù)值格式尚未形成共識。因此,本文認為該類計算中數(shù)值耗散的作用值得進一步研究。

2 數(shù)值方法

2.1 壓氣機流動求解器

本文用的流動求解程序是自主編寫的UNITs,其準確性在跨聲速軸流壓氣機[4]和離心壓氣機[17]中得到了良好驗證。

流動求解采用有限體積法。粘性項離散采用二階中心格式,時間推進采用雙時間步隱式LU-SGS(Lower-upper Symmetric Gauss-Seidel)格式。湍流模式采用基于SSTk-ω湍流模式的IDDES方法[3]。

N-S(Navier-Stokes)方程對流項的離散格式對數(shù)值耗散影響最大。在有限體積框架下,對流項的離散過程通常分為兩步,一是變量重構,二是通量演化。相應的,對流項離散格式由變量重構格式和近似Riemann求解器這兩部分組成。為了對比不同數(shù)值格式的作用,本文采用了多種不同對流項離散格式,重構格式包括三階MUSCL格式[18]、四階MDCD格式[19]和五階WENO格式[20];近似Riemann求解器包括標準Roe格式[21]和自適應耗散Roe格式[22]。

四階MDCD重構和自適應耗散Roe格式組合成自適應耗散格式,與IDDES方法相配合,能夠較準確模擬湍流,將在下文詳細介紹。

2.2 MDCD重構格式

MDCD(Minimized Dispersion and Controllable Dissipation)是孫振生等[19]提出的先進通量重構格式。MDCD在每個方向使用六個模版點進行變量重構,但只取四階名義精度,從而獲得兩個自由參量。這兩個參量分別稱為色散系數(shù)γdisp和耗散系數(shù)γdiss,各自獨立地控制重構格式的色散和耗散特性。經(jīng)過優(yōu)化,MDCD格式可以得到最優(yōu)色散系數(shù),而耗散系數(shù)大小則根據(jù)算例特點進行人為調(diào)整。

為了捕捉激波,孫振生等[19]將WENO格式[20]的子模版加權思想引入MDCD格式,只是線性權用MDCD色散和耗散系數(shù)來表示。從這個角度來看,MDCD可以視為WENO線性權的優(yōu)化技術。MDCD的變量重構計算式為

(1)

式中W為流動變量,下標i+1/2為網(wǎng)格界面指標,上標L代表界面左側,上標(m)為第m個子模版重構結果,ωm為第m個子模版的非線性權。式(1)只給出了界面左側變量重構的計算式,而界面右側變量重構是完全對稱的。

(2)

(3)

(4)

(5)

式中 下標i-1,i和i+1等代表網(wǎng)格中心指標。式(1)中非線性權ωm計算為

(6)

式中Cm為線性權,ISm為光滑因子,具體表達式見文獻[20]。從式(1～6),MDCD格式和WENO格式完全相同,但其中線性權Cm的取值并不相同。MDCD格式中,線性權的計算式為

(7)

(8)

(9)

(10)

如前所述,γdisp為色散系數(shù),其最優(yōu)值為0.0463783;γdiss為耗散系數(shù),根據(jù)算例特點進行調(diào)整,本文取0.012。

2.3 自適應耗散Roe格式

(11)

自適應耗散函數(shù)φ表達式為

(12)

詳細展開后較為復雜,包含多個開關函數(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)

3 結果與討論

3.1 各向同性衰減湍流

各向同性衰減湍流DIT(Decaying Isotropic Turbulence)是最簡單的湍流流動之一。從系綜平均觀點看,該流動的平均速度處處為零,只有脈動速度非零,而且湍流統(tǒng)計量在空間均勻分布。湍動能k輸運方程的對流項、生成項和擴散項均為零,只余下非定常項和耗散項ε,即

?k/?t=-ε

(14)

因此,DIT非常適合于定量測試數(shù)值耗散。本文用IDDES方法對DIT進行了模擬,并比較了不同數(shù)值格式的效果。值得注意的是,因為DIT不存在壁面,所以全計算域都激活LES分支,而RANS分支并不生效。

本文所用程序主要求解可壓縮流動,故選擇初始湍流馬赫數(shù)為0.3的可壓縮DIT作為測試算例。計算域是立方體,每條邊長度均為2π,網(wǎng)格總數(shù)是64×64×64。所有邊界都設置為周期性邊界。本文計算結果與相應DNS計算結果對比。DNS工作由閆博文[25]完成。

不同數(shù)值格式預測的湍動能隨時間衰減曲線如圖1所示?？梢钥闯?數(shù)值耗散從低到高的排列順序為,(1) 四階中心格式; (2) MDCD重構+自適應耗散Roe格式; (3) MDCD重構+標準Roe格式; (4) 五階WENO重構+標準Roe格式; (5) 三階MUSCL重構+標準Roe格式。定量地來看,在t=1時刻,這五種數(shù)值格式預測湍動能分別是DNS結果的89%,84%,69%,50%和32%。即使是耗散最低的中心格式,IDDES預測的湍動能衰減也快于DNS結果,這種現(xiàn)象產(chǎn)生的原因稍后將結合湍流能譜進行分析。

圖1 采用不同數(shù)值格式模擬DIT,湍動能隨時間的衰減曲線(橫軸表示時間,用大渦翻轉時間無量綱化;縱軸表示湍動能,用初始湍動能無量綱化)

針對t=1時刻的流場,用Fourier變換分析湍動能能譜,并和DNS預測的能譜作比較,結果如圖2所示。圖中橫軸是Fourier波數(shù),用K表示,以便與代表湍動能的k作區(qū)分。圖2反映出來的數(shù)值格式耗散高低排列順序,與圖1一致。在K<5的低波數(shù)區(qū)域,各種數(shù)值格式預測的結果差別不大,而且都與DNS結果吻合較好。然而,在中高波數(shù)區(qū)域,各數(shù)值格式的耗散特性出現(xiàn)明顯區(qū)別。三階MUSCL+標準Roe格式的耗散最強,在K>5時其預測的能譜就遠低于DNS結果。這意味著,三階MUSCL+標準Roe格式只能捕捉大尺度相干結構,且會耗散掉絕大多數(shù)的小尺度結構。采用高階重構后,湍流能譜預測結果得到改善,但在中高波數(shù)區(qū)域的耗散仍然過高。五階WENO重構+標準Roe格式在K>7區(qū)域預測的能譜明顯低于DNS結果,而MDCD重構+標準Roe格式在K>9區(qū)域預測的能譜明顯低于DNS結果。這說明,如果只改變重構格式,而不改變近似Riemann求解器,數(shù)值耗散并不能達到LES分支要求的理想值。

分析四階中心格式預測的湍流能譜。在K<11的區(qū)域,中心格式和DNS結果總體吻合較好,略微偏高。在1123的高波數(shù)區(qū)域,中心格式預測的能譜高于DNS結果,即能譜有翹起的趨勢。這種現(xiàn)象稱為能量累積(energy accumulation),主要由色散誤差和混淆誤差引起。如果把四階中心格式直接用于壓氣機流動模擬,因為流場非常復雜,能譜高波數(shù)區(qū)域會出現(xiàn)更加嚴重的能量累積,引起強烈數(shù)值振蕩,導致計算發(fā)散。所以,對于復雜流動,必須引入適當?shù)臄?shù)值耗散來抑制高波數(shù)區(qū)域的能量累積現(xiàn)象。

MDCD重構配合自適應耗散Roe格式的預測結果最好。在中低波數(shù)區(qū)域,該格式的表現(xiàn)和中心格式類似,基本吻合DNS結果。而且,該格式有適度數(shù)值耗散,抑制了高波數(shù)區(qū)域的能量累積現(xiàn)象。

圖2 t=1時刻不同數(shù)值格式預測的DIT湍動能能譜

總體來說,DIT模擬結果表明,對于IDDES的LES分支,三階MUSCL+標準Roe格式的耗散過高,只能分辨大尺度湍流結構,無法完全發(fā)揮LES對湍流的解析能力。如果僅采用高階的重構格式,如WENO或MDCD重構,可以在一定程度上提高數(shù)值分辨率,但對中、高波數(shù)區(qū)域的能量耗散仍然過高。純中心格式能夠比較準確地預測中低波數(shù)的湍流能譜,但在最高波數(shù)區(qū)域會出現(xiàn)能量累積現(xiàn)象,這很容易引起計算發(fā)散。MDCD重構配合自適應耗散Roe格式可以抑制高波數(shù)的能量累積現(xiàn)象,在全波數(shù)范圍內(nèi)都表現(xiàn)出較好的性能。

3.2 跨聲速離心壓氣機

比較兩種不同數(shù)值格式在跨聲速離心壓氣機上的IDDES計算結果。該跨聲速離心壓氣機由清華大學設計[26],主要參數(shù)列入表2。本文模擬了離心葉輪和無葉擴壓器內(nèi)部流動。表2列出的70500 轉/分是60%設計轉速,也是本文計算采用的轉速。

表2 跨聲速離心壓氣機主要參數(shù)

本文主要討論數(shù)值格式的影響,為提高計算效率,計算域只包括單個葉片通道及其對應的無葉擴壓器,如圖3所示?？偩W(wǎng)格數(shù)約340萬,葉高方向使用91個網(wǎng)格,其中26個網(wǎng)格位于葉尖間隙內(nèi)。絕大多數(shù)壁面上的第一層網(wǎng)格大小滿足y+<1,而且Δx+和Δz+均小于60,這符合IDDES湍流模式的一般要求。時間步長取為轉子旋轉周期的1/1440,即每個葉片通過周期內(nèi)包括60個物理時間步。每個物理時間步內(nèi)進行20次子迭代,計算殘差可降低一個數(shù)量級。非定常計算收斂的判斷依據(jù)是:壓氣機總體性能參數(shù),包括質量流量和靜總壓比,時間平均值基本不再變化。

進口邊界給定總溫288 K,總壓101 kPa,速度取為軸向。出口邊界給定均勻分布的靜壓,取值為145 kPa～205 kPa。不同的出口靜壓值對應壓氣機不同的工作點。環(huán)向取旋轉周期性邊界條件。所有固體壁面都給定絕熱無滑移條件。

圖3 跨聲速離心壓氣機計算域和網(wǎng)格

作為一種RANS/LES混合方法,IDDES通過混合函數(shù)fd來控制兩個分支的生效區(qū)域。采用自適應耗散格式計算時,該混合函數(shù)fd的分布如圖4所示。在靠近壁面的區(qū)域,fd趨近于1,RANS分支生效;在遠離壁面的主流區(qū)域,fd趨近于0,LES分支生效。可以看出,fd分布基本符合IDDES方法的設計預期。但同時,圖4也顯示出其模擬壓氣機流動的一個缺陷,即葉尖間隙區(qū)域內(nèi)fd趨近于1。這說明葉尖間隙流動全部采用了RANS模式計算,這很可能會造成葉尖泄漏渦模擬不準。目前有學者[13]正嘗試改進DES類方法,以解決該問題。

圖4 自適應耗散格式計算時,RANS/LES混合函數(shù)fd分布

圖5 數(shù)值格式的自適應耗散函數(shù)分布

CFD預測的壓氣機性能曲線和實驗測量結果對比如圖6所示,圖中橫軸是質量流量,縱軸是靜總壓比?？梢钥闯?兩種格式的模擬結果均低估了壓比,而自適應耗散格式相比于三階迎風格式,結果與實驗吻合得更好。在流量0.33 kg/s附近,兩種數(shù)值格式結果的差異最大,自適應耗散格式預測壓比與實驗結果誤差為1%,三階迎風格式與實驗結果誤差為7.3%;在流量0.29 kg/s附近,兩種數(shù)值格式差異最小,自適應耗散格式預測壓比與實驗結果誤差為4.6%,三階迎風格式與實驗結果誤差為6.8%。在大多數(shù)工作點,自適應耗散格式預測壓比與實驗結果的相對誤差約為4%～5%;三階迎風格式預測的壓比與實驗結果的相對誤差約為6%～7%。

在DIT模擬中,本文發(fā)現(xiàn)自適應耗散格式可以更準確地預測高波數(shù)能譜,這說明其對湍流小尺度結構的分辨率更高。在跨聲速離心壓氣機計算中,自適應耗散格式同樣在小尺度結構分辨率方面表現(xiàn)出了優(yōu)勢。圖7用熵云圖展示了兩種數(shù)值格式預測的葉片尾跡流動,圖8則用渦識別判據(jù)λ2等值面展示了葉片通道內(nèi)的三維渦結構。圖7和圖8都反映出自適應耗散格式比三階迎風格式捕捉了更多的小尺度相干結構。在圖7中,自適應耗散格式可以捕捉到葉片尾跡的非定常脫落渦,而三階迎風格式預測的尾跡渦幾乎呈現(xiàn)定常狀態(tài)。在圖8中,自適應耗散格式可以捕捉到豐富的小尺度渦結構,而三階迎風格式只能捕捉大尺度渦結構。

圖6 離心壓氣機性能曲線

圖7 兩種不同數(shù)值方法預測的葉片尾跡(用熵渲染顏色)

圖8 兩種不同數(shù)值方法預測的三維渦結構(用λ2等值面識別渦,用熵渲染顏色)

4 結 論

通過對比不同數(shù)值格式的計算結果,本文研究了數(shù)值耗散對各向同性衰減湍流和跨聲速離心壓氣機分離渦模擬DES的影響。主要結論如下。

(1) 在各向同性衰減湍流模擬中,三階迎風格式的耗散過高,嚴重低估了中高波數(shù)的湍流能譜。通過采用高階重構,能夠在一定程度上改善能譜預測結果,但中高波數(shù)的耗散仍然偏高。自適應耗散格式可以在整個波數(shù)范圍內(nèi)準確地預測湍流能譜。

(2) 相比于三階迎風格式,自適應耗散格式更準確地預測了離心壓氣機的性能參數(shù)。具體來說,三階迎風格式預測的壓比與實驗結果的相對誤差約為6%～7%,而自適應耗散格式的預測誤差則降低至4%～5%。

(3) 在跨聲速離心壓氣機的IDDES計算中,三階迎風格式會耗散大部分小尺度相干結構,只能捕捉大尺度渦旋。相比之下,自適應耗散格式顯著提高了小尺度相干結構的分辨率。

綜合來看,在采用DES類方法模擬壓氣機流動時,有必要采用數(shù)值耗散較低的離散格式,以充分發(fā)揮DES類方法的精度優(yōu)勢。由MDCD重構和自適應耗散Roe格式組合成的自適應耗散格式模擬結果比高耗散迎風格式更加準確。

后續(xù)研究的可能方向是改良自適應耗散函數(shù),以使其更加適合于壓氣機流動模擬,進一步提升模擬精度。