張 健，唐志共，鄧有奇，周乃春

(中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，四川綿陽(yáng) 621000)

基于混合平面法的跨聲速軸流壓氣機(jī)三維數(shù)值模擬研究

張 健*，唐志共，鄧有奇，周乃春

(中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，四川綿陽(yáng) 621000)

利用自主開(kāi)發(fā)的三維流場(chǎng)求解器MFlow，基于混合平面法技術(shù)對(duì)航空發(fā)動(dòng)機(jī)壓氣機(jī)流場(chǎng)進(jìn)行了數(shù)值模擬。在轉(zhuǎn)/靜子交界面處的不同的展向位置，將流動(dòng)周向平均后進(jìn)行傳遞，使得非定常流動(dòng)問(wèn)題簡(jiǎn)化成為“準(zhǔn)定?！眴?wèn)題，從而將計(jì)算過(guò)程簡(jiǎn)化為只需要對(duì)一個(gè)葉片流道進(jìn)行計(jì)算，大大減小了計(jì)算量。通過(guò)計(jì)算得到了壓氣機(jī)的總體工作特性和幾個(gè)關(guān)鍵氣動(dòng)參數(shù)的展向分布曲線，并將結(jié)果與試驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比分析。計(jì)算結(jié)果和試驗(yàn)數(shù)據(jù)基本吻合，驗(yàn)證了計(jì)算程序的適用性和可靠性。

壓氣機(jī);數(shù)值模擬;混合平面法;氣動(dòng)性能

0 引言

航空發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)十分復(fù)雜，對(duì)于其可靠性、壽命及可維護(hù)性要求極為苛刻［1］，對(duì)于軍用航空發(fā)動(dòng)機(jī)來(lái)說(shuō)，推重比更是衡量其先進(jìn)程度的一個(gè)重要指標(biāo)。為了提高航空發(fā)動(dòng)機(jī)的性能，普遍的做法是提高渦輪的進(jìn)口溫度、增加壓縮比、提高壓氣機(jī)和渦輪的效率［2］。壓氣機(jī)是航空發(fā)動(dòng)機(jī)的重要部件，壓氣機(jī)的氣動(dòng)性能是影響整個(gè)發(fā)動(dòng)機(jī)的效率、推力和運(yùn)轉(zhuǎn)性能的關(guān)鍵因素，因此充分認(rèn)識(shí)壓氣機(jī)內(nèi)部的流動(dòng)機(jī)理，準(zhǔn)確預(yù)測(cè)其氣動(dòng)性能對(duì)發(fā)動(dòng)機(jī)的設(shè)計(jì)有著重要意義。

隨著高性能計(jì)算機(jī)的快速發(fā)展，計(jì)算流體力學(xué)(CFD)技術(shù)被越來(lái)越多的用于航空發(fā)動(dòng)機(jī)壓氣機(jī)的設(shè)計(jì)和分析，其優(yōu)點(diǎn)是周期短、成本低、效率高，是傳統(tǒng)試驗(yàn)方法很好的互補(bǔ)。壓氣機(jī)CFD模擬的一大難點(diǎn)是如何模擬轉(zhuǎn)/靜子葉片排之間相對(duì)運(yùn)動(dòng)引起的非定常流動(dòng)。采用滑移面網(wǎng)格［3］進(jìn)行全流道非定常計(jì)算，上下游葉片排之間利用插值進(jìn)行數(shù)據(jù)傳遞，能夠捕捉到這些非定?，F(xiàn)象，然而缺點(diǎn)是計(jì)算量太大。Denton提出的混合平面方法［4］基本思想很簡(jiǎn)單，轉(zhuǎn)/靜界面上下游的通量在相同展向高度處先通過(guò)周向平均后再進(jìn)行傳遞，將非定常流動(dòng)簡(jiǎn)化為準(zhǔn)定常問(wèn)題，從而只需要對(duì)一個(gè)葉片流道進(jìn)行計(jì)算，大大減小了計(jì)算量。雖然這種方法忽略了轉(zhuǎn)/靜子相對(duì)運(yùn)動(dòng)造成的非定常效應(yīng)，但是依然能夠準(zhǔn)確模擬壓氣機(jī)的性能和一些流動(dòng)特性。因此，本文選取混合平面法，在已有的三維流場(chǎng)求解器MFlow的基礎(chǔ)上，著力于自主開(kāi)發(fā)壓氣機(jī)內(nèi)流模擬程序，達(dá)到對(duì)壓氣機(jī)氣動(dòng)性能的預(yù)測(cè)能力。相比于商業(yè)軟件，利用自主開(kāi)發(fā)的程序具有更好的靈活性，能夠應(yīng)對(duì)各種不同的復(fù)雜問(wèn)題。本文通過(guò)對(duì)NASA的兩個(gè)不同壓氣機(jī)Rotor 67［5］和Stage 35［6］進(jìn)行計(jì)算，并同試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析，檢驗(yàn)了軟件對(duì)壓氣機(jī)性能的模擬能力，能夠?yàn)閴簹鈾C(jī)氣動(dòng)設(shè)計(jì)提供參考。

