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

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

(中國空氣動力研究與發展中心計算空氣動力研究所，四川綿陽 621000)

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

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

(中國空氣動力研究與發展中心計算空氣動力研究所，四川綿陽 621000)

利用自主開發的三維流場求解器MFlow，基于混合平面法技術對航空發動機壓氣機流場進行了數值模擬。在轉/靜子交界面處的不同的展向位置，將流動周向平均后進行傳遞，使得非定常流動問題簡化成為“準定常”問題，從而將計算過程簡化為只需要對一個葉片流道進行計算，大大減小了計算量。通過計算得到了壓氣機的總體工作特性和幾個關鍵氣動參數的展向分布曲線，并將結果與試驗數據進行了對比分析。計算結果和試驗數據基本吻合，驗證了計算程序的適用性和可靠性。

壓氣機;數值模擬;混合平面法;氣動性能

0 引言

航空發動機內部流動十分復雜，對于其可靠性、壽命及可維護性要求極為苛刻［1］，對于軍用航空發動機來說，推重比更是衡量其先進程度的一個重要指標。為了提高航空發動機的性能，普遍的做法是提高渦輪的進口溫度、增加壓縮比、提高壓氣機和渦輪的效率［2］。壓氣機是航空發動機的重要部件，壓氣機的氣動性能是影響整個發動機的效率、推力和運轉性能的關鍵因素，因此充分認識壓氣機內部的流動機理，準確預測其氣動性能對發動機的設計有著重要意義。

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

1 計算模型

1.1 控制方程

壓氣機內部流場的三維可壓縮粘性 Navier-Stokes(N-S)方程在旋轉直角坐標系(x，y，z)下的形式為:

式中，Ω表示控制體的體積，S表示控制體封閉面的面積，U為守恒變量，FI為無粘通量，FV為粘性通量，ST為科里奧利力和離心力引起的源項，具體定義參考文獻［7］。

1.2 數值方法

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

1.3 邊界條件

進出口采用特征邊界條件。在進口，給定總溫T0，總壓P0以及來流速度方向。向進口外發出的黎曼不變量R-由第一個內部邊界點外插計算得到，利用R-計算出絕對速度大小V，再根據速度方向將V分解為各個速度分量［12］:

式中，vd代表邊界內部第一個單元速度矢量，cd為當地聲速，γ為比熱比，n為法向量，Cp為定壓比熱。密度ρ和壓力p利用T0、P0、V之間的代數關系以及等熵關系計算得到。

在出口邊界，展向高度r處的靜壓是給定的。虛擬網格的密度和速度分量直接通過外插得到，而當地靜壓結合式(4)給出的徑向平衡方程和給定位置的靜壓積分得到。

由于利用單個葉片通道的網格進行計算，因此引入了周期邊界條件。對于兩個周期面上的網格，首先需要根據坐標建立一一對應的關系，在計算時周期邊界上的網格類似內部網格，虛擬網格物理量直接通過相對應周期邊界的網格得到。

對于粘性物面，采用無滑移邊界條件，即相對于物面的速度為零。物面分轉動部分和靜止部分，對于轉子和部分輪轂，有固定的轉動速度。物面為絕熱壁，物面溫度直接由內部點外插得到。

1.4 轉/靜界面處理

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

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

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

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

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

式中，正號用于亞聲速流動。進而通過式(7)迅速得到其他幾個平均量:

2 算例與分析

2.1 NASA rotor 67

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

圖2 Rotor 67葉片和計算網格Fig.2 Rotor 67 blade and computational grid

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

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

圖4 給出了最佳效率點附近工況從葉根到葉尖30%和70%兩個不同展向位置計算和試驗對比的相對馬赫數等值線圖。在70%位置，葉片前緣處有一道強的弓形激波，在靠近尾緣附近的流道間有一道強的正激波，激波位置和試驗結果一致。在30%位置，葉片吸力面前端形成一塊小的超聲速區域，和試驗結果相似。計算結果驗證了本文程序的準確性和可靠性。

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

2.2 NASA stage 35

NASA stage 35是多級跨聲速壓氣機中的其中一級，由轉子部分rotor 35和靜子部分stator 35組成，其中轉子一圈共36個葉片，靜子46個葉片。Reid和Moore曾經對stage 35的氣動性能進行了測試并給出了詳細的數據［6］，這些數據為眾多CFD軟件模擬多級壓氣機提供了參考。通過對stage 35的計算，進一步驗證了程序對“混合平面法”應用的正確性。stage 35的計算網格如圖5所示，轉/靜子網格塊只取一個流道，各自生成，同樣周期面網格完全對應，在轉子尾緣和靜子前緣中間位置上下游網格共面，即“混合平面”。

圖5 Stage 35計算網格示意圖Fig.5 Stage 35 computational grid

圖6 給出了計算得到的在葉高中間位置處壓力的等值云圖，可以看到在轉子前緣處有一道強的弓形激波，在靜子葉片15%弦長位置有一道正激波。在上下游的交界面處，當地靜壓并不連續，這是因為“混合平面”法是在此利用周向平均值進行傳遞的，僅僅保證了上下游的平均量的連續性。圖7給出了近失速狀態對熵值進行周向平均后在子午面上的分布云圖，可以看到在轉子葉尖縫隙之后有明顯的熵增，這說明葉尖泄流是造成失速工況下流動損失的很大原因。

圖8給出了stage 35總壓比、總溫比和絕熱效率計算結果，并同試驗數據以及NASA Glenn研究中心的求解器SWIFT［15］的計算結果進行了對比。總體來說，CFD計算的結果和試驗比較一致，對比MFlow和SWIFT的結果，MFlow在某些流量工況下壓比和溫比數值要偏低一些。

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

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

圖8 Stage 35整體性能計算結果和試驗對比Fig.8 Overall performance of stage 35，pressure ratio，temperature ratio and effiency

圖9 給出了最大效率狀態和近失速狀態rotor 35的性能展向分布，計算結果和試驗比較一致。可以看出在40%展向位置以下有一個低的總壓區域，即總壓變化有呈“虧損”趨勢，造成這個現象的原因是“角區失速”。本文的計算結果同樣捕捉到了這個現象。但是總壓并沒有呈現出試驗那樣偏低得非常明顯，分析其原因很可能是因為計算時忽略了靜子的輪轂和轉子之間的間隙泄流導致的。

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

3 結論

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

［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)周盛，韓振學，孟慶國.計算流體力學在吸氣式推進領域中的某些應用［J］.空氣動力學學報，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)李建鋒，呂俊復.渦輪燃燒技術在渦扇航空發動機上的應用分析［J］.空氣動力學學報，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-)，男，湖北鐘祥人，碩士研究生，研究方向:計算空氣動力學.E-mail:lsszhangjian@gmail.com

張健，唐志共，鄧有奇，等.基于混合平面法的跨聲速軸流壓氣機三維數值模擬研究［J］.空氣動力學學報，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.