航空發(fā)動機(jī)葉輪機(jī)流場數(shù)值模擬方法研究

唐曉毅 余雅琪

（中國航發(fā)湖南動力機(jī)械研究所，湖南株洲 412002）

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

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

1.N-S方程求解

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

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

2.Spalart-Allamaras湍流模型方程

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

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

方程左端第二項(xiàng)分別代表湍流方程無粘通量和粘性通量，P、Ds和Df為方程源項(xiàng)表達(dá)式。由于引入了連續(xù)方程，Df有別于標(biāo)準(zhǔn)表達(dá)式，并且本文將其中的交叉偏導(dǎo)項(xiàng)作為粘性項(xiàng)消除了梯度計(jì)算的困難。

實(shí)際應(yīng)用中，S-A方程偶爾存在計(jì)算發(fā)散的現(xiàn)象，綜合分析發(fā)現(xiàn)主要原因?yàn)樵错?xiàng)量值過大，因此對進(jìn)行限制，防止其出現(xiàn)負(fù)值或小值：

對r進(jìn)行限制，防止其值過大：

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

以增強(qiáng)計(jì)算不穩(wěn)定工況點(diǎn)附近湍流方程求解的穩(wěn)定性。變量定義及公式中常量見文獻(xiàn)[6]中值。

3.壁面邊界條件設(shè)置

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

式（9）中Re2的參考長度為靠近壁面第一層網(wǎng)格的法向距離，參考速度為靠近壁面第一個(gè)網(wǎng)格點(diǎn)處的相對速度Vref。計(jì)算得到摩擦系數(shù)Cf之后，可計(jì)算出壁面處的切應(yīng)力：

4.算例驗(yàn)證

NASA Rotor 37為壓比20核心壓氣機(jī)的進(jìn)口級轉(zhuǎn)子，其基元葉型為多圓弧葉型，采用低展弦比、高稠度設(shè)計(jì)的，其葉片進(jìn)口相對馬赫數(shù)均為超聲速，設(shè)計(jì)壓比2.1，絕熱效率高達(dá)0.877[9]，其內(nèi)部流場結(jié)構(gòu)復(fù)雜，成為眾多學(xué)者的研究對象。本文也以其為對象檢驗(yàn)所開發(fā)求解器的準(zhǔn)確性。

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

圖1 NASA Rotor 37計(jì)算網(wǎng)格

圖2 連續(xù)方程殘差收斂史

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

圖3 NASA Rotor 37計(jì)算與試驗(yàn)壓比特性曲線

圖4 NASA Rotor 37計(jì)算與試驗(yàn)效率特性曲線

表1 NASA Rotor 37計(jì)算與試驗(yàn)性能參數(shù)對比

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

圖5 98%堵塞流量工況轉(zhuǎn)子出口50%葉高葉片通道相對馬赫數(shù)分布

5.結(jié)論

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