1 計(jì)算模型

1.1 控制方程

壓氣機(jī)內(nèi)部流場(chǎng)的三維可壓縮粘性 Navier-Stokes(N-S)方程在旋轉(zhuǎn)直角坐標(biāo)系(x，y，z)下的形式為:

式中，Ω表示控制體的體積，S表示控制體封閉面的面積，U為守恒變量，F(xiàn)I為無(wú)粘通量，F(xiàn)V為粘性通量，ST為科里奧利力和離心力引起的源項(xiàng)，具體定義參考文獻(xiàn)［7］。

1.2 數(shù)值方法

對(duì)于式(1)的雷諾平均N-S方程采用基于控制體的有限體積方法離散。采用Roe迎風(fēng)格式［8］計(jì)算通過(guò)控制體單元面的無(wú)粘通量。利用格林-高斯公式計(jì)算控制體單元的梯度，選擇 Venkataknshnan限制器［9］將控制體的原始變量通過(guò)線性重構(gòu)插值到控制體單元面左右，達(dá)到更高的空間精度。本文中粘性項(xiàng)采用中心格式計(jì)算，時(shí)間項(xiàng)利用LU-SGS隱式推進(jìn)方法［10］，并采用局部時(shí)間步長(zhǎng)加速收斂的技術(shù)，提高計(jì)算效率。湍流模型選取Spalart-Allmaras(SA)一方程模型［11］。

1.3 邊界條件

進(jìn)出口采用特征邊界條件。在進(jìn)口，給定總溫T0，總壓P0以及來(lái)流速度方向。向進(jìn)口外發(fā)出的黎曼不變量R-由第一個(gè)內(nèi)部邊界點(diǎn)外插計(jì)算得到，利用R-計(jì)算出絕對(duì)速度大小V，再根據(jù)速度方向?qū)分解為各個(gè)速度分量［12］:

式中，vd代表邊界內(nèi)部第一個(gè)單元速度矢量，cd為當(dāng)?shù)芈曀?，γ為比熱比，n為法向量，Cp為定壓比熱。密度ρ和壓力p利用T0、P0、V之間的代數(shù)關(guān)系以及等熵關(guān)系計(jì)算得到。

在出口邊界，展向高度r處的靜壓是給定的。虛擬網(wǎng)格的密度和速度分量直接通過(guò)外插得到，而當(dāng)?shù)仂o壓結(jié)合式(4)給出的徑向平衡方程和給定位置的靜壓積分得到。

由于利用單個(gè)葉片通道的網(wǎng)格進(jìn)行計(jì)算，因此引入了周期邊界條件。對(duì)于兩個(gè)周期面上的網(wǎng)格，首先需要根據(jù)坐標(biāo)建立一一對(duì)應(yīng)的關(guān)系，在計(jì)算時(shí)周期邊界上的網(wǎng)格類似內(nèi)部網(wǎng)格，虛擬網(wǎng)格物理量直接通過(guò)相對(duì)應(yīng)周期邊界的網(wǎng)格得到。

