航空發動機葉輪機流場數值模擬方法研究

唐曉毅 余雅琪

（中國航發湖南動力機械研究所，湖南株洲 412002）

隨著航空發動機朝著高性能方向發展，對葉輪機設計提出了更高的要求，因而采用CFD技術結合優化設計方法成為重要的設計手段。由于優化設計方法需要計算大量樣本以進行尋優或建立代理模型[1-3]，CFD求解器的計算時間成為設計時間成本中的主要部分。計算準確度高、收斂速度快、魯棒性好的CFD求解器可極大的縮短設計周期，節約成本。Jameson等人提出的JST格式在全三維數值仿真中得到廣泛應用，該格式采用一階后差格式離散N-S方程時間項，中心差分離散N-S方程空間項，加入人工粘性保證計算的穩定性。該格式為二階精度，計算量較小、收斂速度快、魯棒性好，自格式推出后不斷得到改進，在葉輪機械CFD領域得到了廣泛的應用。

本文基于改進的JST格式[4-6]，自主開發了葉輪機數值模擬求解器，為滿足高精度流場快速仿真需求，對Spalart-Allamaras湍流模型方程及滑移壁面函數的應用進行研究，開展了高負荷壓氣機轉子葉片流場模擬，通過與試驗流場、數據對比，驗證了計算方法的正確性。

1.N-S方程求解

本文求解器采用有限體積法求解N-S方程，坐標系其絕對圓柱坐標系（z，θ，r）。其中，考慮轉速的積分形式守恒型N-S方程為：

其中Ω為網格單元控制體，Az,Aθ,Ar分別為控制體表面在相應坐標軸方向的投影面積。守恒變量矢量為Q，ω為旋轉角速度。矢量（E、F和G）、（Ev、Fv和Gv），對應（z，θ，r）方向的對流通量及粘性通量，源項為H。采用基于JST格式的五步混合Rung-Kutta顯示時間推進法求解方程，應用了當地時間步長、隱式殘差光順技術加速收斂[5]。

2.Spalart-Allamaras湍流模型方程

S-A模型由Spalart和Allamaras[6-7]于1992年提出，通過求解中間變量的控制方程計算湍流粘性系數，方程的封閉系數經過豐富的實驗數據調校，擁有較高的準確度，相對于零方程湍流模型（如B-L模型），S-A模型計算得到湍流粘性系數在整個流場內保持連續，對存在分離現象的流場也適用，對網格量要求低，方程收斂性好，本文采用S-A模型計算湍流粘性，在保證湍流粘性計算精度的同時獲得較高的計算效率。

同N-S方程（1）式一致，本文在文獻[6]中S-A方程通過引入連續方程，得到了積分形式的守恒型S-A模型方程：

方程左端第二項分別代表湍流方程無粘通量和粘性通量，P、Ds和Df為方程源項表達式。由于引入了連續方程，Df有別于標準表達式，并且本文將其中的交叉偏導項作為粘性項消除了梯度計算的困難。

實際應用中，S-A方程偶爾存在計算發散的現象，綜合分析發現主要原因為源項量值過大，因此對進行限制，防止其出現負值或小值：

對r進行限制，防止其值過大：

壁面處給定為=0，在方程求解過程中限制μt≤5000μl。為增強S-A方程計算的穩定性，對源項（4）、（5）式采用相同格式進行了隱式處理，如（4）式處理方法為：

以增強計算不穩定工況點附近湍流方程求解的穩定性。變量定義及公式中常量見文獻[6]中值。

3.壁面邊界條件設置

對于無滑移壁面邊界條件，如需準確模擬壁面附近粘性底層的流動，需在靠近壁面處布置大量網格，使壁面附近的無量綱距離Y+<10。為了減少近壁面處對網格量及網格尺寸的要求，本文利用Denton[8]提出的壁面函數對壁面附近粘性底層流動進行模擬，假設近壁面第一層網格處在層流底層或者湍流對數律區外，壁面處流體摩擦系數可通過式（9）求出：

式（9）中Re2的參考長度為靠近壁面第一層網格的法向距離，參考速度為靠近壁面第一個網格點處的相對速度Vref。計算得到摩擦系數Cf之后，可計算出壁面處的切應力：

4.算例驗證

NASA Rotor 37為壓比20核心壓氣機的進口級轉子，其基元葉型為多圓弧葉型，采用低展弦比、高稠度設計的，其葉片進口相對馬赫數均為超聲速，設計壓比2.1，絕熱效率高達0.877[9]，其內部流場結構復雜，成為眾多學者的研究對象。本文也以其為對象檢驗所開發求解器的準確性。

在經過網格無關性驗證后，最終選定網格數為：119×39×45（流向×周向×展向），總體網格單元數約20萬，圖3為本文采用的Rotor 37的計算網格，在其葉頂間隙設置了3層網格，葉片通道內共55個網格點，在葉片前、后緣處加密網格以增強對前、后緣圓弧的模擬，如圖1所示。經過約1300個時間步后連續方程殘差相對值下降了4個數量級，進出口流量平衡，達到收斂標準，殘差收斂史見圖2。

圖1 NASA Rotor 37計算網格

圖2 連續方程殘差收斂史

數值計算得到的特性曲線如圖3、圖4所示。從圖中可以看出，特性曲線的變化趨勢與試驗值一致。在數值上，壓比特性計算值略高于對應試驗值，而效率特性預測值總體偏低，偏差約為2%，顯示出較高的精度。表1給出了計算與試驗性能參數對比，可見計算精度不低于3%。

圖3 NASA Rotor 37計算與試驗壓比特性曲線

圖4 NASA Rotor 37計算與試驗效率特性曲線

表1 NASA Rotor 37計算與試驗性能參數對比

S-A湍流模型能夠有效模擬存在分離的流動，并且能得出滿意的湍流粘性分布和模擬尾跡，為驗證所開發求解器對尾跡模擬的可靠性，圖5為Rotor 37轉子出口測量站50%葉高相對馬赫數分布，由圖中可以看出，計算結果在馬赫數數值大小、尾跡大小和尾跡寬度等流動細節方面均與試驗測量值一致，進一步表明了本文所開發的求解器能夠準確模擬分離以及尾跡區的復雜流動，并有較高精度。

圖5 98%堵塞流量工況轉子出口50%葉高葉片通道相對馬赫數分布

5.結論

本文給出了圓柱坐標系下采用JST格式的S-A湍流方程求解數值方法，以及滑移壁面函數的應用。對NASA Rotor 37轉子開展了數值模擬，通過與試驗數據對比，計算結果顯示求解器對高負荷葉輪機特性預測精度較高，較好的捕捉了葉片尾跡，驗證了本文數值程序的正確性和精度，表明本文所采用的方法對跨聲速轉子流場有較高的預測精度，能夠滿足工程精度要求。