鐘藝凱 楊金廣 張 敏 楊 帥 劉 艷 王春雪

（1.大連理工大學能源與動力學院；2.北京動力機械研究所）

0 引言

隨著計算機技術的快速發(fā)展和數(shù)值計算方法的日益完善，計算流體力學（CFD）技術在葉輪機械中的有著越來越廣泛的應用，其低廉的成本、高效的性能，使葉輪機械在設計研究過程中減少了大量的實驗成本[1]。如宋國興等[2]對軸流壓氣機進氣旋流畸變進行仿真研究，陳振毅[3]等采用數(shù)值模擬方法研究了軸流壓氣機近失速工況下軸向間隙對徑向流的影響，都獲得了準確的結(jié)果，由此可見數(shù)值模擬在葉輪機械中的應用是準確可靠的。對于葉輪機械內(nèi)部流動的數(shù)值模擬，主要是采用求解雷諾平均N-S方程（RANS）對其進行數(shù)值模擬。由于雷諾應力項的存在使得RANS方程不封閉，科學家根據(jù)湍流運動的基本規(guī)律發(fā)展了使方程組封閉的湍流模型，促進了湍流理論應用的發(fā)展。自1970年以來，湍流模型的研究得到了快速發(fā)展，先后建立了零方程、一方程、兩方程及雷諾應力模型等湍流模型。雖然這些湍流模型能很好地捕捉到一些流動現(xiàn)象，但是其適用性都有一定的限制，其在葉輪機復雜流動數(shù)值分析中的適用性、精度以及穩(wěn)定性等尚無定論，有必要進一步研究[4]。

本文以Durhum靜子葉柵和NASA rotor37壓氣機轉(zhuǎn)子為研究對象，基于所開發(fā)的CFL3D_Turbo數(shù)值計算平臺，研究了B-L，S-A，Chien k-ε，SST k-ω和k-ε-Rt五種湍流模型計算的結(jié)果，重點考察了不同湍流模型的數(shù)值計算精度、流場細節(jié)捕捉程度以及湍流模型計算穩(wěn)定性的影響。

1 求解器開發(fā)

本文求解器基于經(jīng)過廣泛驗證的CFL3D開發(fā)而成，所采用的控制方程組為直角絕對坐標系下的三維雷諾平均Navier-Stokes方程組，其控制方程[5-6]可以表示為：

式中，t為時間；?為守恒變量；為無粘通量；為粘性通量。

控制方程組空間離散格式采用通量差分分裂的Roe[7]格式進行離散，時間項采用近似因子分解法進行迭代求解，采用當?shù)貢r間步長、多重網(wǎng)格以及殘差光順技術加速收斂。

本文對CFL3D經(jīng)過適應性改進，形成專用于葉輪機械計算的CFL3D_Turbo數(shù)值計算平臺。主要改進包括：1）開發(fā)了匹配的前處理模塊。目前采用的策略是讀取商用軟件的網(wǎng)格生成結(jié)果，這里選用NUMECA軟件包中Autogrid5模塊。通過讀取Autogrid5生成的網(wǎng)格文件，對其進行邊界條件分析和坐標系轉(zhuǎn)化，另存為CFL3D_Turbo識別的網(wǎng)格文件格式；2）改進了原求解器。主要有增加旋轉(zhuǎn)域模擬能力，添加適合于圓柱內(nèi)流計算的邊界條件等；3）開發(fā)了較為完備的后處理軟件。

2 湍流模型

2.1 B-L模型

B-L模型[8]是在C-S模型的基礎上改進而得到的一種代數(shù)模型，它不僅適用于邊界層的計算，還能應用于N-S方程中，因其應用簡單且有相當?shù)木?，所以其在工程領域中應用廣泛[9]。下面給出其湍流粘性系數(shù)計算式：

其中，y為到壁面法向距離；ycrossover為內(nèi)、外層粘性相等時y的最小值。其中內(nèi)層粘性計算公式為：

式中，l為混合長度；ω為旋度。

外層粘性計算公式為：

式中，K為Clauser常數(shù)；CCP為附加常數(shù)。

2.2 S-A模型

S-A模型[10]是一種基于經(jīng)驗和量綱分析，建立了求解湍流粘性的輸運方程，主要是針對簡單流動而逐步發(fā)展起來的湍流模型，它對壓力梯度邊界層也有較好的預測。其輸運方程如下：