對(duì)于粘性物面，采用無(wú)滑移邊界條件，即相對(duì)于物面的速度為零。物面分轉(zhuǎn)動(dòng)部分和靜止部分，對(duì)于轉(zhuǎn)子和部分輪轂，有固定的轉(zhuǎn)動(dòng)速度。物面為絕熱壁，物面溫度直接由內(nèi)部點(diǎn)外插得到。

1.4 轉(zhuǎn)/靜界面處理

對(duì)于轉(zhuǎn)/靜界面，也是一種特殊的邊界條件。本文對(duì)轉(zhuǎn)/靜界面的處理采用混合平面法，如圖1所示。在混合面處，對(duì)上下游變量沿圓周方向(θ方向)分別進(jìn)行周向平均，將平均量作為邊界條件的值進(jìn)行傳遞，在每一步的計(jì)算迭代中，上游混合面網(wǎng)格的平均值傳遞給下游虛擬網(wǎng)格，同理下游混合面網(wǎng)格的平均值傳遞給上游虛擬網(wǎng)格。最終收斂時(shí)，混合面上下游的平均值將趨向一致，達(dá)到質(zhì)量、動(dòng)量和能量守恒的目的。

對(duì)于平均方法的選取，有多種不同的選擇，包括面積平均、質(zhì)量平均等等。本文采用通量平均策略(通常稱為Mixed-out Average)［13］，以保證質(zhì)量、動(dòng)量和能量的守恒。為了求解平均量，首先計(jì)算出無(wú)粘通量的面積平均，再根據(jù)式(5)的關(guān)系式反解出構(gòu)成這些通量的原始變量。

圖1 混合平面法示意圖Fig.1 Mixing plane method

式中，I代表通量，A代表面積，ε代表單位法向分量，x、θ、r分別代表軸向方向、圓周方向和徑向方向，U為速度在法向的投影，e為能量，上標(biāo)的橫線代表是平均量。

式(5)給出了一個(gè)關(guān)于p珋的二次方程，其解為:

式中，正號(hào)用于亞聲速流動(dòng)。進(jìn)而通過(guò)式(7)迅速得到其他幾個(gè)平均量:

2 算例與分析

2.1 NASA rotor 67

NASA rotor 67是一個(gè)跨聲速壓氣機(jī)轉(zhuǎn)子葉片，其在33.25 kg/s的流量下設(shè)計(jì)壓比為1.63，一圈共有22個(gè)葉片，設(shè)計(jì)轉(zhuǎn)速為16 043 r/min，葉尖轉(zhuǎn)速達(dá)到429 m/s，葉尖的相對(duì)馬赫數(shù)為1.38。Strazisar和Wood等對(duì)這個(gè)轉(zhuǎn)子進(jìn)行過(guò)測(cè)試，得到了詳盡的試驗(yàn)數(shù)據(jù)［5］，Chima等也對(duì)rotor 67進(jìn)行過(guò)數(shù)值計(jì)算［14］。本文首先利用MFlow對(duì)這個(gè)單獨(dú)的葉片進(jìn)行計(jì)算，并與試驗(yàn)數(shù)據(jù)對(duì)比，對(duì)計(jì)算程序的控制方程、邊界條件、湍流模型等進(jìn)行了驗(yàn)證。計(jì)算網(wǎng)格采用全六面體網(wǎng)格，忽略了葉尖間隙，總網(wǎng)格量在70萬(wàn)左右，如圖2所示。

圖2 Rotor 67葉片和計(jì)算網(wǎng)格Fig.2 Rotor 67 blade and computational grid

通過(guò)調(diào)節(jié)出口背壓來(lái)計(jì)算得到不同流量下總壓比和絕熱效率的工況曲線，并且同試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，如圖3所示，其中對(duì)流量除以壅塞流量進(jìn)行了歸一化。試驗(yàn)測(cè)得的壅塞流量為34.96 kg/s，而計(jì)算得到的壅塞流量為34.68 kg/s，之間相差0.8%。計(jì)算得到的失速流量在最大流量的92%左右，與試驗(yàn)相符。計(jì)算得到的總壓比和絕熱效率變化趨勢(shì)和試驗(yàn)吻合較好，除了在近失速狀態(tài)時(shí)計(jì)算的壓比要偏低一些。

圖3 Rotor 67總體性能隨流量變化Fig.3 Overall performance of rotor 67

圖4 給出了最佳效率點(diǎn)附近工況從葉根到葉尖30%和70%兩個(gè)不同展向位置計(jì)算和試驗(yàn)對(duì)比的相對(duì)馬赫數(shù)等值線圖。在70%位置，葉片前緣處有一道強(qiáng)的弓形激波，在靠近尾緣附近的流道間有一道強(qiáng)的正激波，激波位置和試驗(yàn)結(jié)果一致。在30%位置，葉片吸力面前端形成一塊小的超聲速區(qū)域，和試驗(yàn)結(jié)果相似。計(jì)算結(jié)果驗(yàn)證了本文程序的準(zhǔn)確性和可靠性。

圖4 不同展向位置相對(duì)馬赫數(shù)分布計(jì)算與試驗(yàn)對(duì)比Fig.4 Blade to blade Mach number countors at different span，CFD(left)and experiment(right)

2.2 NASA stage 35

NASA stage 35是多級(jí)跨聲速壓氣機(jī)中的其中一級(jí)，由轉(zhuǎn)子部分rotor 35和靜子部分stator 35組成，其中轉(zhuǎn)子一圈共36個(gè)葉片，靜子46個(gè)葉片。Reid和Moore曾經(jīng)對(duì)stage 35的氣動(dòng)性能進(jìn)行了測(cè)試并給出了詳細(xì)的數(shù)據(jù)［6］，這些數(shù)據(jù)為眾多CFD軟件模擬多級(jí)壓氣機(jī)提供了參考。通過(guò)對(duì)stage 35的計(jì)算，進(jìn)一步驗(yàn)證了程序?qū)Α盎旌掀矫娣ā睉?yīng)用的正確性。stage 35的計(jì)算網(wǎng)格如圖5所示，轉(zhuǎn)/靜子網(wǎng)格塊只取一個(gè)流道，各自生成，同樣周期面網(wǎng)格完全對(duì)應(yīng)，在轉(zhuǎn)子尾緣和靜子前緣中間位置上下游網(wǎng)格共面，即“混合平面”。

圖5 Stage 35計(jì)算網(wǎng)格示意圖Fig.5 Stage 35 computational grid

圖6 給出了計(jì)算得到的在葉高中間位置處壓力的等值云圖，可以看到在轉(zhuǎn)子前緣處有一道強(qiáng)的弓形激波，在靜子葉片15%弦長(zhǎng)位置有一道正激波。在上下游的交界面處，當(dāng)?shù)仂o壓并不連續(xù)，這是因?yàn)椤盎旌掀矫妗狈ㄊ窃诖死弥芟蚱骄颠M(jìn)行傳遞的，僅僅保證了上下游的平均量的連續(xù)性。圖7給出了近失速狀態(tài)對(duì)熵值進(jìn)行周向平均后在子午面上的分布云圖，可以看到在轉(zhuǎn)子葉尖縫隙之后有明顯的熵增，這說(shuō)明葉尖泄流是造成失速工況下流動(dòng)損失的很大原因。