式中，G?為湍流粘度生成項；Y?為湍流粘度耗散項；S?為源項。湍流動力粘性系數(shù)的計算公式為：

2.3 Chien k-ε模型

Chienk-ε模型[11]是一種低雷諾數(shù)模型，將泰勒級數(shù)展開技術用來處理固體壁附近的動能及其耗散率，同時它具有良好的經(jīng)濟性和計算精度，主要應用于各種管道流動和邊界層流動。其輸運方程為：

式中，P為湍流生成項；C1，C2，f1和f2為常數(shù)。湍流動力粘性系數(shù)計算公式為：

2.4 SST k-ω模型

SST k-ω模型[12]是Menter在標準k-ω模型的基礎上進行改進和發(fā)展出來的，它考慮了近壁區(qū)逆壓梯度邊界層中的主剪切應力的傳遞，并將Bradshaw提出的主剪切應力與湍流動能成比例的假設引入到渦粘性的定義中，這使得它在近壁區(qū)逆壓梯度和分離流動的計算有更好的預測精度，能夠有效的預測逆壓梯度條件下的流體分離的開始點和分離區(qū)的大小。其輸運方程為：

式中右側(cè)前三項分別為生成項、耗散項以及擴散項，ω方程中第四項為交叉擴散項。湍流動力粘性系數(shù)的計算公式為：

2.5 k-ε-Rt模型

k-ε-Rt模型[13]是由湍動能k方程、湍流耗散率ε方程以及無阻尼渦粘度Rt方程三個方程共同組成。這個模型在湍流耗散率ε方程中加入了一個額外的源項，旨在提高非平衡流區(qū)域的耗散率ε水平，這樣可以降低動能和長度尺度以此來改善對逆壓梯度流的預測，其具有很強的數(shù)值魯棒性且易于使用。下面給出它的輸運方程：

式中，P為湍流生成項；E為耗散率ε方程中的附加源項；Cε1，Cε2，Cε3，C1，C2，C3，f1，f2為常數(shù)。其中湍流粘性系數(shù)計算公式為：

3 算例與分析

3.1 Durham壓氣機葉柵

Durham壓氣機葉柵為低速壓氣機平面葉柵，其采用可控擴散葉片。在模擬中進口總溫為293.15K，進口總壓為107 000Pa，進口氣流角為37°；出口給定輪轂處靜壓97 000Pa，其余位置處靜壓由簡單徑向平衡方程得到。本算例計算網(wǎng)格如下圖1所示。

圖1 Durham葉柵網(wǎng)格Fig.1 Durham cascade grid

首先對比CFL3D_Turbo和NUMECA都采用S-A模型計算的出口處參數(shù)分布，證明二者得到了相同的結(jié)果，如圖2所示，確認了求解器的正確性。出口總壓的吻合證明二者計算得到的葉柵損失特性基本相同；而靜壓分布的一致性則表明氣流偏轉(zhuǎn)能力預測的正確性。

圖2 出口邊界條件對比Fig.2 Comparison of outlet boundary conditions

1）流場分析

將B-L模型，S-A模型，Chien k-ε模型，SST k-ω模型和k-ε-Rt模型湍流模型應用到Durham壓氣機葉柵數(shù)值計算中，得到了馬赫數(shù)、靜壓、溫度等在流場中的分布情況。由于缺少Durham壓氣機葉柵實驗流場數(shù)據(jù)，所以將各模型計算得到的流場與NUMECA計算結(jié)果進行對比。圖3給出Durham壓氣機葉柵50%葉高處S1流面的馬赫數(shù)分布，然后對各湍流模型計算得到的流場進行分析和對比，可以看出，氣流以較高的速度進入通道中且進口到葉柵前緣部分馬赫數(shù)分布均勻，五種模型對這一區(qū)域的馬赫數(shù)預測基本一致，且與NUMECA計算結(jié)果相同。當氣流流經(jīng)尾緣處，出現(xiàn)流動分離，形成低速區(qū)，產(chǎn)生尾跡。從計尾跡區(qū)域可以看出，在數(shù)值上基本一致，k-ε-Rt模型計算得到的尾跡區(qū)域偏大，與NUMECA的結(jié)果較為吻合。