圖8給出了stage 35總壓比、總溫比和絕熱效率計(jì)算結(jié)果，并同試驗(yàn)數(shù)據(jù)以及NASA Glenn研究中心的求解器SWIFT［15］的計(jì)算結(jié)果進(jìn)行了對(duì)比?？傮w來(lái)說(shuō)，CFD計(jì)算的結(jié)果和試驗(yàn)比較一致，對(duì)比MFlow和SWIFT的結(jié)果，MFlow在某些流量工況下壓比和溫比數(shù)值要偏低一些。

圖6 Stage 35一半葉高位置靜壓計(jì)算結(jié)果Fig.6 Static pressure at mid-span of stage 35

圖7 周向平均熵值云圖Fig.7 Pitch averaged entropy contours

圖8 Stage 35整體性能計(jì)算結(jié)果和試驗(yàn)對(duì)比Fig.8 Overall performance of stage 35，pressure ratio，temperature ratio and effiency

圖9 給出了最大效率狀態(tài)和近失速狀態(tài)rotor 35的性能展向分布，計(jì)算結(jié)果和試驗(yàn)比較一致。可以看出在40%展向位置以下有一個(gè)低的總壓區(qū)域，即總壓變化有呈“虧損”趨勢(shì)，造成這個(gè)現(xiàn)象的原因是“角區(qū)失速”。本文的計(jì)算結(jié)果同樣捕捉到了這個(gè)現(xiàn)象。但是總壓并沒(méi)有呈現(xiàn)出試驗(yàn)?zāi)菢悠偷梅浅Ｃ黠@，分析其原因很可能是因?yàn)橛?jì)算時(shí)忽略了靜子的輪轂和轉(zhuǎn)子之間的間隙泄流導(dǎo)致的。

圖9 Stage 35最佳效率(PE，Peak Efficiency)狀態(tài)和近失速(NS，Near Stall)狀態(tài)出口展向性能分布Fig.9 Spanwise profiles of exit at peak effiency(PE，Peak Efficiency)state and near stall(NS，Near Stall)state

3 結(jié)論

本文采用了基于周向平均的混合面法，對(duì)航空發(fā)動(dòng)機(jī)壓氣機(jī)轉(zhuǎn)/靜子相對(duì)運(yùn)動(dòng)進(jìn)行了模擬，使得復(fù)雜的非定常計(jì)算簡(jiǎn)化為準(zhǔn)定常計(jì)算，大大減小了計(jì)算量。雖然由于方法的局限性導(dǎo)致在混合面處無(wú)法保證物理量的連續(xù)性，但是不影響一些關(guān)鍵氣動(dòng)性能參數(shù)的計(jì)算，具有較好的工程應(yīng)用價(jià)值。通過(guò)和試驗(yàn)數(shù)據(jù)的對(duì)比，驗(yàn)證了本文程序的適用性和可靠性，可以為壓氣機(jī)設(shè)計(jì)進(jìn)行氣動(dòng)性能預(yù)測(cè)。但是，在某些方面模擬精度仍存在不足，例如在靠近輪轂附近的總壓虧損，計(jì)算得到的總壓沿展向分布的曲線在這一區(qū)域雖然也有減弱的趨勢(shì)，但遠(yuǎn)沒(méi)有試驗(yàn)結(jié)果明顯。因此，計(jì)算格式、湍流模型、網(wǎng)格等因素對(duì)計(jì)算結(jié)果的影響還需要在今后工作中進(jìn)一步研究。

［1］Zhou Sheng，Han Zhenxue，Meng Qingguo.Some applications of computational fluid dynamics in the field of air-breathing propulsion［J］.Acta Aerodynamica Sinica，1998，16(1):97-107.(in Chinese)周盛，韓振學(xué)，孟慶國(guó).計(jì)算流體力學(xué)在吸氣式推進(jìn)領(lǐng)域中的某些應(yīng)用［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1998，16(1):97-107.

［2］Li Jianfeng，Lyu Junfu.Performance study of technology of combustion inside turbine using in fanjet［J］.Acta Aerodynamica Sinica，2010，28(3):356-362.(in Chinese)李建鋒，呂俊復(fù).渦輪燃燒技術(shù)在渦扇航空發(fā)動(dòng)機(jī)上的應(yīng)用分析［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28(3):356-362.