圖3 50%葉高S1流面馬赫數(shù)云圖Fig.3 Mach number contour at 50%span

2）壓力系數(shù)分布

為了定量對比不同湍流模型對葉輪機械數(shù)值模擬結(jié)果的影響，將各湍流模型計算得到的Durham壓氣機葉柵表面的靜壓力進行處理，得到葉柵表面壓力系數(shù)及其分布，將處理后的結(jié)果與NUMECA的計算結(jié)果和實驗數(shù)據(jù)進行對比，葉片表面壓力系數(shù)計算公式為：

式中，p為葉片表面靜壓；p01為進口中截面處總壓；p2為出口中截面處靜壓。

從圖4中可以看出，本文所計算得到的結(jié)果與實驗數(shù)值以及NUMECA計算得到的結(jié)果分布趨勢基本一致，數(shù)值吻合度也很高。在20%葉高處靠近前緣的壓力側(cè)，各湍流模型與實驗數(shù)值的分布差距較大且過低預測了前緣處的壓力系數(shù)，而在吸力側(cè)SST k-ω和k-ε-Rt湍流模型的模擬數(shù)值與實驗數(shù)值完全吻合。同時在壓力側(cè)SST k-ω和k-ε-Rt模型得到的模擬數(shù)值與實驗數(shù)值吻合度很高，其對壓力側(cè)的流動預測具有很高的精度；而在吸力側(cè)各模型的模擬數(shù)值與實驗數(shù)值分布也是基本一致的，其中k-ε-Rt模型計算得到的結(jié)果與實驗數(shù)值有很好的吻合度，在葉柵后半段的壓力系數(shù)分布中，Chien k-ε和SST k-ω模型與實驗數(shù)值有著很好的吻合度。在50%葉高處，各湍流模型計算得到的壓力系數(shù)分布與實驗數(shù)值分布趨勢是一致的，k-ε-Rt模型預測的表面壓力系數(shù)在壓力側(cè)和吸力側(cè)與實驗數(shù)值都最為吻合，而在壓力側(cè)靠近前緣位置處S-A模型計算得到的結(jié)果與實驗數(shù)值偏差最大，在壓力側(cè)和吸力側(cè)后半段Chien k-ε模型預測的數(shù)值與實驗數(shù)值誤差相對較大。在70%葉高處，各模型在葉柵前緣處的壓力系數(shù)預測相對較差，而在壓力側(cè)B-L，S-A以及k-ε-Rt模型所預測的數(shù)值與實驗數(shù)值十分吻合；Chien k-ε模型計算得到的葉柵表面壓力系數(shù)分布誤差相對較大。在90%葉高處，各模型預測的壓力系數(shù)值與實驗數(shù)值的分布趨勢基本一致，S-A模型對壓力側(cè)的計算結(jié)果相較于其它四種模型更接近于實驗數(shù)值，而k-ε-Rt模型在吸力側(cè)預測的數(shù)值與實驗數(shù)值更為吻合，但是各模型對靠近前緣處的靜壓系數(shù)預測都偏小。在95%葉高處，各模型計算得到的壓力系數(shù)與實驗數(shù)值的分布趨勢是基本一致的，但是數(shù)值上有一定的誤差，同時從圖中可以看出NUMECA計算得到的結(jié)果與實驗數(shù)值更為吻合，而SST k-ω和k-ε-Rt模型得到結(jié)果的誤差相比大一些，這可能是由于未給定進口端壁邊界層的緣故。

圖4 表面壓力系數(shù)對比Fig.4 Comparison of surface pressure coefficient

圖5給出了5種湍流模型的收斂史對比。由圖可以看出，在300步之前，Chien k-ε模型的收斂曲線振蕩較大，而其余四種收斂過程相對平穩(wěn)；而在300～500步之間，S-A，SST k-ω和B-L模型殘差下降速度更快；在500～1 000步之間，S-A和SST k-ω模型的殘差收斂效果更好，計算穩(wěn)定性更好，B-L和k-ε-Rt模型的計算穩(wěn)定性相較于Chien k-ε模型稍好，但相較于S-A模型殘差收斂效果較差，而Chien k-ε模型收斂曲線振蕩較大，計算穩(wěn)定性和殘差收斂相對較差。

圖5 Durham壓氣機葉柵收斂史Fig.5 Convergence history of the Durham compressor cascade case

3.2 NASA rotor37壓氣機轉(zhuǎn)子

NASA rotor37軸流壓氣機轉(zhuǎn)子的進口處于超聲速狀態(tài)，其轉(zhuǎn)速為17188.7r/min，葉頂周向速度為454.14m/s，進口總壓為101 325.0Pa，進口總溫為288.1K，葉片數(shù)為36。計算網(wǎng)格如下圖6所示。計算考慮了高度為0.356mm的葉尖間隙。

圖6 Rotor37網(wǎng)格Fig.6 Grid used in the Rotor 37 grid

3.2.1 結(jié)果及分析

1）網(wǎng)格無關性驗證

網(wǎng)格疏密程度是數(shù)值模擬中一個很重要的問題，因為不同的網(wǎng)格數(shù)量可能會對計算結(jié)果和收斂速度造成一定的影響。為了既保證數(shù)值計算結(jié)果的精度，又能減小計算量、加快收斂速度，因此需要進行網(wǎng)格無關性驗證的工作。選擇不同的網(wǎng)格數(shù)量對NASA rotor37進行數(shù)值模擬，通過對計算結(jié)果的分析和比較，找到合適的網(wǎng)格數(shù)量。

本文采用四種不同的數(shù)量網(wǎng)格進行網(wǎng)格無關性研究，分別是方案1網(wǎng)格點的數(shù)量為35萬、方案2網(wǎng)格點的數(shù)量為45萬、方案3網(wǎng)格點的數(shù)量為55萬和方案4網(wǎng)格點的數(shù)量為65萬，將四種網(wǎng)格分別進行數(shù)值計算，采用S-A湍流模型對方程組進行封閉求解，對計算結(jié)果進行處理，將不同網(wǎng)格的總體性能與實驗數(shù)值進行對比。

從圖7中可以看出，隨著計算網(wǎng)格數(shù)量的增加，計算得到的總壓比和絕熱效率都會有所提高，當網(wǎng)格點的數(shù)量大于55萬時，可以看到總壓比和絕熱效率都不會發(fā)生明顯的變化，因此認為網(wǎng)格的數(shù)目達到了網(wǎng)格無關性的要求。同時為了減少計算所需時間，提高效率，因此選擇網(wǎng)格數(shù)量為55萬的方案3來進行后續(xù)的數(shù)值計算。

圖7 網(wǎng)格無關性驗證Fig.7 Grid independence verification

2）近最高效率點分析

為了研究湍流模型在NASA rotor37數(shù)值模擬中對流場特性的影響，在近最高效率點，從計算結(jié)果中分別提取了50%，70%和95%葉高處的S1流面的相對馬赫數(shù)分布圖，并將其與實驗數(shù)據(jù)進行對比，從對流場信息的捕捉能力來分析不同湍流模型的模擬能力。

從圖8中可以看出，在50%葉高處，氣流以超聲速狀態(tài)進入通道，在到達前緣處的相對馬赫數(shù)為1.3左右，這與實驗數(shù)值基本一致，對進口處產(chǎn)生的激波位置和強度的預測與實驗數(shù)值基本一致，但對激波的寬度計算相較于實驗數(shù)值偏大。而在通道中間位置，氣流由超聲速變成了亞聲速狀態(tài)，產(chǎn)生了通道內(nèi)的激波，B-L和Chien k-ε模型的計算得到的激波位置、強度以及寬度較實驗數(shù)值偏大，其余湍流模型對激波位置的預測較為準確，S-A模型計算得到的激波強度相較于實驗數(shù)值偏大，而SST k-ω和k-ε-Rt模型計算得到的激波強度和寬度與實驗數(shù)值基本一致。當氣流流經(jīng)吸力側(cè)尾部時，產(chǎn)生了氣流分離現(xiàn)象，形成了尾跡區(qū)域，但各模型的預測結(jié)果與實驗數(shù)值相比偏小，同時對尾緣流動分離現(xiàn)象的捕捉也不夠明顯。

圖8 50%葉高S1流面相對馬赫數(shù)云圖Fig.8 Relative Mach number contour at 50%span