［3］Rai M.Three-dimensional Navier-Stokes simulations of turbine rotor-stator interaction.I-Methodology［J］.AIAA Journal of Propulsion and Power，1989，5(3):305-311.

［4］Denton J D.The calculation of 3-D viscous flow through multistage turbomachines［J］.Journal of Turbomachinery，1992，114:18-26.

［5］Strazisar A J，Wood J R，Hathaway M D，et al.Laser anemometer measurements in a transonic axial flow fan rotor［R］.NASA TP-2879，1989.

［6］Reid L，Moore R D.Performance of a single-stage axial-flow transonic compressor with rotor and stator aspect ratios of 1.19 and 1.26，respectively，and with design pressure ratio of 1.82［R］.NASA TP-1338，1978.

［7］Chen J P，Ghosh A R，Sreenivas K，et al.Comparison of computations using Navier-Stokes equations in rotating and fixed coordinates for flow through turbo machinery［R］.AIAA-97-0878，1997.

［8］Luo H，Baum J D，Lohner R.An improved finite volume scheme for compressible flows on unstructured grids［R］.AIAA-95-0348，1995.

［9］Venkatakrishnan V.On the accuracy of limiters and convergence to steady state solutions［R］.AIAA-93-0880，1993.

［10］Rieger H，Jameson A.Solution of steady 3-D compressible Euler and Navier-Stokes equations by an implicit L U scheme［R］.AIAA-88-0619，1988.

［11］Spalart P R，Allmaras S R.A one equation turbulence model for aerodynamic flows［R］.AIAA-92-0439，1992.

［12］Chima R V.Viscous Three-dimensional calculations of transonic fan performance［R］.NASA TM-103800，1991.

［13］Chima R V.Calculation of multistage turbomachinery using steady characteristic boundary conditions［R］.NASA TM-206613，1998.

［14］Chima R V.Viscous three-dimensional calculations of transonic fan performance［R］.NASA TM-103800，1991.

［15］Chima R V.SWIFT code assessment for two similar transonic compressors［R］.NASA TM-2009-215520，2009.

Three-dimensional numerical simulations for transonic axial compressors using a mixing plane method

Zhang Jian*，Tang Zhigong，Deng Youqi，Zhou Naichun
(Computational Aerodynamics Institute of China Aerodynamics Research and Development Center，Mianyang 621000，China)

An in-house three-dimensional unstructured-grid flow solver based on a mixing plane method was used to simulate the flow field of the compressors.Pitchwise averaging was used at the rotor/ stator interface for information transport，making the unsteady problem to become a“qusi-steady”one.Thereafter the computational gird needs only one flow passage，which greatly reduced the amount of the calculation.Two different transonic compressors，NASA rotor 67 and stage 35，were tested in this paper.The overall performance maps and spanwise profiles of several aerodyamic parameters were obtained from the calculation，and the results were compared to the experiment data.The results presented here including pressure ratio，temperature ratio and adiabatic efficiency are in good agreement with the experiment data，indicating the applicability and reliability of the method.

turbomachine;numerical simulation;mixing plane;aerodynamic performance

V211.3

:Adoi:10.7638/kqdlxxb-2014.0019

0258-1825(2015)05-0631-05

2014-03-25;

:2014-06-16

張健*(1989-)，男，湖北鐘祥人，碩士研究生，研究方向:計(jì)算空氣動(dòng)力學(xué).E-mail:lsszhangjian@gmail.com

張健，唐志共，鄧有奇，等.基于混合平面法的跨聲速軸流壓氣機(jī)三維數(shù)值模擬研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33(5):631-635.

10.7638/kqdlxxb-2014.0019 Zhang J，Tang Z G，Deng Y Q，et al.Three-dimensional numerical simulations for transonic axial compressors using a mixing plane method［J］.Acta Aerodynamica Sinica，2015，33(5):631-635.