從圖9中可以看出，在70%葉高處，各模型對進口處激波位置的預測基本準確，S-A、SST k-ω和k-ε-Rt模型計算的激波寬度和強度與實驗數(shù)值更為吻合。而對于通道內(nèi)激波的預測，B-L和Chien k-ε模型的結(jié)果與實驗數(shù)值相差較大，激波位置和強度的預測不夠準確，其余三種模型對激波位置的預測與實驗數(shù)據(jù)基本一致，k-ε-Rt和SST k-ω模型計算的激波強度和寬度與實驗數(shù)值基本一致。

圖9 70%葉高S1流面相對馬赫數(shù)云圖Fig.9 Relative Mach number contour at 70%span

從圖10中可以看出，在95%葉高處，氣流以超聲速狀態(tài)進入進口段，到達葉片前緣處產(chǎn)生了一道激波，SA，SST k-ω和k-ε-Rt模型計算得到的激波寬度和強度與實驗數(shù)值相對來說更為接近；而對于通道內(nèi)激波的預測，S-A模型的預測結(jié)果與實驗數(shù)值相差較大，激波位置的預測不準，而其余三種模型對激波位置的預測與實驗數(shù)據(jù)基本一致，同時各湍流模型對于尾跡的預測與實驗數(shù)值大致相同。

圖10 95%葉高S1流面相對馬赫數(shù)云圖Fig.10 Relative Mach number contour at 95%span

為了進一步評估湍流模型的精度，在近最高效率點，對徑向參數(shù)分布進行對比，如圖11～圖14所示。從總壓比對比圖（圖11）中可以看出，各湍流模型計算得到總壓比分布趨勢與實驗數(shù)值基本一致。在0%～20%葉高之間，B-L模型計算過低預測了總壓比，其它四種湍流模型計算出來的數(shù)值較實驗數(shù)值偏大，而Chien kε模型的結(jié)果相較于其它模型與實驗數(shù)值吻合度更高，SST k-ω模型得到的結(jié)果誤差最大為3.7%；在20%～40%葉高之間B-L模型計算得到的總壓比與實驗數(shù)值幾乎完全吻合，而k-ε-Rt模型誤差相對較大，誤差約為1.8%；在40%～80%，B-L模型得到的結(jié)果與實驗數(shù)值誤差最大為1.92%；在80%～100%葉高處，各種模型的預測結(jié)果大致相當，但是在97.1%葉高處，所有模型都過高預測了總壓比，可能是受到機匣邊界層的影響造成的。

圖11 總壓比的比較Fig.11 Comparison of total pressure ratio

從總溫比對比圖（圖12）中可以看出，在0～20%葉高之間，Chien k-ε和SST k-ω模型計算得到的總溫比與實驗數(shù)值更為接近，其余三種模型計算精度相對較差；在20%～100%葉高之間，k-ε-Rt模型計算的總溫比無論是數(shù)值上還是分布趨勢都與實驗數(shù)值具有很高的吻合度，其計算精度最高，而在20%～50%葉高之間，Chien kε模型的計算結(jié)果與實驗數(shù)值誤差最大為1.27%，在70%～100%葉高之間，B-L模型相較于其余四種模型數(shù)值誤差更大，約為1.26%。

圖12 總溫比的比較Fig.12 Comparison of total temperature ratio

從絕熱效率對比圖（圖13）中可以看出，在0～10%葉高之間，所有湍流模型的計算結(jié)果都比實驗數(shù)值偏小，其中k-ε-Rt模型與實驗數(shù)值相對更吻合，而SST kω模型的預測值誤差最大為4%；在10%～30%葉高之間，S-A模型得到的結(jié)果與實驗數(shù)值吻合度更高，而Chien k-ε模型的數(shù)值誤差最大為2.2%；在30%～60%葉高之間，所有模型得到的結(jié)果與實驗數(shù)值的分布趨勢不太一致，而S-A模型計算得到的結(jié)果偏差最大為2.83%；在60%～100%葉高之間，S-A模型的計算得到的結(jié)果和分布趨勢相較于其它模型與實驗數(shù)值更為接近，SST k-ω模型的預測結(jié)果與實驗數(shù)值偏差最大為6.5%。從軸向絕對速度對比圖（圖14）中可以看出，所有湍流模型計算得到的進口軸向絕對速度與實驗數(shù)值都有很高的吻合度，在靠近輪轂，B-L模型的數(shù)值偏差較大為4.5%，在靠近機匣處，SST k-ω模型的數(shù)值偏差相對較大為3.9%。

圖13 絕熱效率的對比Fig.13 Comparison of adiabatic efficiency

圖14 軸向速度對比Fig.14 Axial velocity comparison

圖15為近最高效率計算的收斂史。從圖中可以看出，在450步之前，k-ε-Rt模型出現(xiàn)了五次殘差較為劇烈的跳躍，而其余四種模型的收斂曲線都較為光滑，收斂速度也基本相同；而在450～1 200步之間，S-A模型的收斂曲線最為平穩(wěn)，殘差下降速度也相對較快，而其它四種模型收斂曲線都有較大的振蕩，數(shù)值穩(wěn)定性較差，Chien k-ε殘差收斂速度最慢；在1 200～1 900步之間，SST k-ω和k-ε-Rt模型收斂曲線都出現(xiàn)了不同程度的振蕩，計算穩(wěn)定性相對較差，而B-L和S-A模型的計算穩(wěn)定性和收斂性更好。

圖15 Rotor37收斂史Fig.15 Convergence history of the Rotor 37 case

3）變工況性能分析

從總壓比特性對比圖（圖16）中可以看出，在93%～95%阻塞流量之間，B-L模型計算得到的總壓比與實驗數(shù)值最為吻合，而SST k-ω模型的計算結(jié)果與實驗數(shù)值誤差最大為1.6%；在95%～100%阻塞流量之間，k-ε-Rt模型的計算結(jié)果與實驗數(shù)值最為吻合，而在接近阻塞流量處Chien k-ε和SST k-ω模型的計算結(jié)果與實驗數(shù)值相差較大，最大數(shù)值誤差為2.4%。

圖16 總壓比特性對比Fig.16 Comparison of total pressure ratio characteristics

從絕熱效率特性對比圖（圖17）中可以看出，在92%～95%阻塞流量之間，S-A模型計算得到絕熱效率與實驗數(shù)值最為接近，而SST k-ω模型的計算結(jié)果與實驗數(shù)值相比誤差最大為2.07%。在95%～99%之間，S-A和k-ε-Rt模型的計算結(jié)果相較于其它三種模型與實驗數(shù)值更為接近，Chien k-ε模型則過高預測了絕熱效率，與實驗誤差最大為1.5%。同時在靠近阻塞流量處，SST k-ω和Chien k-ε模型的計算結(jié)果與實驗誤差最小，而S-A模型在這里的預測結(jié)果最差，與實驗數(shù)值誤差最大為1.56%。

圖17 絕熱效率特性對比Fig.17 Comparison of adiabatic efficiency characteristics

4 結(jié)論

本文以Durham壓氣機葉柵和NASA rotor37壓氣機轉(zhuǎn)子為研究對象，將五種不同的湍流模型分別應用于兩個算例的數(shù)值模擬中，對計算得到的性能參數(shù)和流場進行了對比分析，得到以下結(jié)論：

1）通過對計算得到的葉片表面壓力系數(shù)分布、絕熱效率特性曲線、總壓比特性曲線、總溫比、總壓比等模擬結(jié)果和實驗值進行對比分析，Chien k-ε和B-L模型的計算精度最差，k-ε-Rt模型的計算結(jié)果與實驗數(shù)值符合較好，計算精度相對較高，能夠應用于低聲速葉珊的數(shù)值計算。

2）通過對流場的分析對比可以看出，B-L和Chien k-ε模型對流場細節(jié)捕捉不夠準確，而k-ε-Rt和SST kω模型對于激波的預測更為準確，能夠較好地反映真實的流動情況，可以應用于跨聲速軸流壓氣機的數(shù)值模擬，但是靠近尾緣處的氣流分離現(xiàn)象的模擬結(jié)果不夠明顯。

3）從各湍流模型的計算收斂史可以看出，Chien k-ε模型的收斂性最差，k-ε-Rt模型的計算穩(wěn)定性較差，而S-A模型的計算穩(wěn)定性以及殘差收斂性都表現(xiàn)得更好。